台风
最近开学了，有些跨专业的师弟师妹不太会绘制气象的二维图，今天以最近台风"摩羯"为例子。
这次使用ERA5再分析数据为例。
Python
plot_maxmin_points寻找局地最大值以及最小值点
from datetime import datetime
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
import matplotlib.pyplot as plt
import matplotlib.colors as colors
from matplotlib import patheffects
import metpy.calc as mpcalc
from metpy.units import units
import numpy as np
from scipy.ndimage import gaussian_filter, maximum_filter, minimum_filter
import xarray as xr
# Plotting Max and Min data points
def plot_maxmin_points(ax, lon, lat, data, extrema, nsize, symbol, color='k',
plotValue=True, transform=None):
outline_effect = [patheffects.withStroke(linewidth=2.5, foreground='w')]
if (extrema == 'max'):
data_ext = maximum_filter(data, nsize, mode='nearest')
elif (extrema == 'min'):
data_ext = minimum_filter(data, nsize, mode='nearest')
else:
raise ValueError('Value for hilo must be either max or min')
mxy, mxx = np.where(data_ext == np.array(data))
lon, lat = np.meshgrid(lon, lat)
for i in range(len(mxy)):
A = ax.text(lon[mxy[i], mxx[i]], lat[mxy[i], mxx[i]], symbol, color=color, size=24,
clip_on=True, clip_box=ax.bbox, horizontalalignment='center', verticalalignment='center',
transform=transform)
A.set_path_effects(outline_effect)
a = np.array(data[mxy[i], mxx[i]])
a_trunc = np.trunc(a)
B = ax.text(lon[mxy[i], mxx[i]], lat[mxy[i], mxx[i]],
'\n' + str(a_trunc),
color=color, size=12,
clip_on=True,
clip_box=ax.bbox, fontweight='bold',
horizontalalignment='center', verticalalignment='top', transform=transform)
B.set_path_effects(outline_effect)
读取ERA5数据，对单位进行转换。
# Open ERA5 dataset; change name accordingly to yours
ds = xr.open_dataset(r'G:\EdgeDownLoad/era5_700hPa_20240904_0908.nc').sel(latitude = slice(60, 0), longitude = slice(80, 180)).metpy.parse_cf()
# Open another ERA5 dataset (optional); change name accordingly to yours
ds1 = xr.open_dataset(r'G:\EdgeDownLoad/era5_single_20240904_0908.nc').sel(latitude = slice(60, 0), longitude = slice(80, 180)).metpy.parse_cf()
# Grab lat/lon values
lats = ds.latitude
lons = ds.longitude
# Select and grab data
w = ds['w']
geop = ds['z']
geop_w_units = geop * units('m ** 2 / s ** 2')
z = mpcalc.geopotential_to_height(geop_w_units)
rain = ds1['mtpr']
# Select and grab 700-hPa geopotential heights, wind components, and vertical velocity
# Then, smooth with gaussian_filter
w_700 = gaussian_filter(w.sel(pressure_level=700).data[0], sigma=1.485) * units('Pa/s')
w_700ub = w_700.to('microbar/second')
hght_700 = gaussian_filter(z.sel(pressure_level=700).data[0], sigma=3.0) * units('m')
hght_700dm = hght_700.to(units.decameter)
# Add precipitation and convert to mm/hr. Smooth with gaussian_filter
precip = gaussian_filter(rain.data[0], sigma=0.5) * units('kg/m^2/s')
mmhr1 = precip * (1 * units('m^3')) / (1000 * units('kg'))
mmhr2 = mmhr1 * (1000 * units('mm')) / (1 * units.meters)
mmhr = mmhr2.to('mm/hr') #Conversion
# Create a clean datetime object for plotting based on time of Geopotential heights
vtime = datetime.strptime(str(ds1.valid_time.data[1].astype('datetime64[ms]')),
'%Y-%m-%dT%H:%M:%S.%f')
进行绘图。
# Set up the projection that will be used for plotting
mapcrs = ccrs.PlateCarree(central_longitude=180)
# Set up the projection of the data; if lat/lon then PlateCarree is what you want
datacrs = ccrs.PlateCarree(central_longitude=0)
# Start the figure and create plot axes with proper projection
fig = plt.figure(1, figsize=(14, 12), dpi = 300)
ax = plt.subplot(111, projection=datacrs)
ax.set_extent([100, 140, 5, 30], datacrs)
# Add geopolitical boundaries for map reference
land = cfeature.NaturalEarthFeature('physical', 'land', '50m')
ax.add_feature(land, color='lightgray', zorder=-1)
ax.coastlines()
# Highlight boundary on VT
clevs_700_vtt = np.arange(-10, 8, 2)
vt = ax.contour(lons, lats, w_700ub, clevs_700_vtt, cmap='seismic_r',
transform=datacrs, alpha=0.9, linestyles='--')
# Plot 700-hPa Geopotential Heights in meters
clevs_700_hght = np.arange(305, 316, 1)
cs = ax.contour(lons, lats, hght_700dm, clevs_700_hght, colors='blue',
transform=datacrs)
precip_colors = [
"#bde9bf", # 0.01 - 0.02 mm hr-1 1
"#adddb0", # 0.02 - 0.03 mm hr-1 2
"#9ed0a0", # 0.03 - 0.04 mm hr-1 3
"#8ec491", # 0.04 - 0.05 mm hr-1 4
"#7fb882", # 0.05 - 0.06 mm hr-1 5
"#70ac74", # 0.06 - 0.07 mm hr-1 6
"#60a065", # 0.07 - 0.08 mm hr-1 7
"#519457", # 0.08 - 0.09 mm hr-1 8
"#418849", # 0.90 - 1.00 mm hr-1 9
"#307c3c", # 1.00 - 1.20 mm hr-1 10
"#1c712e", # 1.20 - 1.40 mm hr-1 11
"#f7f370", # 1.40 - 1.60 mm hr-1 12
"#fbdf65", # 1.60 - 1.80 mm hr-1 13
"#fecb5a", # 1.80 - 2.00 mm hr-1 14
"#ffb650", # 2 - 3 mm hr-1 15
"#ffa146", # 3 - 4 mm hr-1 16
"#ff8b3c", # 4 - 5 mm hr-1 17
"#ff8b3c", # 5 - 6 mm hr-1 18
"#ff8b5c", # 6 - 7 mm hr-1 19
"#ff8b5c" # 7 - 8 mm hr-1 19
]
precip_colormap = colors.ListedColormap(precip_colors)
clev_precip = np.concatenate((np.arange(0.1, 1, .1), np.arange(1, 2, .2), np.arange(2, 9, 1)))
norm = colors.BoundaryNorm(clev_precip, 20)
# Plot contours of Average Total Precipitation Rate
cf = ax.contourf(lons, lats, mmhr, clev_precip, norm=norm, cmap=precip_colormap, transform=datacrs)
cbar = plt.colorbar(cf, ticks=clev_precip, orientation='vertical', pad=0.01, aspect=35, shrink=0.45)
cbar.ax.tick_params(which='major', length=10, width=1.5, top=True, right=True, direction='in')
# Plot Highs and Lows (Min/Max) values of Geopotential Height
plot_maxmin_points(ax, lons, lats, hght_700dm, 'max', 40, symbol='H', color='b', transform=datacrs)
plot_maxmin_points(ax, lons, lats, hght_700dm, 'min', 40, symbol='L', color='r', transform=datacrs)
# Set X and Y-ticks for Latitude and Longitude Coordinates
lat_formatter = cticker.LatitudeFormatter()
long_formatter = cticker.LongitudeFormatter()
ax.yaxis.set_major_formatter(lat_formatter)
ax.xaxis.set_major_formatter(long_formatter)
ax.set_xticks(np.arange(100, 141, 10)) #Hide end points of longitude tickmarks
ax.set_yticks(np.arange(5, 30 + 1, 5)) #Hide end points of latitude tickmarks
ax.minorticks_on()
ax.tick_params(which='major', length=10, width=1.5, top=True, right=True, direction='in')
ax.tick_params(which='minor', length=5, width=1, top=True, right=True, direction='in')
# Adjust image and show
plt.subplots_adjust(bottom=0, top=1)
plt.rcParams['axes.linewidth'] = 2
plt.rcParams['font.size'] = 10
plt.rcParams['figure.dpi'] = 500
plt.rcParams['figure.facecolor'] = 'white'
plt.rcParams['mathtext.fontset'] = 'dejavuserif'
plt.rcParams['font.family'] = 'STIXGeneral'
# Make some nice titles for the plot (one right, one left)
plt.title('ERA5: Z$_{g}$(h)$_{700}$ (dam), $\omega$$_{700}$ (Pa s$^{-1}$), Mean Precip. (mm hr$^{-1}$)', loc='left', fontsize = 15)
plt.title('Time: {} UTC'.format(vtime), loc='right', fontsize = 15)
plt.show()
完整代码
from datetime import datetime
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
import matplotlib.pyplot as plt
import matplotlib.colors as colors
from matplotlib import patheffects
import metpy.calc as mpcalc
from metpy.units import units
import numpy as np
from scipy.ndimage import gaussian_filter, maximum_filter, minimum_filter
import xarray as xr
# Plotting Max and Min data points
def plot_maxmin_points(ax, lon, lat, data, extrema, nsize, symbol, color='k',
plotValue=True, transform=None):
outline_effect = [patheffects.withStroke(linewidth=2.5, foreground='w')]
if (extrema == 'max'):
data_ext = maximum_filter(data, nsize, mode='nearest')
elif (extrema == 'min'):
data_ext = minimum_filter(data, nsize, mode='nearest')
else:
raise ValueError('Value for hilo must be either max or min')
mxy, mxx = np.where(data_ext == np.array(data))
lon, lat = np.meshgrid(lon, lat)
for i in range(len(mxy)):
A = ax.text(lon[mxy[i], mxx[i]], lat[mxy[i], mxx[i]], symbol, color=color, size=24,
clip_on=True, clip_box=ax.bbox, horizontalalignment='center', verticalalignment='center',
transform=transform)
A.set_path_effects(outline_effect)
a = np.array(data[mxy[i], mxx[i]])
a_trunc = np.trunc(a)
B = ax.text(lon[mxy[i], mxx[i]], lat[mxy[i], mxx[i]],
'\n' + str(a_trunc),
color=color, size=12,
clip_on=True,
clip_box=ax.bbox, fontweight='bold',
horizontalalignment='center', verticalalignment='top', transform=transform)
B.set_path_effects(outline_effect)
# Open ERA5 dataset; change name accordingly to yours
ds = xr.open_dataset(r'G:\EdgeDownLoad/era5_700hPa_20240904_0908.nc').sel(latitude = slice(60, 0), longitude = slice(80, 180)).metpy.parse_cf()
# Open another ERA5 dataset (optional); change name accordingly to yours
ds1 = xr.open_dataset(r'G:\EdgeDownLoad/era5_single_20240904_0908.nc').sel(latitude = slice(60, 0), longitude = slice(80, 180)).metpy.parse_cf()
# Grab lat/lon values
lats = ds.latitude
lons = ds.longitude
# Select and grab data
w = ds['w']
geop = ds['z']
geop_w_units = geop * units('m ** 2 / s ** 2')
z = mpcalc.geopotential_to_height(geop_w_units)
rain = ds1['mtpr']
# Select and grab 700-hPa geopotential heights, wind components, and vertical velocity
# Then, smooth with gaussian_filter
w_700 = gaussian_filter(w.sel(pressure_level=700).data[0], sigma=1.485) * units('Pa/s')
w_700ub = w_700.to('microbar/second')
hght_700 = gaussian_filter(z.sel(pressure_level=700).data[0], sigma=3.0) * units('m')
hght_700dm = hght_700.to(units.decameter)
# Add precipitation and convert to mm/hr. Smooth with gaussian_filter
precip = gaussian_filter(rain.data[0], sigma=0.5) * units('kg/m^2/s')
mmhr1 = precip * (1 * units('m^3')) / (1000 * units('kg'))
mmhr2 = mmhr1 * (1000 * units('mm')) / (1 * units.meters)
mmhr = mmhr2.to('mm/hr') #Conversion
# Create a clean datetime object for plotting based on time of Geopotential heights
vtime = datetime.strptime(str(ds1.valid_time.data[1].astype('datetime64[ms]')),
'%Y-%m-%dT%H:%M:%S.%f')
# Set up the projection that will be used for plotting
mapcrs = ccrs.PlateCarree(central_longitude=180)
# Set up the projection of the data; if lat/lon then PlateCarree is what you want
datacrs = ccrs.PlateCarree(central_longitude=0)
# Start the figure and create plot axes with proper projection
fig = plt.figure(1, figsize=(14, 12), dpi = 300)
ax = plt.subplot(111, projection=datacrs)
ax.set_extent([100, 140, 5, 30], datacrs)
# Add geopolitical boundaries for map reference
land = cfeature.NaturalEarthFeature('physical', 'land', '50m')
ax.add_feature(land, color='lightgray', zorder=-1)
ax.coastlines()
# Highlight boundary on VT
clevs_700_vtt = np.arange(-10, 8, 2)
vt = ax.contour(lons, lats, w_700ub, clevs_700_vtt, cmap='seismic_r',
transform=datacrs, alpha=0.9, linestyles='--')
# Plot 700-hPa Geopotential Heights in meters
clevs_700_hght = np.arange(305, 316, 1)
cs = ax.contour(lons, lats, hght_700dm, clevs_700_hght, colors='blue',
transform=datacrs)
precip_colors = [
"#bde9bf", # 0.01 - 0.02 mm hr-1 1
"#adddb0", # 0.02 - 0.03 mm hr-1 2
"#9ed0a0", # 0.03 - 0.04 mm hr-1 3
"#8ec491", # 0.04 - 0.05 mm hr-1 4
"#7fb882", # 0.05 - 0.06 mm hr-1 5
"#70ac74", # 0.06 - 0.07 mm hr-1 6
"#60a065", # 0.07 - 0.08 mm hr-1 7
"#519457", # 0.08 - 0.09 mm hr-1 8
"#418849", # 0.90 - 1.00 mm hr-1 9
"#307c3c", # 1.00 - 1.20 mm hr-1 10
"#1c712e", # 1.20 - 1.40 mm hr-1 11
"#f7f370", # 1.40 - 1.60 mm hr-1 12
"#fbdf65", # 1.60 - 1.80 mm hr-1 13
"#fecb5a", # 1.80 - 2.00 mm hr-1 14
"#ffb650", # 2 - 3 mm hr-1 15
"#ffa146", # 3 - 4 mm hr-1 16
"#ff8b3c", # 4 - 5 mm hr-1 17
"#ff8b3c", # 5 - 6 mm hr-1 18
"#ff8b5c", # 6 - 7 mm hr-1 19
"#ff8b5c" # 7 - 8 mm hr-1 19
]
precip_colormap = colors.ListedColormap(precip_colors)
clev_precip = np.concatenate((np.arange(0.1, 1, .1), np.arange(1, 2, .2), np.arange(2, 9, 1)))
norm = colors.BoundaryNorm(clev_precip, 20)
# Plot contours of Average Total Precipitation Rate
cf = ax.contourf(lons, lats, mmhr, clev_precip, norm=norm, cmap=precip_colormap, transform=datacrs)
cbar = plt.colorbar(cf, ticks=clev_precip, orientation='vertical', pad=0.01, aspect=35, shrink=0.45)
cbar.ax.tick_params(which='major', length=10, width=1.5, top=True, right=True, direction='in')
# Plot Highs and Lows (Min/Max) values of Geopotential Height
plot_maxmin_points(ax, lons, lats, hght_700dm, 'max', 40, symbol='H', color='b', transform=datacrs)
plot_maxmin_points(ax, lons, lats, hght_700dm, 'min', 40, symbol='L', color='r', transform=datacrs)
# Set X and Y-ticks for Latitude and Longitude Coordinates
lat_formatter = cticker.LatitudeFormatter()
long_formatter = cticker.LongitudeFormatter()
ax.yaxis.set_major_formatter(lat_formatter)
ax.xaxis.set_major_formatter(long_formatter)
ax.set_xticks(np.arange(100, 141, 10)) #Hide end points of longitude tickmarks
ax.set_yticks(np.arange(5, 30 + 1, 5)) #Hide end points of latitude tickmarks
ax.minorticks_on()
ax.tick_params(which='major', length=10, width=1.5, top=True, right=True, direction='in')
ax.tick_params(which='minor', length=5, width=1, top=True, right=True, direction='in')
# Adjust image and show
plt.subplots_adjust(bottom=0, top=1)
plt.rcParams['axes.linewidth'] = 2
plt.rcParams['font.size'] = 10
plt.rcParams['figure.dpi'] = 500
plt.rcParams['figure.facecolor'] = 'white'
plt.rcParams['mathtext.fontset'] = 'dejavuserif'
plt.rcParams['font.family'] = 'STIXGeneral'
# Make some nice titles for the plot (one right, one left)
plt.title('ERA5: Z$_{g}$(h)$_{700}$ (dam), $\omega$$_{700}$ (Pa s$^{-1}$), Mean Precip. (mm hr$^{-1}$)', loc='left', fontsize = 15)
plt.title('Time: {} UTC'.format(vtime), loc='right', fontsize = 15)
plt.show()
往期回顾
Python | MJO | 位相图 Python | 半球间经度转换 Python | 大气科学 | 偏相关 37 个 CMIP6 模型对流参数化方案 Python | 解决 cmaps 库报错的问题 Python | 批量下载 NCEP2 再分析数据 Python | 北大西洋涛动 | NAO 指数 | EOF
