目的
気圧場が与えられたとき、COL (Cut-Off Law，ちょっと台風に似ている)を形成してしているような成分を分離したい。つまり気圧場の極小値をとってきたいな〜というのが今回のモチベーションです。厳密に閉鎖領域をもとめるのがとてもめんど……メモリ容量的に厳しかったので、上手いこと取ってこれないかなと思っていました。
フィールド全体にスムージングをかけて、極大値の逆バージョン的なのりで極小を求めようとしたのですけれどこれが案外うまく行かなくて。あと２次元で極値を求めるコードがほとんどネットの海に転がっていなかったので少々苦戦しました。最終的にガウシアンフィルタ→極小値フィルタで解決しました。下にpythonコードを載せておきます。ご活用ください。
もとのコードをきれいにして貼っているのでもしかするとどこかエラーが出るかもしれませんが、まあご愛嬌ということで……。
ところでコードとスクリプトとプログラムの堺って年々曖昧になっていってません？ 動的言語が登場したのでさもありなんといった感じではありますが……
あ、このコード自体は別に気圧場じゃなくても閾値さえちゃんとしておけば動きます。
python スクリプト
import numpy as np
from scipy.ndimage import gaussian_filter, minimum_filter
def local_minimum_pressure\
(field, lsmask=None, footprint=np.ones((5,5)), sigma=2, foot = 1000):
'''
<input>
field ::: pressure field (2d array-like object)
lsmask ::: land-sea mask. 0 for land, 1 for sea.
footprint ::: smoothing filter size.
sigma ::: smoothing constant. 2 is standard deviation for gaussian kernel
foot ::: minimum value for detection.
<output>
peaks_index ::: index for local minimum
'''
field = np.array(field)
assert field.ndin == 2, \
'Input array dimension is %d. Please confirm the dimension to be 2.'%field.ndim
if lsmask == None:
lsmask = np.ones(field.shape)
assert lsmask.shape == field.shape, \
'array shape does not match for field with shape str(field.shape) \
and lsmask with str(lsmask.shape).'
filtered_field = gaussian_filter(field, sigma=sigma)
local_min = minimum_filter(filtered_field, footprint=footprint, mode='constant')
local_min = np.ma.masked_where(field > foot, local_min)
local_min_mask = (filtered_field == local_min)
peaks_index = np.where((local_min_mask == True))
return peaks_index
プロット例
ECMWF（欧州中期予報センター）のERA5再解析データから気圧場および陸水マスクをダウンロードしてきた前提で、プロットのためのコードを載せておきます。
#(continuum from previous code)
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from netCDF4 import Dataset
fname1 = 'pressure.nc' # please insert appropriate file name
fname2 = 'lsmask.nc'
nc1 = Dataset(fname1, 'r')
nc2 = Dataset(fname2, 'r')
lon = nc1['longitude'][:]
lat = nc1['latitude'][:]
field = nc1['sp'][:]/100 # 'sp' means surface pressure. convert [Pa] to [hPa]
lsmask = nc2['lsmask'][:]
### ! caution !
### lsmask and pressure field resolution is different in the case of ECMWF.
### in this code, to solve this problem, reduced the resolution of pressure field.
peaks_index = local_minimum_pressure(field[::2,::2], lsmask=lsmask,\
footprint=np.ones((5,5)), sigma=2, foot = 1000)
LON, LAT = np.meshgrid(lon[::2], lat[::2])
LON = np.ma.masked_where(lsmask == 0, LON)
LAT = np.ma.masked_where(lsmask == 0, LAT)
# gentrate figure
fig = plt.figure(figsize = (8, 4))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
cont = ax.contourf(LON, LAT, field[::2,::2] * lsmask,\
levels = np.linspace(950,1030,100), cmap='rainbow')
fig.colorbar(cont, ax = ax, orientation='vertical', ticks = np.arange(950, 1050, 10),\
shrink=0.7)
cs = ax.contour(LON, LAT, field[::2,::2] * lsmask, levels = np.arange(950,1050,10),\
colors='black', linewidths=1, linestyles = ['-', '--'])
ax.clabel(cs, fontsize=15, colors="black")
ax.scatter(LON[peaks_index[0], peaks_index[1]], LAT[peaks_index[0], peaks_index[1]],\
size=100, marker='+', color='red', zorder=10)
fig.show()
結果
我ながら結構いい感じですね
参考
・ガウシアンフィルター
・極小値フィルター
・python で２次元場の極値をもとめる
・ERA5再解析データ
Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J., Nicolas, J., Peubey, C., Radu, R., Rozum, I., Schepers, D., Simmons, A., Soci, C., Dee, D., Thépaut, J-N. (2023): ERA5 hourly data on single levels from 1940 to present. Copernicus Climate Change Service (C3S) Climate Data Store (CDS), DOI: 10.24381/cds.adbb2d47 (Accessed on 26-Jun-2024)
・Cut-Off Lawについて（どうやら1800年代のスペインあたりから文献が見られるらしいのだが、あいにくそこまでの言語能力がないのでwikipediaのページを貼っておく）