
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import xarray as xr
import cartopy.crs as ccrs
import cartopy.util as cutil
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
#from pandas.plotting import register_matplotlib_converters
#register_matplotlib_converters()

Warning: ecCodes 2.21.0 or higher is recommended. You are running version 2.17.0


#月平均東西風データを指定
file = "../../DB_2020/data/ncep.reanalysis.derived/pressure/uwnd.mon.mean.nc" # --> 【パスを適宜変更すること！！！！】 

#データの読み込み
data = xr.open_dataset(file)

#データ構造の確認
print(data)

#必要な変数（"uwnd"）の切りだし
uwnd = data["uwnd"]

<xarray.Dataset>
Dimensions:  (level: 17, lat: 73, lon: 144, time: 855)
Coordinates:
  * level    (level) float32 1e+03 925.0 850.0 700.0 ... 50.0 30.0 20.0 10.0
  * lat      (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
  * lon      (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * time     (time) datetime64[ns] 1948-01-01 1948-02-01 ... 2019-03-01
Data variables:
    uwnd     (time, level, lat, lon) float32 ...
Attributes:
    title:          monthly mean u wind from the NCEP Reanalysis
    description:    Data from NCEP initialized reanalysis (4x/day).  These ar...
    platform:       Model
    Conventions:    COARDS
    NCO:            20121012
    history:        Mon Jul  5 22:36:33 1999: ncrcat uwnd.mon.mean.nc /Datase...
    References:     http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reana...
    dataset_title:  NCEP-NCAR Reanalysis 1


# (1) 京都上空（135E, 35N)における風速の時間変動

# データの切り出し（2010-2018年, 500 hPaレベル）
uwnd_1D = uwnd.sel(lon = 135, lat = 35, time=slice('2010','2018'), level=500)

# 描画
uwnd_1D.plot()
plt.show()


# (1) 季節変化と年平均を図示する発展版

# 図の作成（1行2列のサブプロット）
fig, axes = plt.subplots(1, 2, figsize=(10, 5))

# 月情報を取得（1〜12月）
month_key = uwnd_1D["time"].dt.month

# 季節平均（各月ごとの平均）を計算してプロット
uwnd_1D.groupby(month_key).mean(dim="time").plot(ax=axes[0])
axes[0].set_title("Monthly Mean")

# 年平均（年ごとに平均）を計算してプロット
uwnd_1D.resample(time="A-DEC").mean(dim="time").plot(ax=axes[1])
axes[1].set_title("Inter-annual variaton (annual-mean)")

# 表示
plt.tight_layout()
plt.show()


# (2) データの切りだし・平均
uwnd_mean_t = uwnd.sel(time=slice('2010','2018'), level=500).mean(dim="time")

# カラーマップのレベル（等高線の値の範囲と間隔）
clevels = np.arange(-20, 50, 5)

# 図の作成
fig = plt.figure(figsize=(8, 5))

# 地図投影（PlateCarree：正距円筒図法、中心経度を180度に設定）
ax = fig.add_subplot(111, projection=ccrs.PlateCarree(central_longitude=180))

# 等高線塗りつぶしプロット
uwnd_mean_t.plot.contourf(
    ax=ax,
    transform=ccrs.PlateCarree(),  # データの座標系を指定（通常PlateCarreeでOK）
    levels=clevels,
    cmap='rainbow',
    cbar_kwargs={'label': 'U-velocity (m s$^{-1}$)'}
)

# 緯度・経度の目盛り設定
xticks = np.arange(0, 361, 30)
yticks = np.arange(-90, 91, 30)
ax.set_xticks(xticks, crs=ccrs.PlateCarree())
ax.set_yticks(yticks, crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))
ax.yaxis.set_major_formatter(LatitudeFormatter())

# 描画範囲の設定（経度90〜180度、緯度20〜70度）
ax.set_extent([90, 180, 20, 70], crs=ccrs.PlateCarree())

# 地図のオプション：海岸線とグリッド線
ax.coastlines()
ax.gridlines(draw_labels=False)

# 表示
plt.show()


lon = uwnd_mean_t["lon"]
lat = uwnd_mean_t["lat"]

# 描画設定
fig = plt.figure(figsize=(5, 5))

# projectionのオプションを変更
ax = fig.add_subplot(111, projection=ccrs.NearsidePerspective(central_latitude=90, central_longitude=140))

#データをサイクリックに(360Eの飛びを埋める)
data_cyclic, lon_cyclic = cutil.add_cyclic_point(uwnd_mean_t.values, coord=lon.values)

#↑add_cyclic_pointはNumPy配列に適用できるので、.valuesとしてデータ配列を渡している

# 塗りつぶし図
contour = ax.contourf(
    lon_cyclic,
    lat,
    data_cyclic,
    transform=ccrs.PlateCarree(),
    cmap="jet",
    levels=np.linspace(-10, 30, 9)
)

ax.coastlines()
ax.gridlines(draw_labels=False)
ax.set_global()

plt.colorbar(
    contour,
    label = 'U-velocity (m s$^{-1}$)',
)

plt.show()


data_mean_xt =  # ここを各自埋めよ


## ここから描画-----

# コンターレベルを設定
clevels = np.arange(-20,80,10) # 塗りつぶし用
cnt_levels = np.arange(-20,80,10) #コンター用

fig,ax = plt.subplots(figsize=(8,6))

# ベタ塗り
data_mean_xt.plot.contourf(ax=ax,cmap='rainbow',levels=clevels,yincrease=False,yscale="log")
# Y軸は上下反転させ,かつlog-scaleに(気圧は高度と共に指数関数的に減少するので)

# コンター
ctr = data_mean_xt.plot.contour(ax=ax,levels=cnt_levels,colors='black',yincrease=False,yscale="log")

# コンターラベル
ctr.clabel(fmt='%3i',fontsize=12,colors='black')

# 軸の書式設定
ax.tick_params(labelsize=15)
ax.set_ylabel('Pressure (hPa)',size=15)
ax.set_xlabel('Latitude ($^\circ$N)',size=15) # $...$はTexでの数式の表現

plt.show()


# 指定時刻・高度のデータを取得して次元を縮小
data_cut = uwnd.sel(time="2010-01", level=500).squeeze()

### 方法 (1): Xarrayを用いた簡便な偏差計算 ###
anom_xr = data_cut - data_cut.mean(dim="lon")
### 方法 (1) ここまで ###

### 方法 (2): NumPy配列に変換して手動で処理する場合 ###

# NumPy配列に変換
data_v = data_cut.values

# データの形状を確認し、同じサイズの空配列を作成
ny, nx = data_v.shape
anom_v = np.empty((ny, nx))

# 経度方向の平均を計算
mean_v = data_v.mean(axis=1)

# 平均を拡張して"ブロードキャスト"し、偏差を計算
anom_v = data_v - mean_v[:, np.newaxis]

# ↑ 以下のように書くとエラーになる：
# anom_v = data_v - data_v.mean(axis=1)
# 理由：配列の形が (73,144) と (73,) で異なり、計算できないため。
# mean_v[:, np.newaxis] により (73,1) の形に変換し、経度方向に平均値を並べることで対応。
# この処理を「ブロードキャスト」という。

# XarrayのDataArrayとして再構築（座標・次元情報を保持）
anom_xr_new = xr.DataArray(anom_v, coords=data_cut.coords, dims=data_cut.dims)
### 方法 (2) ここまで ###

### 結果をプロット（地図描画は省略） ###
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
anom_xr.plot.contourf(ax=axes[0])
anom_xr_new.plot.contourf(ax=axes[1])
plt.show()


def read_csv(START_YEAR,END_YEAR,STATION):
#START_YEAR = 1960
#END_YEAR = 2023
    DATA_DIR = "../../DB_2020/data/"+STATION+"/"+STATION+"_" #<---【注】"data/AMeDAS/Kyoto/Kyoto_"に変更のこと
    ENCODING = "shift-jis"

    # データを格納するリスト
    dataframes = []

    for year in range(START_YEAR, END_YEAR + 1):
        file_path = DATA_DIR + str(year) + ".csv"

        df = pd.read_csv(
            file_path,
            encoding=ENCODING,
            header=2,
            index_col=0,
            parse_dates=[0]
        )

        # 余計なヘッダー行を除去
        df = df.iloc[2:, :]

        dataframes.append(df)

    # すべての年のデータを結合（時間軸方向に）
    dfm = pd.concat(dataframes, axis=0)
    
    return dfm

yst = 2018; yen = 2019
#注：演習では yst = 2021; yen = 2023と変更してください（６時間間隔の再解析データが2021年以降しかないため）

dfm = read_csv(yst,yen,"Maizuru")


snow = dfm["降雪量合計(cm)"]

mask = snow > 0

#print(snow[mask])


filesSLP = []; filesAIR = []
for yy in range(yst,yen+1):
    cyr = str(yy) # 文字列化
    fileSLP = "../../DB_2020/surface/slp."+ cyr  +".nc"#<---【注】パスを変更のこと
    fileAIR = "../../DB_2020/pressure/air."+ cyr  +".nc"#<---【注】パスを変更のこと
    filesSLP.append(fileSLP)
    filesAIR.append(fileAIR)


slp = xr.open_mfdataset(filesSLP,combine="by_coords")["slp"].sel(lon=slice(110,160),lat=slice(60,20)) 
# xr.open_mfdataset([file1,file2,...],combine="by_coords")とすると、複数のファイルを纏めて開くことができます
#（ただし、時刻が連続的に並んでいる場合のみ）

air = xr.open_mfdataset(filesAIR,combine="by_coords")["air"].sel(lon=slice(110,160),lat=slice(60,20),level=850) 

#日平均値に
slp_daily = slp.resample(time="D").mean(dim="time")
air_daily = air.resample(time="D").mean(dim="time")

#舞鶴で積雪があった日のデータを抜き出す
slp_snow = slp_daily[mask.values,:,:].mean(dim="time")
air_snow = air_daily[mask.values,:,:].mean(dim="time")
## maskを再解析データの時間軸（次元0）に適用；slp_daily[mask]とやってもうまくいかない


# 描画
fig = plt.figure(figsize=(8,5))

# (1) 領域を用意（地図投影オプションをつけて）

# PlateCaree: 正距円筒図法, central_longitude: 図の中心の経度
ax = fig.add_subplot(111,projection=ccrs.PlateCarree(central_longitude=135))

(air_snow-273.15).plot.contourf(ax=ax,levels=np.linspace(-30,20,11),cmap="rainbow",transform=ccrs.PlateCarree()) #unit: K -> degrees C
ctr = (slp_snow/100).plot.contour(ax=ax,levels=np.linspace(980,1050,15),transform=ccrs.PlateCarree()) #unit: Pa -> hPa

#ctr.clabel(fmt='%3i',fontsize=12,colors='black')

ax.coastlines()
plt.show()


fileu = "../../DB_2020/data/ncep.reanalysis.derived/pressure/uwnd.mon.mean.nc" #要変更
filev = "../../DB_2020/data/ncep.reanalysis.derived/pressure/vwnd.mon.mean.nc" #要変更

# ファイルの読み込み
U = xr.open_dataset(fileu)
V = xr.open_dataset(filev)

p_cut = 500; time_cut = '2017-01'
Usub = U['uwnd'].sel(level=p_cut,time=time_cut).squeeze()
Vsub = V['vwnd'].sel(level=p_cut,time=time_cut).squeeze()

# low-resolution data
# (経度-緯度方向は適当に間引いて読み込む; そうしないと、描いたベクトルがギチギチになるので）
dx = 4; dy = 2
U_lr = Usub[::dy,::dx].values #[::c]は "c毎にデータを取得する"]の意味, ".values"はndarrayへの変換（軸情報を落として単なる数値配列に）
V_lr = Vsub[::dy,::dx].values

# 軸もそれに合わせて間引く
lon_lr = Usub['lon'][::dx].values; lat_lr = Vsub['lat'][::dy].values


# ここから描画
fig = plt.figure(figsize=(10,7))
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
#ax = plt.axes(projection=ccrs.NearsidePerspective(central_latitude=30,central_longitude=140))

# ベクトル場を描く
Q = ax.quiver(lon_lr,lat_lr,U_lr,V_lr,color='blue',transform=ccrs.PlateCarree())

# 単位ベクトルを描く
ax.quiverkey(Q,0.9,0.8,20,label='20 m s$^{-1}$',labelpos='E',coordinates="figure")

# 目盛り
xticks = np.arange(13)*30
yticks = -90 + np.arange(7)*30
ax.set_xticks(xticks,crs=ccrs.PlateCarree())
ax.set_yticks(yticks,crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)


ax.coastlines()
#ax.stock_img() # 地形を描く
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)

