Go to Qiita Advent Calendar 2024 Top
LoginSignup
1
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

お題は不問！Qiita Engineer Festa 2024で記事投稿！
Qiita Engineer Festa20242024年7月17日まで開催中！
@aokikenichi(青木 健一)

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

Last updated at Posted at 2024-07-14

Pythonによる気象・気候データ解析I: Pythonの基礎・気候値と偏差・回帰相関分析

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

A

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

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

B

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

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

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

image.png

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

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

C

\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

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 
# 平均値との差
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

1
2
2

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?