![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
富山県のホタルイカの身投げは、産卵行動や環境要因と密接に関連する現象であり、 その量の分析と予測は、生態学的理解の深化や気候変動の指標、さらに地域観光や漁 業資源管理における重要な意義を持っている。
1 過去のホタルイカ身投げデータと、その日の天気、気温、月齢、風などのデータを収集
2 身投げが起きやすい条件の分析、学習
3 ホタルイカの身投げ量の予測モデルの作成
4 身投げ量の予報を提供するWebサイトの作成
1. ホタルイカ掲示板というサイト（https://rara.jp/hotaruika-toyama/）の口コミからスクレイピングで口コミ内容を取得
2. その口コミ内容から、ホタルイカの身投げ量を5段階に分類し、0〜4の整数で数値化
3. その値の日毎の平均をその日の身投げ量とし、身投げ量、日付をjsonファイルに保存
2015年〜2025年まで10年分、それぞれ2〜5月の4ヶ月分の約1300日分のデータを収集した。過去の身投げ量以外の以下に記載したデータを収集するときは、この日付の日のデータを収集する。
気象庁の公式サイト（https://www.data.jma.go.jp/risk/obsdl/#）に、過去の天気、風、気温、降水量などのデータがあり、指定した場所、期間を指定すると、過去の天気などを CSV 形式でダウンロードできる.それをjsonデータに変換する。
API経由で潮位、潮の種類などが取得でき、その取得したデータをjsonファイルに保存する。（https://tide736.net/#google_vignette）
月齢は、西暦と日付から計算できる.日付ごとの月齢をjsonファイルにまとめる。
これらのjsonファイルを日付をキーにして一つのjsonファイルにまとめる。
目的変数（予測したい値）：身投げ量（avg_amount）
説明変数（目的変数に影響を与える値）：下記特徴量
データ分割：1220 日分のデータのうち85％を学習用、残りの15％をテスト用として使用
年、月、日、曜日、年の第何週
月齢は 0～29.53 日で周期的に変化するため、sin と cos で周期性を表現。
月齢 sin は「上弦の月」で最大（1）、「下弦の月」で最小（-1）。
月齢 cos は「新月」で最大（1）、満月で最小（-1）。
時間帯別（10-13 時、14-17 時、18-21 時、22-0 時、1-4 時）の平均気温も算出。
時間帯別の降水量合計も算出。降水の有無をバイナリ変数で表現。
潮型（若潮、長潮、小潮、中潮、大潮）を数値変換。
夜間（21 時～5 時）の満潮・干潮高（平均 / 最大 / 最小）。
気象・潮汐変数の 1 日前、2 日前の値。
目的変数（avg_amount）の 1 日前、2 日前、3 日間平均。
import json import pandas as pd import numpy as np from sklearn.model_selection import TimeSeriesSplit, GridSearchCV from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score from sklearn.preprocessing import LabelEncoder, MinMaxScaler import matplotlib.pyplot as plt import seaborn as sns import lightgbm as lgb import holidays
# イベント日リスト（YYYY-MM-DD形式で実際のイベント日を記入）
EVENT_DATES = [
"2023-03-25", "2023-04-10", "2023-04-20"
]
jp_holidays = holidays.country_holidays("JP")
DIRECTION_MAP = {
'北': 0, '北北東': 1, '北東': 2, '東北東': 3, '東': 4, '東南東': 5, '南東': 6, '南南東': 7,
'南': 8, '南南西': 9, '南西': 10, '西南西': 11, '西': 12, '西北西': 13, '北西': 14, '北北西': 15
}
TIDE_TYPE_MAP = {'若潮': 0, '長潮': 1, '小潮': 2, '中潮': 3, '大潮': 4}
def direction_to_radian(direction):
idx = DIRECTION_MAP.get(direction)
if idx is None:
return np.nan
return 2 * np.pi * idx / 16
def mean_wind_direction(directions):
radians = np.array([direction_to_radian(d) for d in directions if d in DIRECTION_MAP])
if len(radians) == 0:
return np.nan
sin_sum = np.sin(radians).mean()
cos_sum = np.cos(radians).mean()
mean_rad = np.arctan2(sin_sum, cos_sum)
if mean_rad < 0:
mean_rad += 2 * np.pi
return int(np.round(mean_rad * 16 / (2 * np.pi))) % 16
def safe_float(x):
try:
return float(x)
except Exception:
return np.nan
def is_wind_night_time(time_str):
try:
hour = int(time_str.split(' ')[1].split(':')[0])
return (20 <= hour <= 23) or (0 <= hour <= 4)
except Exception:
return False
def is_temp_precip_night_time(time_str):
try:
hour = int(time_str.split(' ')[1].split(':')[0])
return (10 <= hour <= 23) or (0 <= hour <= 4)
except Exception:
return False
def is_night_tide_time(time_str):
try:
hour = int(time_str.split(':')[0])
return (21 <= hour <= 23) or (0 <= hour <= 5)
except Exception:
return False
def get_time_period(hour):
if 10 <= hour <= 13:
return '10_13'
elif 14 <= hour <= 17:
return '14_17'
elif 18 <= hour <= 21:
return '18_21'
elif hour in [22, 23, 0]:
return '22_0'
elif 1 <= hour <= 4:
return '1_4'
else:
return None
all_rows = []
# JSONファイルのパスを適切に設定してください
with open('avg_amount_moon_temperature_wind_precipitation_tide_2.json', encoding='utf-8') as f:
data = json.load(f)
for entry in data:
date = entry['date']
avg_amount = safe_float(entry['avg_amount'])
moon_age = safe_float(entry['moon_age'])
try:
date_dt = pd.to_datetime(date, format='%Y/%m/%d')
year = date_dt.year
month = date_dt.month
day = date_dt.day
weekday = date_dt.weekday()
date_str_ymd = date_dt.strftime('%Y-%m-%d')
week_of_year = date_dt.isocalendar().week
week_of_month = (date_dt.day - 1) // 7 + 1
day_of_year = date_dt.dayofyear
except Exception:
year = np.nan
month = np.nan
day = np.nan
weekday = np.nan
date_str_ymd = None
week_of_year = np.nan
week_of_month = np.nan
day_of_year = np.nan
is_weekend = 1 if weekday in [5, 6] else 0
is_holiday = 1 if date_str_ymd in jp_holidays else 0
is_event = 1 if date_str_ymd in EVENT_DATES else 0
temp_periods = {'10_13': [], '14_17': [], '18_21': [], '22_0': [], '1_4': []}
prec_periods = {'10_13': [], '14_17': [], '18_21': [], '22_0': [], '1_4': []}
temps = []
for d in entry.get('temperature_data', []):
temp = safe_float(d.get('average_temperature(℃)'))
time_str = d.get('time', '')
if temp is not None and temp != '' and time_str:
try:
hour = int(time_str.split(' ')[1].split(':')[0])
except Exception:
continue
if is_temp_precip_night_time(time_str):
temps.append(temp)
period = get_time_period(hour)
if period:
temp_periods[period].append(temp)
temperature_mean = np.nanmean(temps) if temps else np.nan
temperature_max = np.nanmax(temps) if temps else np.nan
temperature_min = np.nanmin(temps) if temps else np.nan
temperature_std = np.nanstd(temps) if temps else np.nan
precs = []
for d in entry.get('precipitation_data', []):
prec = safe_float(d.get('precipitation_totals(mm)'))
time_str = d.get('time', '')
if prec is not None and prec != '' and time_str:
try:
hour = int(time_str.split(' ')[1].split(':')[0])
except Exception:
continue
if is_temp_precip_night_time(time_str):
precs.append(prec)
period = get_time_period(hour)
if period:
prec_periods[period].append(prec)
precipitation_sum = np.nansum(precs) if precs else 0
precipitation_binary = 1 if precipitation_sum > 0 else 0
temp_period_features = {}
prec_period_features = {}
for p in temp_periods:
temp_period_features[f'temperature_mean_{p}'] = np.nanmean(temp_periods[p]) if temp_periods[p] else np.nan
for p in prec_periods:
prec_period_features[f'precipitation_sum_{p}'] = np.nansum(prec_periods[p]) if prec_periods[p] else 0
wind_speeds = [safe_float(d.get('average_wind_speed_m_s')) for d in entry.get('wind_data', [])
if d.get('average_wind_speed_m_s') not in [None, '']
and is_wind_night_time(d.get('time',''))]
wind_speed_mean = np.nanmean(wind_speeds) if wind_speeds else np.nan
wind_speed_max = np.nanmax(wind_speeds) if wind_speeds else np.nan
wind_speed_min = np.nanmin(wind_speeds) if wind_speeds else np.nan
wind_speed_std = np.nanstd(wind_speeds) if wind_speeds else np.nan
wind_dirs = [d.get('wind_direction') for d in entry.get('wind_data', [])
if d.get('wind_direction', '') != '' and is_wind_night_time(d.get('time',''))]
wind_direction_mean = mean_wind_direction(wind_dirs) if wind_dirs else np.nan
tide_type_str = entry.get('tide_data_night', {}).get('tide_type', None)
tide_type_num = TIDE_TYPE_MAP.get(tide_type_str, np.nan)
high_tides = entry.get('tide_data_night', {}).get('high_tide', [])
night_high_tides = [
safe_float(ht.get('tide_height_cm')) for ht in high_tides
if ht.get('time') and is_night_tide_time(ht['time'])
]
high_tide_height_mean = np.nanmean(night_high_tides) if night_high_tides else np.nan
high_tide_height_max = np.nanmax(night_high_tides) if night_high_tides else np.nan
high_tide_height_min = np.nanmin(night_high_tides) if night_high_tides else np.nan
low_tides = entry.get('tide_data_night', {}).get('low_tide', [])
night_low_tides = [
safe_float(lt.get('tide_height_cm')) for lt in low_tides
if lt.get('time') and is_night_tide_time(lt['time'])
]
low_tide_height_mean = np.nanmean(night_low_tides) if night_low_tides else np.nan
low_tide_height_max = np.nanmax(night_low_tides) if night_low_tides else np.nan
low_tide_height_min = np.nanmin(night_low_tides) if night_low_tides else np.nan
all_rows.append({
'date': date,
'year': year,
'month': month,
'day': day,
'avg_amount': avg_amount,
'moon_age': moon_age,
'weekday': weekday,
'week_of_year': week_of_year,
'week_of_month': week_of_month,
'day_of_year': day_of_year,
'is_weekend': is_weekend,
'is_holiday': is_holiday,
'is_event': is_event,
'temperature_mean': temperature_mean,
'temperature_max': temperature_max,
'temperature_min': temperature_min,
'temperature_std': temperature_std,
'precipitation_sum': precipitation_sum,
'precipitation_binary': precipitation_binary,
'wind_speed_mean': wind_speed_mean,
'wind_speed_max': wind_speed_max,
'wind_speed_min': wind_speed_min,
'wind_speed_std': wind_speed_std,
'wind_direction_mean': wind_direction_mean,
'tide_type_num': tide_type_num,
'high_tide_height_mean': high_tide_height_mean,
'high_tide_height_max': high_tide_height_max,
'high_tide_height_min': high_tide_height_min,
'low_tide_height_mean': low_tide_height_mean,
'low_tide_height_max': low_tide_height_max,
'low_tide_height_min': low_tide_height_min,
**temp_period_features,
**prec_period_features
})
df = pd.DataFrame(all_rows) df['date'] = pd.to_datetime(df['date'], format='%Y/%m/%d') df = df[df['date'].dt.month >= 2]
df = df.sort_values('date').reset_index(drop=True)
df['moon_age_sin'] = np.sin(2 * np.pi * df['moon_age'] / 29.53) df['moon_age_cos'] = np.cos(2 * np.pi * df['moon_age'] / 29.53)
df['moon_age'] = df['moon_age'].ffill() df['year'] = df['year'].fillna(df['year'].mode()[0]) df['month'] = df['month'].fillna(df['month'].mode()[0]) df['day'] = df['day'].fillna(df['day'].mode()[0]) df['week_of_year'] = df['week_of_year'].fillna(df['week_of_year'].mode()[0]) df['week_of_month'] = df['week_of_month'].fillna(df['week_of_month'].mode()[0]) df['day_of_year'] = df['day_of_year'].fillna(df['day_of_year'].mode()[0])
for c in ['temperature_mean','temperature_max','temperature_min','temperature_std']:
df[c] = df[c].ffill()
df['precipitation_sum'] = df['precipitation_sum'].fillna(0)
df['precipitation_binary'] = df['precipitation_binary'].fillna(0)
for c in ['wind_speed_mean','wind_speed_max','wind_speed_min','wind_speed_std']:
df[c] = df[c].ffill()
df['wind_direction_mean'] = df['wind_direction_mean'].ffill()
df['tide_type_num'] = df['tide_type_num'].fillna(df['tide_type_num'].mode()[0])
for c in [
'high_tide_height_mean','high_tide_height_max','high_tide_height_min',
'low_tide_height_mean','low_tide_height_max','low_tide_height_min'
]:
df[c] = df[c].fillna(0)
for c in ['weekday','is_weekend','is_holiday','is_event']:
df[c] = df[c].fillna(0)
for p in ['10_13', '14_17', '18_21', '22_0', '1_4']:
df[f'temperature_mean_{p}'] = df[f'temperature_mean_{p}'].ffill()
df[f'precipitation_sum_{p}'] = df[f'precipitation_sum_{p}'].fillna(0)
le = LabelEncoder() df['wind_direction_encoded'] = le.fit_transform(df['wind_direction_mean'])
for win in [3]:
df[f'temperature_mean_{win}d_mean'] = df['temperature_mean'].rolling(window=win, min_periods=1).mean().shift(1)
df[f'precipitation_sum_{win}d_sum'] = df['precipitation_sum'].rolling(window=win, min_periods=1).sum().shift(1)
df['avg_amount_lag1'] = df['avg_amount'].shift(1)
### 変更点 ### 潮関連のデータをラグ特徴量のリストから除外
lag_cols = [
'moon_age', 'moon_age_sin', 'moon_age_cos', 'temperature_mean', 'temperature_max', 'temperature_min',
'temperature_std', 'precipitation_sum', 'precipitation_binary', 'wind_speed_mean', 'wind_speed_max', 'wind_speed_min',
'wind_speed_std', 'wind_direction_encoded',
'day_of_year',
]
for col in lag_cols:
df[f'{col}_lag1'] = df[col].shift(1)
df[f'{col}_lag2'] = df[col].shift(2)
df = df.dropna().reset_index(drop=True)
features = [
'year', 'month',
'week_of_month',
'day_of_year',
'moon_age', 'moon_age_sin', 'moon_age_cos',
'weekday', 'is_weekend', 'is_holiday', 'is_event',
'temperature_mean', 'temperature_max', 'temperature_min', 'temperature_std',
'precipitation_sum', 'precipitation_binary',
'wind_speed_mean', 'wind_speed_max', 'wind_speed_min', 'wind_speed_std',
'wind_direction_encoded',
'tide_type_num',
'high_tide_height_mean','high_tide_height_max','high_tide_height_min',
'low_tide_height_mean','low_tide_height_max','low_tide_height_min',
'avg_amount_lag1',
'temperature_mean_3d_mean',
'precipitation_sum_3d_sum'
]
features += [f'temperature_mean_{p}' for p in ['10_13', '14_17', '18_21', '22_0', '1_4']]
features += [f'precipitation_sum_{p}' for p in ['10_13', '14_17', '18_21', '22_0', '1_4']]
### 変更点 ### lag_colsに対応して特徴量リストを更新
features += [f"{col}_lag1" for col in lag_cols]
features += [f"{col}_lag2" for col in lag_cols]
X = df[features] y = df['avg_amount'].values.reshape(-1, 1)
scaler_X = MinMaxScaler() scaler_y = MinMaxScaler() X_scaled = scaler_X.fit_transform(X) y_scaled = scaler_y.fit_transform(y)
tscv = TimeSeriesSplit(n_splits=5)
param_grid = {
'num_leaves': [15, 31, 63],
'learning_rate': [0.01, 0.05, 0.1],
'n_estimators': [100, 300, 500],
'random_state': [42]
}
lgb_model = lgb.LGBMRegressor()
gs = GridSearchCV(lgb_model, param_grid, scoring='r2', cv=tscv, n_jobs=-1, verbose=1)
gs.fit(X_scaled, y_scaled.ravel())
best_model = gs.best_estimator_
best_params = gs.best_params_
print("Best Params:", best_params)
feature_importance = pd.DataFrame({
'feature': features,
'importance': best_model.feature_importances_
}).sort_values('importance', ascending=False)
print(f"\n重要度0の特徴量: {feature_importance[feature_importance['importance'] == 0]['feature'].tolist()}\n")
### 変更点 ### 重要度が10以上の特徴量を選択
selected_features = feature_importance[feature_importance['importance'] >= 10]['feature'].tolist()
print(f"選択された特徴量の数: {len(selected_features)}")
print(f"選択された特徴量: {selected_features}\n")
# 選択された特徴量に目的変数そのものが含まれていないか確認
if 'avg_amount' in selected_features:
print("警告：選択された特徴量に 'avg_amount' が含まれています。除外します。")
selected_features.remove('avg_amount')
X_sel = df[selected_features] scaler_X_sel = MinMaxScaler() X_sel_scaled = scaler_X_sel.fit_transform(X_sel)
test_size = int(len(df) * 0.15) X_train, X_test = X_sel_scaled[:-test_size], X_sel_scaled[-test_size:] y_train, y_test = y_scaled[:-test_size], y_scaled[-test_size:]
best_model.fit(X_train, y_train.ravel()) y_pred = best_model.predict(X_test)
mse = mean_squared_error(y_test, y_pred)
mae = mean_absolute_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
### 変更点 ### 評価指標のタイトルを更新
print("\n=== モデル評価（正規化後, 重要度10以上の特徴量のみ） ===")
print("平均二乗誤差 (MSE):", mse)
print("平均絶対誤差 (MAE):", mae)
print("決定係数 (R²):", r2)
feature_importance_sel = pd.DataFrame({
'feature': selected_features,
'importance': best_model.feature_importances_
}).sort_values('importance', ascending=False)
### 変更点 ### 特徴量重要度のタイトルを更新
print("\n特徴量の重要度（重要度10以上の特徴量のみ）:\n", feature_importance_sel)
plt.figure(figsize=(10, 12))
sns.barplot(x='importance', y='feature', data=feature_importance_sel)
### 変更点 ### グラフのタイトルを更新
plt.title('Feature Importance (LightGBM, Importance >= 10)')
plt.tight_layout()
plt.show()
plt.figure(figsize=(8, 6))
plt.plot(y_test, label='True (normalized)')
plt.plot(y_pred, label='Pred (normalized)')
plt.legend()
### 変更点 ### グラフのタイトルを更新
plt.title('Prediction vs True (Normalized, Importance >= 10)')
plt.show()
y_test_inv = scaler_y.inverse_transform(y_test.reshape(-1, 1))
y_pred_inv = scaler_y.inverse_transform(y_pred.reshape(-1, 1))
plt.figure(figsize=(8, 6))
plt.plot(y_test_inv, label='True (original scale)')
plt.plot(y_pred_inv, label='Pred (original scale)')
plt.legend()
### 変更点 ### グラフのタイトルを更新
plt.title('Prediction vs True (Original Scale, Importance >= 10)')
plt.show()
from joblib import dump
# すでに学習済みのモデル・スケーラー・特徴量リストがある前提 dump(best_model, 'best_model.joblib') dump(scaler_X_sel, 'scaler_X_sel.joblib') dump(scaler_y, 'scaler_y.joblib')
import pickle
with open('selected_features.pkl', 'wb') as f:
pickle.dump(selected_features, f)
・平均二乗誤差 (MSE): 0.0179
・平均絶対誤差 (MAE): 0.0990
・決定係数 (R²): 0.4191
Webサイトのリポジトリ：https://github.com/yuchi1128/hotaruika-bakuwaki-forecast
ターミナルで
# コードをローカルに持ってくる git clone https://github.com/yuchi1128/hotaruika-bakuwaki-forecast.git # Dockerコンテナの起動（自動でフロントエンド、バックエンド、データベースのサ ーバーが起動されます） docker-compose up -d --build
・ホタルイカの身投げ量を 1 週間後まで予測して表示
・その日の天気や月齢、潮の満ち引きなどの情報を表示
・ホタルイカ掬いでの注意点を掲載
・口コミ機能（good・bad、返信）
・フロントエンド：TypeScript（Next.js）
・バックエンド：Go
・データベース：PostgreSQL
・インフラ：Docker
・ホスティング：GCP（Cloud Run）
予測APIのリポジトリ：https://github.com/yuchi1128/hotaruika-prediction-api
・日付を指定しないで GET リクエストを送るとその日から 1 週間のホタルイカの身投げ量を予測して json 形式で返す
・日付を指定して GET リクエストを送るとその日のみのホタルイカの身投げ量を予測してjson で返す
・Python（FastAPI）
・Docker
・GCP（Cloud Run）