#plt.savefig('test_vector.png') #画面出力だとユニットベクトルが見えないが、保存データにはちゃんと描かれている

plt.show()

C:\Users\takas\Anaconda3\lib\site-packages\cartopy\mpl\geoaxes.py:1841: UserWarning: Some vectors at source domain corners may not have been transformed correctly
  u, v = self.projection.transform_vectors(t, x, y, u, v)

3. 全球大気再解析データ（４次元データ）の解析¶

3.0 はじめに¶

3.1 データの読み込み：NetCDFデータを読み込む¶

3.2 データを切り出す¶

3.2 統計解析・時系列解析¶

3.4 描画（ラインプロット・コンタープロット）¶

$\ \ \ \ $ 3.4.1 二次元プロットの方法¶

$\ \ \ \ $ 3.4.2 地図描画¶

3.5 【応用例】コンポジット解析¶

3.6 【参考】ベクトル場プロット¶

3.0 はじめに¶

XArrayについて¶

...さらに凝った解析や描画を行うには...¶

【練習課題 （１）】 2010-2018年の月平均東西風データ（500 hPa 気圧面～高度 5km）を用いて、¶

(a) 京都上空（135E, 35N)における時間変動¶

(b) 日本領域(500hPa)の東西風の分布（緯度-経度断面）¶

