Pythonによる気象・気候データ解析I 4章章末問題

[Pythonによる気象・気候データ解析I: Pythonの基礎・気候値と偏差・回帰相関分析](https://amzn.to/3yvtPCG)

<a href="https://amzn.to/3yvtPCG" target="_blank"><img src="https://m.media-amazon.com/images/I/71u7jHR9KgL._SY522_.jpg" width=150 /></a>


の章末問題を一つずつ解いていきます。

# A

```python
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
import japanize_matplotlib

tokyo_temp = np.genfromtxt('drive/MyDrive/Tokyo_temp.csv', delimiter=',', usecols=(0, 1, 2))
y = tokyo_temp[:, 0]
m = tokyo_temp[:, 1]
temp = tokyo_temp[:, 2]

y_start = 1960
y_end = 1989
temp = temp[(y_start <= y)*(y <= y_end)]
m = m[(y_start <= y)*(y <= y_end)]
y = y[(y_start <= y)*(y <= y_end)]

month = np.arange(1, 13, 1)
temp_clim = np.zeros(12)

for mm in range(1, 13):
  temp_clim[mm - 1] = np.nanmean(temp[m == mm], 0)
plt.plot(month, temp_clim)
plt.xticks(month)
plt.show()
```

まず単純に平均値のグラフを描いた。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/2a0d5632-d794-03d4-e6d4-517745da919a.png)

```python
mon = np.arange(y_start, y_end+1, 1/12)
tmt = temp.shape

tempa = np.zeros((tmt))
varsum = 0

for yy in range(y_start, y_end):
  for mm in range(1, 13):
    tempa[(y==yy)*(m==mm)] = temp[(y==yy)*(m==mm)] - temp_clim[mm-1]
    varsum += (temp[(y==yy)*(m==mm)] - temp_clim[mm-1])**2

std = np.sqrt(varsum/tmt)

temp_winter = tempa.copy()
temp_summer = tempa.copy()
temp_winter[(m >= 3)*(m <= 11)] = np.nan # 12-2以外はNaN
temp_summer[(m >= 9) | (m <= 5)] = np.nan # 6-8以外はNaN

plt.figure(figsize=(20, 6))
plt.plot(mon, tempa, c = 'black', lw=1)
plt.plot(mon, temp_summer, c = 'red', lw=3, label='夏季')
plt.plot(mon, temp_winter, c = 'blue', lw=3, label='冬季')

plt.hlines(0.4, y_start, y_end, color = 'red', linestyle='solid')
plt.hlines(0.3, y_start, y_end, color = 'skyblue', linestyle='dashed')
plt.hlines(0, y_start, y_end, color = 'black')
plt.hlines(-0.3, y_start, y_end, color = 'pink', linestyle='dashed')
plt.hlines(-0.4, y_start, y_end, color = 'blue', linestyle='solid')

# 年はWikipediaから
mousyo = [1961, 1967, 1973, 1978, 1984, 1985, 1990]
mousyo = list(map(lambda x: x + 0.5, mousyo)) # 夏季にplotのためxをずらす

reika = [1965, 1969, 1974, 1976, 1977, 1980, 1981, 1982, 1983, 1986, 1987, 1988, 1989]
reika = list(map(lambda x: x + 0.5, reika)) # 夏季にplotのためxをずらす

danto = [1960, 1969, 1972, 1973, 1979, 1987, 1988, 1989]
kanto = [1961, 1963, 1967, 1968, 1970, 1971, 1974, 1975, 1977, 1978, 1980, 1981, 1984, 1985, 1986]


plt.scatter(mousyo, [3] * len(mousyo), color='red', label='猛暑')
plt.scatter(reika, [-3] * len(reika), color='skyblue', label='冷夏')
plt.scatter(danto, [3] * len(danto), color='pink', label='暖冬')
plt.scatter(kanto, [-3] * len(kanto), color='blue', label='寒冬')

plt.xlim(y_start, y_end)
plt.legend(loc='upper left', bbox_to_anchor=(1, 1))

plt.xlabel('Year')
plt.ylabel('Temperature[℃]')
plt.title('Temperature of Tokyo')
plt.show()
```
- 赤線が上に出ているのが夏季に平年値より気温が高い→猛暑
- 赤線が下に出ているのが夏季に平年値より気温が低い→冷夏
- 青線が上に出ているのが冬季に平年値より気温が高い→暖冬
- 青線が下に出ているのが冬季に平年値より気温が低い→寒冬
- 上記とWikipediaの（日本全体での）猛暑、冷夏、暖冬、寒冬を比較
    - 日本全体と東京のみの気温の比較であり、若干ズレるが概ね妥当なところかと思う。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/116c261b-2624-2a6a-6a95-5ccbc64b7d62.png)

# B

