
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)

<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


uwnd = data["uwnd"] #変数 "uwnd"の切り出し


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


#時間平均をとる
uwnd_mean_t = uwnd_sub.mean(dim=('time'))


### 最後に地図上に描画しましょう。

clevels = np.arange(-20,50,5) # 塗りつぶし用

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

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

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

# (2) 図を描く
uwnd_mean_t.plot.contourf(ax=ax,transform=ccrs.PlateCarree(),levels=clevels,cmap='rainbow',cbar_kwargs={'label':'U-velocity (m s$^{-1}$)'})
# transform = にはデータ自身の座標系を指定する; 大抵の場合ccrs.PlateCarree()としておけば問題ない

# (3) 目盛り
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)


# (4) 描画領域の設定
ax.set_extent([90,180,20,70],crs=ccrs.PlateCarree())


# (5) オプション
ax.coastlines() #海岸線
ax.gridlines(draw_labels=False) #罫線：ラベルはすでに上で描いたので"False"

plt.show()


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

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

# NearsidePerspective: 宇宙から見た絵
ax = fig.add_subplot(111,projection=ccrs.NearsidePerspective(central_latitude=40,central_longitude=140))

# (2) 図を描く
uwnd_mean_t.plot.contourf(ax=ax,transform=ccrs.PlateCarree(),levels=clevels,cmap='rainbow',\
                          cbar_kwargs={'label':'U-velocity (m s$^{-1}$)'})

### 目盛りはかけない。また全球データを使うので、範囲指定もできない###

# (5) オプション
ax.coastlines() #海岸線
ax.gridlines(draw_labels=False) #罫線：ラベルはすでに上で描いたので"False"

ax.set_global() # 全球データをつかうときのオマジナイ

ax.set_title('U-wind at 500 hPa')

plt.show()


uwnd_mean_xt = # ここを各自埋めよ


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

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

fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)

# ベタ塗り
uwnd_mean_xt.plot.contourf(ax=ax,cmap='rainbow',levels=clevels)


# コンター
ctr = uwnd_mean_xt.plot.contour(ax=ax,levels=cnt_levels,colors='black',yincrease=False,yscale="log")
# Y軸は上下反転させ,かつlog-scaleに(気圧は高度と共に指数関数的に減少するので)

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()


yst = 2018; yen = 2019 
#（yst = 2021; yen = 2023）に変更してください。　
#6時間間隔の再解析データは2021年以降のものしか（VDI上には）ないため、AMeDASの解析期間もそれに合わせる
dlist = []
for yy in range(yst,yen+1,1): 
    cyr = str(yy) # 文字列化
    file = "../../DB_2020/data/Maizuru/Maizuru_" + cyr  +".csv" #<---【注】パスを変更のこと
    
    df = pd.read_csv(file,encoding = "shift-jis",header=2,index_col=0,parse_dates=[0])
    # "encoding=":文字コードを指定、"header=2":2行目をヘッダー行（変数情報）に指定、"index_col=0": 0列目をインデックス行（時刻）に指定
    # ""parse_dates =[0]": 「日時」の型として読み込む（これにより、年・月・曜日といった情報を簡単に取り出せる） 
    
    df = df.iloc[2:,:] # ヘッダーの余計な部分（空白行など）を除いて上書き
    # これで一年分のデータのかたまりができた

    dlist.append(df) # listに要素(df)を追加していくことで

# pandas.concatで、データのかたまり（41年分）を結合する
dfm = pd.concat(dlist,axis=0) # axis=0 (時間軸方向）に結合


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

mask = snow > 0

print(snow[mask])

年月日
2018-01-13     1.0
2018-01-24     9.0
2018-01-25    31.0
2018-01-26    20.0
2018-01-27     1.0
2018-01-29     1.0
2018-01-30    15.0
2018-02-01     1.0
2018-02-04    13.0
2018-02-07     1.0
2018-02-08     1.0
2018-02-11     4.0
2018-02-12     3.0
2018-12-28     6.0
2018-12-29     3.0
2018-12-30     2.0
2019-01-26     2.0
2019-01-27     1.0
2019-02-11     1.0
Name: 降雪量合計(cm), dtype: float64


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年平均の日本上空(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 地図描画¶

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

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

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

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

課題¶

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

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

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

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

3.6【参考】ベクトル場¶

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

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