を図示する。¶

3.1 データの読み込み：NetCDFデータを読み込む、変数を切り出す¶

3.2 データを切り出す¶

3.2.1 ラベルで範囲指定¶

(1) 軸の名前と値を指定した切り出し: .sel(), .isel() （お勧め！）¶

(2) 各軸のラベル・整数のインデックス値を用いた切り出し：.loc[ ]、[ ]¶

3.2.2 条件指定して切り出し¶

◆ 時刻データを元に条件指定する際の注意点¶

【補足】 .squeeze() メソッド: 大きさ１の次元を削除¶

3.3 統計解析・時系列解析¶

◆ 統計解析¶

◆ 時系列データ解析¶

・resampleの例（年平均データの作成）：¶

・groupbyの例（各月のコンポジット）：¶

◆ さらに凝った解析を行うには？？¶

3.4 データの二次元描画¶

3.4.1 二次元描画：.plot.contourf（塗りつぶし）, .plot.contour（コンター）¶

【コンターラベルの作成】¶

3.4.2 地図描画¶

◆ 作成の手順（通常の描画と違う点を中心に）¶

別の地図投影も試してみよう（宇宙から見た図）。¶

ここでは新たに：¶

・360Eのデータの飛びを埋める¶

・XArrayに頼らない描画方法を使ってみます¶

