WRF后处理 | contour图的绘制

绘制流程：

绘制WRF（Weather Research and Forecasting）模型的等值线（contour）图需要使用WRF数据和wrf-python库、Matplotlib库或其他绘图库。以下是一般步骤来绘制WRF模型的等值线图：

1. 获取WRF数据：首先，需要获取WRF模型的输出数据文件。通常，WRF输出文件以NetCDF格式存储，可以使用Python的netCDF4库来读取数据。
2. 提取所需的变量：可以使用getvarWRF输出文件中提取您希望绘制的变量，例如温度、降水量、风速等。
3. 创建网格：如果您的WRF输出数据是经纬度坐标的，您可以使用latlon_coords函数来创建网格，以便在等值线图中绘制。
4. 绘制等值线图：使用Matplotlib的contour函数或contourf函数来绘制等值线图。contour函数绘制轮廓线，contourf函数绘制填充轮廓线，具体取决于您的需求。

Basemap 和 Cartopy:

在Python中，要绘制WRF（Weather Research and Forecasting）模型的地图，还有两个常用的库：Basemap和Cartopy。这两个库都可以用于制作地图，但它们具有不同的特点和用法。

1. Basemap：
   • Basemap是Matplotlib的一部分，它提供了一种相对简单的方式来绘制地图和地理数据。
   • Basemap支持多种地图投影。
   • Basemap可以轻松添加海岸线、国界、河流等地理特征。
2. Cartopy：
   • Cartopy是一个较新的库，它提供了更灵活的地图绘制选项和更多地图投影。
   • Cartopy的设计更现代，支持自定义地图特征和坐标系。
   • Cartopy还提供了更好的地理数据集成，以便更容易地处理地理数据。

总的来说，如果您需要更高级的地图制作功能，或者需要更多地理数据处理能力，那么Cartopy可能是更好的选择。但如果您只需要简单的地图，Basemap可能更容易上手。根据您的需求和个人偏好，您可以选择使用其中一个库来制作WRF模型的地图可视化。

示例：

Basemap：

from netCDF4 import Dataset
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from mpl_toolkits.basemap import Basemap
from pylab import mpl  
import numpy as np
from wrf import to_np, getvar, smooth2d, get_basemap, latlon_coords
import numpy as np
mpl.rcParams['font.sans-serif'] = ['SimHei']

#choose the area
# def select_and_return_index_2d(arr2d, min_val, max_val):
#   selected = arr2d[(arr2d >= min_val) & (arr2d <= max_val)]
#   index = np.where((arr2d >= min_val) & (arr2d <= max_val))
#   return selected, index

# Open the NetCDF file
ncfile = Dataset("G:\wrfout_d01_2023-07-01_06_00_00")

# Get the sea level pressure
slp = getvar(ncfile, "slp")

# Smooth the sea level pressure since it tends to be noisy near the
# mountains
smooth_slp = smooth2d(slp, 3, cenweight=4)

# Get the latitude and longitude points
lats, lons = latlon_coords(slp)

# Get the basemap object
bm = get_basemap(slp,llcrnrlat=18, urcrnrlat=23,llcrnrlon=114, urcrnrlon=119,resolution='h')

levels=np.linspace(1040,1056,10)
# Create a figure
fig = plt.figure(figsize=(12,9))
ax = fig.add_subplot(111)

# Add geographic outlines
bm.drawcoastlines(linewidth=0.25)
bm.drawstates(linewidth=0.25)
bm.drawcountries(linewidth=0.25)
# bm.drawlsmask(land_color='gray',ocean_color='white')  #陆地，海洋的颜色

# Convert the lats and lons to x and y.  Make sure you convert the lats and
# lons to numpy arrays via to_np, or basemap crashes with an undefined
# RuntimeError.
x, y = bm(to_np(lons), to_np(lats))

# Draw the contours and filled contours
bm.contour(x, y, to_np(smooth_slp), levels=levels,colors="black")
bm.contourf(x, y, to_np(smooth_slp), levels=levels, cmap=get_cmap("jet"))

bm.drawparallels( np.arange(18.,24.,1.),   labels=[1,0,0,0], color='k', fontsize=12 )
bm.drawmeridians( np.arange(114.,120.,1.), labels=[0,0,0,1], color='k' , fontsize=12 )

# Add a color bar
plt.colorbar(shrink=1)

plt.title("Sea Level Pressure (hPa)")

plt.show()

Cartopy：

from netCDF4 import Dataset
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
import cartopy.feature as cfeature
import cmaps
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
import numpy as np
import cartopy.crs as ccrs

from wrf import (to_np, getvar, smooth2d, get_cartopy, cartopy_xlim,
                 cartopy_ylim, latlon_coords)

# Open the NetCDF file
ncfile = Dataset("G:\wrfout_d01_2023-07-01_06_00_00")

# Get the sea level pressure
slp = getvar(ncfile, "slp")

# Smooth the sea level pressure since it tends to be noisy near the
# mountains
smooth_slp = smooth2d(slp, 3, cenweight=4)

# Get the latitude and longitude points
lats, lons = latlon_coords(slp)

# Get the cartopy mapping object
cart_proj = get_cartopy(slp)
levels=np.linspace(1040,1056,10)

# Create a figure
fig = plt.figure(figsize=(12,6))
# Set the GeoAxes to the projection used by WRF
ax = plt.axes(projection=cart_proj)


# Make the contour outlines and filled contours for the smoothed sea level
# pressure.
plt.contour(to_np(lons), to_np(lats), to_np(smooth_slp), levels=levels, colors="black",
            transform=crs.PlateCarree())
mesh=plt.contourf(to_np(lons), to_np(lats), to_np(smooth_slp), levels=levels,
             transform=crs.PlateCarree(),
             cmap=cmaps.MPL_jet)

ax.add_feature(cfeature.LAND)#添加陆地
ax.add_feature(cfeature.COASTLINE,lw = 0.3)#添加海岸线
# ax.add_feature(cfeature.OCEAN)#添加海洋

# Add a color bar
colorbar=plt.colorbar(ax=ax, shrink=.98)

# 获取颜色条相关联的轴对象
# cax = colorbar.ax

# # 设置颜色条刻度标签的字体大小
# cax.tick_params(labelsize=5)  # 将字体大小设置为 1
# Set the map bounds
ax.set_xlim(cartopy_xlim(smooth_slp))
ax.set_ylim(cartopy_ylim(smooth_slp))

mesh = ax.gridlines(draw_labels=True, linestyle='--', linewidth=0.6, alpha=0.5, x_inline=False, y_inline=False, color='k')
mesh.top_labels=False
mesh.right_labels=False
mesh.xformatter = LONGITUDE_FORMATTER
mesh.yformatter = LATITUDE_FORMATTER
mesh.xlocator = mticker.FixedLocator(np.arange(114,120,1))
mesh.ylocator = mticker.FixedLocator(np.arange(18,24,1))
mesh.xlabel_style={'size':12}
mesh.ylabel_style={'size':12}

ax.set_extent([114,119,18,23],crs=ccrs.PlateCarree()) # 小范围

plt.title("Sea Level Pressure (hPa)")

plt.show()