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.出来る図
