雨量予想をプロットしてみた

xrayとcartopyで雨量予想をプロットする。ピークの値を表示してそれっぽくプロットする方法を調べてみた。

ピークの検出

以下のFAQが参考になった。

stackoverflow.com

隣接する格子でのピークを参照するためにScipyのndimage.filters.maximum_filterを使った。

maximum_filterで隣接格子の最大値でフィルターしたデータを用意して、元データと比較すればピークの格子を特定できる。

import numpy as np
import numpy.random as random
from scipy.ndimage.filters import maximum_filter

# ランダムデータ
data = random.rand(25,25)
# 各格子を隣接５格子のうち最大の値でフィルターして置き換える
local_max = maximum_filter(data, footprint=np.ones((5,5)), mode='constant')
# 元のデータとフィルターしたデータで値が同じ格子（ピーク格子）以外をマスクする
peak = np.ma.array(data, mask=~(data==local_max))

import matplotlib.pyplot as plt
plt.subplot(2,2,1)
plt.pcolormesh(data, vmin=0, vmax=1)
plt.title('Original')
plt.colorbar()

plt.subplot(2,2,2)
plt.pcolormesh(local_max, vmin=0, vmax=1)
plt.title('Filtered')
plt.colorbar()

plt.subplot(2,2,3)
plt.pcolormesh(peak, vmin=0, vmax=1)
plt.title('Peak')
plt.colorbar()

plt.show()

プロットしてみる

メソモデルの１時間雨量の予想をプロットしてみた。

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.patheffects as PathEffects

import xray
import numpy as np

## select data
# １時間雨量から東日本の範囲をxray.dataArrayで取得
ncfile = 'http://opendap.rish.kyoto-u.ac.jp/opendap/hyrax/jmadata/gpv/latest/MSM-S.nc'
ds = xray.open_dataset(ncfile)
darray = ds['r1h'].sel(lat=slice(38,32),lon=slice(133,143))[18,:,:]

## for annotation
# 数値プロットのために座標を2dのnumpy.ndarrayにしておく
ylab, xlab = darray.dims
xval = darray[xlab].values
yval = darray[ylab].values
xval, yval = np.meshgrid(xval, yval)

## detect peak
# 隣接15格子のピークかつR1=1mm以上の座標を算出
from scipy.ndimage.filters import maximum_filter

zval = darray.to_masked_array(copy=False)
local_max = maximum_filter(zval, footprint=np.ones((15,15)), mode='constant')
local_max[0,:]  = np.ma.masked
local_max[-1,:] = np.ma.masked
local_max[:,0]  = np.ma.masked
local_max[:,-1] = np.ma.masked
mask = (zval == local_max) & (zval >=1)

## for colormap
colors = ['#FFFFFF', '#00FFFF', '#00B0FF', '#0070FF', '#228B22', '#00FF00', '#FFFF00', '#FF8000', '#FF0000', '#FF00FF']
levels = [1,3,5,10,20,30,40,50,80]

## plot
plt.figure(figsize=(8,5))
ax = plt.axes(projection=ccrs.Mercator())
transform = ccrs.PlateCarree()._as_mpl_transform(ax)
darray.plot(ax=ax, levels=levels, colors=colors, extend='both', transform=transform)
for x, y, val in zip(xval[mask].flat, yval[mask].flat, zval[mask].flat):
    val = '{}'.format(int(val))
    ax.plot(x, y, '+', color='k', mew=1, ms=6, transform=transform)
    ax.annotate(val, xy=(x,y), xytext=(1, -1), xycoords=transform,
                textcoords='offset points',  ha='left', va='top', fontsize=10,
                path_effects=[PathEffects.withStroke(linewidth=2, foreground='w')])
ax.coastlines('10m')
plt.show()

メモ

地図座標への変換はCartopyを使っている。matplotlibのannotateを使って変換する場合は次のようにすればいいらしい。

ax = plt.axes(projection=ccrs.Mercator())
transform = ccrs.PlateCarree()._as_mpl_transform(ax)
ax.annotate(val, xy=(x,y), xycoords=transform)

こちらのFAQを参考にした。 stackoverflow.com