【練習課題（2）】2010-2018年データを用いて、時間・東西平均した東西風を緯度-高度（気圧）断面に描いてみよ（一行目の処理を埋めよ）¶

【練習課題（3）】各月（例えば、１月と７月）の平均値を用いて（１）と同様の図を描いてみよ。（「月」の情報の取り出し方は、3.2.2を参照）¶

より凝った解析をしたいとき¶

【練習課題（4）】2010年01月の東西風（500 hPa面）について、東西平均からの偏差を計算して図示せよ。¶

課題¶

【課題 3-1】6時間間隔のジオポテンシャルハイトデータ（2019年01月、北緯40度）を用いてホフメラー図を描いてみよ（＋2 点）。高度によって(例：500 hPaと10 hPa）でどのように特徴が変るだろうか。¶

【課題 3-4】赤道上空東西風を用いて時間-高度（気圧）断面図を作成し、東西風準二年周期変動が見られることを確かめよ。（＋2 点）¶

【課題3-5】（+5 点）¶

地上気圧の年平均時系列を作成し、これらを用いてダーウィンを基準とした一点相関図を作成して南方振動の空間構造を調べよ（参考：伊藤・見延 (2010) 図6.8）。¶

(1) まずは、ダーウィン(135$^\circ$E, 12.5$^\circ$S) と タヒチ(210$^\circ$E, 17.5$^\circ$S)の時系列を描いて相関係数を計算してみよ。¶

(2) 続いて、全格子点についてダーウィンとの相関係数を計算し、それを地図上に二次元マッピングすることで一点相関図を作成する。¶

3.5 【応用例】コンポジット解析¶

【練習課題（5）】舞鶴で積雪が観測された日の「地上気圧場」と「850hPa面気温場」のコンポジット図を作成せよ。舞鶴のデータは日別のAMeDASデータ、地上気圧・気温場は6時間間隔の再解析データを用いることにする。¶

3.6【参考】ベクトル場¶

3.　全球大気再解析データ（４次元データ）の解析¶

【練習課題（１）】 2010-2018年の月平均東西風データ（500 hPa 気圧面～高度 5km）を用いて、¶

(1) 軸の名前と値を指定した切り出し: .sel(), .isel()　（お勧め！）¶

(1) まずは、ダーウィン(135$^\circ$E, 12.5$^\circ$S) とタヒチ(210$^\circ$E, 17.5$^\circ$S)の時系列を描いて相関係数を計算してみよ。¶