1. 必要なデータ
- モデルアウトプットのトラックデータ（本記事ではcsv）
- モデルアウトプット（本記事ではscaleのnetcdfファイル）
2.コード
python
##########################################################################
#
# This python program is makeing figure of Shade; Cd, Contour; MSLP by SCALE netcdf output.
#
# Hiroaki YOSHIOKA, yoshioka-hiroaki-sn@ynu.ac.jp
# Research fellow at Yokohama National University, Typhoon science and technology Research Center
#
# --Released History--
# Jan,12,2024 - 1st released ; Making distribution of Cd and MSLP.
# ; You can not handle ensemble members(s_lat,e_lat,s_lon,e_lon).
# ; Latitude and longitude settings are not implemented.
# ???,??,???? - 2nd released ;
#
##########################################################################
import os
import numpy as np
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from matplotlib.colors import LinearSegmentedColormap
import shapely.geometry as sgeom
import metpy.calc as mpcalc
import scipy.ndimage as ndimage
import tropycal.tracks as tracks
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
from cartopy.vector_transform import vector_scalar_to_grid
##########################################################################
# Data set and TC case
Model_name = 'scale' # scale only
TC_name = 'faxai' #
# Input part; need parameter, file name, file path and so on.
netcdf_dir_path = '/mnt/data1/model/scale/faxai/fnl/flow/'
case_names = ('CTL','Cd1030')
legends = ('CTL','x3.0-R200km') #case_names #('CTL','x3.0-R025km')
s_parameter = 'diff_cd'
c_parameter = 'cdr-C_track'
n_time = 24
s_lat = 0
e_lat = 0
s_lon = 0
e_lon = 0
input_dir_path = './trackdata_flow/'
input_csv_name = Model_name+'-'+TC_name
input_csv_type = '.csv'
parameter = 'mslp'
n_track_time = 73
# Best Track plot from IBTrACS
best_track_plot_swich = False # True or false
# Output part;
experiment_name = 'cd_value' # ' '
output_dir_path = './PI20240501/'
output_parameter = 'S_'+s_parameter+'-C_'+c_parameter # use for sub directry name
#output_fig_name = Model_name+'-'+TC_name+'-'+experiment_name+'-'+output_parameter+'-'+case_names[1]+'-'
output_fig_type = 'tif'
dpi = 75
##### Prepare output directory #####
output_dir = os.path.join(output_dir_path,output_parameter)
os.makedirs(output_dir, exist_ok=True)
##### Setting def #####
def get_value_time(nc_file_path):
# Read NetCDF file
nc_data = xr.open_dataset(nc_file_path)
# Get time
time_var = nc_data.time
# Get dimesion from time
[tmt] = time_var.shape
time = time_var[:]
# Close NetCDF file
nc_data.close()
return tmt,time
##### Prepare dimensions ########################################################
# Number of ensemble members
n_cases = len(case_names)
# Time information
time_var = (get_value_time(netcdf_dir_path+"/"+str(case_names[0])+"/"+str(case_names[0])+".pe000000.nc")[1])
nt = (get_value_time(netcdf_dir_path+"/"+str(case_names[0])+"/"+str(case_names[0])+".pe000000.nc")[0])
print(nt)
##### Main ######################################################################
ds_00 = xr.open_dataset(netcdf_dir_path+"/"+str(case_names[0])+"/"+str(case_names[0])+".pe000000.nc")
print("Now reading data is "+netcdf_dir_path+"/"+str(case_names[0])+"/"+str(case_names[0])+".pe000000.nc")
for ie in range(1,n_cases):
ds_01 = xr.open_dataset(netcdf_dir_path+"/"+str(case_names[ie])+"/"+str(case_names[ie])+".pe000000.nc")
print("Now reading data is "+netcdf_dir_path+"/"+str(case_names[ie])+"/"+str(case_names[ie])+".pe000000.nc")
output_fig_name = Model_name+'-'+TC_name+'-'+experiment_name+'-'+output_parameter+'-'+case_names[ie]+'-'
##### Tack data
# Prepare variables
lon = np.zeros((n_cases,n_track_time),dtype='float32')
lat = np.zeros((n_cases,n_track_time),dtype='float32')
xpt = np.zeros((n_cases,n_track_time))
ypt = np.zeros((n_cases,n_track_time))
print("Now reading trackdata is "+input_dir_path+"/"+input_csv_name+"-"+case_names[ie]+"-"+parameter+"-trackdata"+input_csv_type+".")
trackdata = pd.read_csv(input_dir_path+"/"+input_csv_name+"-"+case_names[ie]+"-"+parameter+"-trackdata"+input_csv_type,index_col=0)
print(trackdata)
lon[ie,:] = trackdata['lon'].tolist()
lat[ie,:] = trackdata['lat'].tolist()
if ie == 1:
time = trackdata['time'].tolist()
for jt in range(0,nt):
##### Prepare plot field
plt.rcParams["font.size"] = 18
fig = plt.figure(figsize=[10, 10], constrained_layout=True)
ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=140.0,
central_latitude=30,
standard_parallels=(30,40)))#,
ax.coastlines(resolution='50m', lw=0.5)
ax.set_title("INIT; "+time_var[0].dt.strftime("%Y-%m-%d %HUTC").values+", VALUE; "+time_var[jt].dt.strftime("%Y-%m-%d %HUTC").values)
ax.set_extent([130, 151, 20, 41], crs=ccrs.PlateCarree()) # setting draw area
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, x_inline=False, y_inline=False, linewidth=0.33, color='k',alpha=0.5)
gl.right_labels = gl.top_labels = False
gl.xlabel_style = {'rotation': 20, 'fontsize': 18}
gl.ylabel_style = {'fontsize': 18}
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlocator = mticker.FixedLocator(np.arange(130,151,5))
gl.ylocator = mticker.FixedLocator(np.arange(20,41,5))
##### Get variavles
cd_00 = ds_00.USER_CD_VALUE.isel(time=jt)#n_time)
cdr_00 = ds_00.USER_CD_RATIO.isel(time=jt)#n_time)
cd_01 = ds_01.USER_CD_VALUE.isel(time=jt)#n_time)
cdr_01 = ds_01.USER_CD_RATIO.isel(time=jt)#n_time)
##### Setting shade
slevels = np.arange(0.0002,0.0031,0.0002)
sdplot = ax.contourf(ds_00.lon, ds_00.lat,
cd_01*cdr_01-cd_00*cdr_00,
transform=ccrs.PlateCarree(),
levels = slevels,
extend='max',
cmap='YlOrRd',
zorder=0)
fig.colorbar(sdplot, ax=ax, extendrect=False, shrink=0.75, ticks=slevels)
##### Setting contour
clevels = np.max(cdr_01)
colors = ['black']
cnplot = ax.contour(ds_00.lon, ds_00.lat,
cdr_01,
levels=clevels,
transform=ccrs.PlateCarree(),
linewidths=5,
extend='both',
colors=colors,
zorder=1 )
#ax.clabel(cnplot, inline=True)
##### Tack data
ax.plot(lon[ie,:],lat[ie,:],linestyle='-',linewidth=2, color='black', label=legends[ie],transform =ccrs.PlateCarree())
ax.scatter(lon[ie,::6],lat[ie,::6], marker='o',s=30, color='black',transform=ccrs.PlateCarree())
if 'best_track_plot_swich':
print('Best Track writing True!')
if jt == 0:
ibtracs = tracks.TrackDataset(basin='all',source='ibtracs',ibtracs_mode='jtwc_neumann',catarina=True)
storm = ibtracs.get_storm(('faxai',2019))
print(storm)
best_lon_t = storm.lon
best_lat_t = storm.lat
index_6hours = pd.date_range(start=storm.time[0], end=storm.time[len(best_lon_t)-1], freq='6H')
index_1hours = pd.date_range(start=storm.time[0], end=storm.time[len(best_lon_t)-1], freq='H')
best_lon_6hourly = pd.Series(best_lon_t, index=index_6hours)
best_lat_6hourly = pd.Series(best_lat_t, index=index_6hours)
best_lon_hourly = best_lon_6hourly.reindex(index_1hours)
best_lat_hourly = best_lat_6hourly.reindex(index_1hours)
start_time_dt = pd.Timestamp(time[0])
end_time_dt = pd.Timestamp(time[-1])
best_lon_selected = best_lon_hourly.loc[start_time_dt:end_time_dt]
best_lat_selected = best_lat_hourly.loc[start_time_dt:end_time_dt]
plt.scatter(best_lon_selected, best_lat_selected, marker='o', c='gold', s=200, ec="k", linewidths=2, label='Best Track', transform=ccrs.PlateCarree())
plt.legend(loc = 3, fontsize=20) #loc:legend place ("best"1:右上, 2:左上, ...)
##### Save figure
plt.savefig(output_dir+'/'+output_fig_name+time_var[jt].dt.strftime("%Y%m%d%H").values+'.'+output_fig_type, format=output_fig_type, dpi=dpi)
plt.close()
#plt.show()
print("Finished !!!")
3.出来る図