```python
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
import japanize_matplotlib

loadfile = 'drive/MyDrive/sst_OISST.npz'
sst_dataset = np.load(loadfile)
sst = sst_dataset['sst']
lon2 = sst_dataset['lon2']
lat2 = sst_dataset['lat2']
y = sst_dataset['y']
m = sst_dataset['m']

[imt, jmt, tmt] = sst.shape

sst_clim = np.zeros((imt, jmt, 12))

for mm in range(1, 13):
  sst_clim[:, :, mm - 1] = np.mean(sst[:, :, m == 12], 2)

sst_1982_1999 = sst[:, :, (1982 <= y) * (y <= 1999) * (m == 12)]
sst_2000_2019 = sst[:, :, (2000 <= y) * (y <= 2019) * (m == 12)]

sst_1982_1999_12mean = np.mean(sst_1982_1999[:, :, :], 2)
sst_2000_2019_12mean = np.mean(sst_2000_2019[:, :, :], 2)

vmin = -5
vmax = 35
vint = 5


cm = plt.get_cmap('seismic')
cs = plt.contourf(lon2, lat2, sst_1982_1999_12mean,
                    cmap=cm, norm=Normalize(vmin=vmin, vmax=vmax),
                    levels=np.arange(vmin, vmax + vint, vint), extend='both')
plt.colorbar(cs)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
clim_title = '1982_1999平均 12月'
plt.title(clim_title)
plt.xlim(0, 360)
plt.ylim(-90, 90)
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/57814162-8b68-74d3-68a9-91d5aad3cc6f.png)

```python
cm = plt.get_cmap('seismic')
cs = plt.contourf(lon2, lat2, sst_2000_2019_12mean,
                    cmap=cm, norm=Normalize(vmin=vmin, vmax=vmax),
                    levels=np.arange(vmin, vmax + vint, vint), extend='both')
plt.colorbar(cs)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
clim_title = '2000_2019平均 12月'
plt.title(clim_title)
plt.xlim(0, 360)
plt.ylim(-90, 90)
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/d769e79a-fbb1-2fc2-66f8-24a7c649327e.png)

```python
sst_diff_12mean = sst_2000_2019_12mean - sst_1982_1999_12mean
plt.hist(sst_diff_12mean.flatten(), bins=50)
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/8a64879f-026f-1f2a-f308-532c30edc577.png)

```python
vmin = -2.0
vmax = 2.0
vint = 0.2

sst_diff_12mean = sst_2000_2019_12mean - sst_1982_1999_12mean

cm = plt.get_cmap('seismic')
cs = plt.contourf(lon2, lat2, sst_diff_12mean,
                    cmap=cm, norm=Normalize(vmin=vmin, vmax=vmax),
                    levels=np.arange(vmin, vmax + vint, vint), extend='both')
plt.colorbar(cs)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
clim_title = '2000_2019平均 と 1982_1999平均の差分 12月'
plt.title(clim_title)
plt.xlim(0, 360)
plt.ylim(-90, 90)
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/5a1627ef-5e7b-5619-b439-331951de1181.png)

- 二期間のそれぞれのグラフから、
    - 温度が高い海域が広がっているように見える
- ヒストグラムから、
    - 全体に温度が上がっている
- 二期間の差分のグラフから、
    - 全体に温度が上がっている海域が多いが下がっている海域もある。
        - 大陸周辺で特に上がっている海域と下がっている海域が見られる。

# C

```math
\begin{align}
偏差&=x-\bar{x} \\
偏差の平均&=\frac{1}{n}\sum_{i=1}^n{(x_i-\bar{x})} \\
ここで、\bar{x}&=\frac{1}{n}\sum_{i=1}^n{x_i} だから\\
偏差の平均&=\frac{1}{n}\sum_{i=1}^n{(x_i-\frac{1}{n}\sum_{i=1}^n{x_i})}\\
&=\frac{1}{n}{\{(x_1+x_2+\cdots+x_n)-(x_1+x_2+\cdots+x_n)\}} \\
&=0
\end{align}
```

# D
- 下記データを用いた
    - https://www.jnto.go.jp/statistics/data/visitors-statistics/
        - 国籍/月別 訪日外客数（2003年～2024年）（Excel）

```python
import pandas as pd
import matplotlib.pyplot as plt

df = pd.DataFrame()
# 各シートから「総数」の各月の値のみを抽出。行が年数、列が月数に転置してデータフレーム化
for year in range(2003, 2024):
  df_org = pd.read_excel('drive/MyDrive/since2003_visitor_arrivals_May_2024.xlsx', sheet_name=str(year), skiprows=3)
  df[year] = df_org.loc[0, ['1月', '2月', '3月', '4月', '5月', '6月', '7月', '8月', '9月', '10月', '11月', '12月']]
df = df.T

# 月毎（列毎）の平均値
df_mean = df.mean()

# 月毎の平均値のプロット
df_mean.plot()
plt.show()
```

- 1月、4月、7月、8月、10月、12月、特に7月が多い
    - 四季が目立つ時期に多い？
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/28f98551-5d20-48f8-bb7f-779287dcd4a1.png)

```python
# 平均値との差
df_july = df['7月'] - df_mean['7月']

plt.plot(np.arange(2003, 2024), df_july)
plt.title('7月の訪日外客数の平年値との差')
plt.xticks(np.arange(2003, 2024, 2))
plt.show()
```
- 7月に着目して、平年値との差を各年でプロットした
    - 2013年頃から訪日外客数が急増した
        - 2012-2013年頃のクールジャパン戦略、2013-2015年の急速な円安もあると思われる
    - 2020-2022年はコロナ禍により急落したが2023年には回復が見られる

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/116049/01e5c698-48ec-10b7-175c-5eb18d27ad5b.png)

---

- [Pythonによる気象・気候データ解析I 1章章末問題](https://qiita.com/aokikenichi/items/755925f2f16e4f9a022c)
- [Pythonによる気象・気候データ解析I 2章章末問題](https://qiita.com/aokikenichi/items/dc653c07059b7ac63c0c)
- [Pythonによる気象・気候データ解析I 3章章末問題](https://qiita.com/aokikenichi/items/4be2cb4a8ebb4648fb4d)
- [Pythonによる気象・気候データ解析I 5章章末問題](https://qiita.com/aokikenichi/items/538c0ec08d47624d9475)
- [Pythonによる気象・気候データ解析II 1章章末問題](https://qiita.com/aokikenichi/items/3dd97d989cdd31043e2f)