巨大地震のきっかけは月なのか？->月っぽい

from pathlib import Path
from kaggle.api.kaggle_api_extended import KaggleApi
import numpy as np
import pandas as pd
from plotnine import *

DATASET = 'usgs/earthquake-database'
FILENAME = 'database.csv'
DATA_DIR = Path().cwd()

# download dataset
api = KaggleApi()
api.authenticate()
api.dataset_download_file(DATASET, FILENAME, path=DATA_DIR)

# read data file
earthquake = pd.read_csv(
      DATA_DIR.joinpath(FILENAME + '.zip'),
      parse_dates={'datetime': ['Date', 'Time']},
      encoding='utf-8',
      error_bad_lines=False
)

for idx in [3378,7512,20650]:
    earthquake.at[idx, "datetime"] = earthquake.at[idx, "datetime"].split(" ")[0]
earthquake["datetime"] = pd.to_datetime(earthquake["datetime"], utc=True)
earthquake = earthquake.set_index(["datetime"])
earthquake = earthquake.query('Depth < 80')

#!pip install jplephem astropy
from astropy.coordinates import EarthLocation, get_body, AltAz, Longitude, Latitude, Angle
from astropy.time import Time
from astropy.coordinates import solar_system_ephemeris
solar_system_ephemeris.set('de432s')

# 各地震の震央
pos = EarthLocation(Longitude(earthquake["Longitude"], unit="deg"), Latitude(earthquake["Latitude"], unit="deg"))

# 位置を求めたい時刻一覧（地震発生時刻のリスト）
dts = Time(earthquake.index)
# 月の位置（赤道座標から地平座標へ変換）
Mpos = get_body("moon", dts).transform_to(AltAz(location=pos))
# 太陽の位置（赤道座標から地平座標へ変換）
Spos = get_body("sun", dts).transform_to(AltAz(location=pos))
# 月―太陽間の位相角（rad)
SM_angle = Mpos.position_angle(Spos).value
# Lunar Ageに相当する数値として、内積を用いる。1なら新月(日食)、-1なら満月(月食)となる
earthquake["inner"] = [np.cos(rad) for rad in SM_angle]
#earthquake["inner"] = Inner
earthquake["moon_dist"] = Mpos.distance.value
earthquake["sun_dist"] = Spos.distance.value
earthquake["moon_direction"] = Mpos.az.degree

# 散布図
ggplot(earthquake, aes(x='Longitude', y='Latitude', color='inner', size='Magnitude')) + geom_point() + labs(color='Lunar Phase')

# タイル状に要約
tmp = earthquake.assign(
    Lon_binned=lambda d: pd.cut(d.Longitude, bins=np.arange(-180, 180, 10), right=True, labels=np.arange(-175, 175, 10)).astype(float),
    Lat_binned=lambda d: pd.cut(d.Latitude, bins=np.arange(-90, 90, 10), right=True, labels=np.arange(-85, 85, 10)).astype(float)
)
ggplot(tmp,
       aes(x='Lon_binned', y='Lat_binned', fill='inner', width=10, height=10)
      ) + geom_tile() + labs(x='Longitude', y='Latitude', fill='Lunar Phase')

kumamoto = earthquake["2016-4-10":"2016-4-17"]
kumamoto = kumamoto[(kumamoto["Latitude"]>31)&(kumamoto["Latitude"]<33)]
kumamoto.head()

Latitude    Longitude   Type    Depth   Depth Error Depth Seismic Stations  Magnitude   Magnitude Type  Magnitude Error Magnitude Seismic Stations  ... Source  Location Source Magnitude Source    Status  p_angle inner   moon_dist   sun_dist    moon_az moon_alt
datetime                                                                                    
2016-04-14 12:26:36+00:00   32.7880 130.7042    Earthquake  9.0 1.7 NaN 6.2 MWW NaN NaN ... US  US  US  Reviewed    136.528213  -0.725713   3.819877e+08    1.500738e+11    255.215900  49.755287
2016-04-14 15:03:47+00:00   32.6973 130.7204    Earthquake  8.0 1.7 NaN 6.0 MWW NaN NaN ... US  US  US  Reviewed    138.516836  -0.749150   3.854937e+08    1.500797e+11    277.541979  17.660531
2016-04-15 16:25:06+00:00   32.7906 130.7543    Earthquake  10.0    1.7 NaN 7.0 MWW NaN NaN ... US  US  US  Reviewed    133.309973  -0.685945   3.916934e+08    1.501228e+11    279.246432  9.941127
2016-04-15 16:45:57+00:00   32.9241 130.8091    Earthquake  10.0    0.9 NaN 5.7 MWW NaN NaN ... US  US  US  Reviewed    132.340430  -0.673534   3.922257e+08    1.501232e+11    281.839741  5.725034