import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from scipy.ndimage import gaussian_filter # 地形平滑(可选)
def plot_simple_contour(wrfout_file, target_lon, target_lat, radius_km=55,
contour_interval=20, smooth_terrain=True):
# 1. 读取WRF数据
nc = Dataset(wrfout_file, 'r')
try:
lats = nc.variables['XLAT'][:]
lons = nc.variables['XLONG'][:]
hgt = nc.variables['HGT'][:]
# 处理时间维度(取第一时次)
if len(lats.shape) == 3:
lats = lats[0, :, :]
lons = lons[0, :, :]
hgt = hgt[0, :, :]
# 2. 智能截取目标网格(按公里数动态计算范围)
deg2km = 111 # 1度≈111公里
n = max(10, int(radius_km / deg2km * max(lats.shape))) # 保证至少10个网格
# 找目标点在网格中的索引
lat_diff = np.abs(lats - target_lat)
lon_diff = np.abs(lons - target_lon)
min_idx = np.argmin(lat_diff + lon_diff)
sy, sx = np.unravel_index(min_idx, lats.shape)
# 边界保护(防止索引越界)
sy_min = max(0, sy - n)
sy_max = min(lats.shape[0] - 1, sy + n)
sx_min = max(0, sx - n)
sx_max = min(lats.shape[1] - 1, sx + n)
# 截取目标区域
target_lats = lats[sy_min:sy_max, sx_min:sx_max]
target_lons = lons[sy_min:sy_max, sx_min:sx_max]
target_hgt = hgt[sy_min:sy_max, sx_min:sx_max]
# 数据校验
if target_hgt.size == 0:
raise ValueError("目标区域无数据!请检查经纬度或半径设置")
# 3. 地形平滑(可选)
if smooth_terrain:
target_hgt = gaussian_filter(target_hgt, sigma=1.2) # 轻度平滑
# 4. 创建地图投影和绘图
# 使用PlateCarree投影(经纬度投影)
projection = ccrs.PlateCarree()
# 创建图形和轴
fig = plt.figure(figsize=(10, 8), dpi=150)
ax = plt.subplot(111, projection=projection)
# 5. 添加地图底图
# 添加陆地、海洋、湖泊和海岸线特征
ax.add_feature(cfeature.LAND, facecolor='#f5f5dc', edgecolor='none') # 浅米色陆地
ax.add_feature(cfeature.OCEAN, facecolor='#e6f2ff', edgecolor='none') # 浅蓝色海洋
ax.add_feature(cfeature.LAKES, facecolor='#e6f2ff', edgecolor='none') # 湖泊
ax.add_feature(cfeature.COASTLINE, edgecolor='#666666', linewidth=0.8) # 海岸线
ax.add_feature(cfeature.BORDERS, linestyle=':', edgecolor='#999999', linewidth=0.8) # 国界线
# 6. 设置地图范围
# 计算目标区域的边界
lon_min, lon_max = target_lons.min(), target_lons.max()
lat_min, lat_max = target_lats.min(), target_lats.max()
# 扩展10%的范围以获得更好的上下文
lon_center = (lon_min + lon_max) / 2
lat_center = (lat_min + lat_max) / 2
lon_extent = (lon_max - lon_min) * 0.1
lat_extent = (lat_max - lat_min) * 0.1
# 设置地图范围
ax.set_extent([
lon_min - lon_extent,
lon_max + lon_extent,
lat_min - lat_extent,
lat_max + lat_extent
], crs=projection)
# 7. 添加网格线和标签
gl = ax.gridlines(crs=projection, draw_labels=True,
linewidth=0.5, color='gray', alpha=0.5, linestyle='--')
gl.top_labels = False
gl.right_labels = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'size': 9, 'color': 'black'}
gl.ylabel_style = {'size': 9, 'color': 'black'}
# 8. 绘制等高线(在地图底图上)
# 计算合理等高线范围(用百分位数避免极端值)
hgt_min = np.percentile(target_hgt, 5)
hgt_max = np.percentile(target_hgt, 95)
levels = np.arange(hgt_min, hgt_max + contour_interval, contour_interval)
# 等高线(深棕色线条)
contour = ax.contour(
target_lons, target_lats, target_hgt,
levels=levels,
colors='#8B4513', # 深棕色
linewidths=1.5,
transform=projection,
zorder=2
)
# 9. 等高线标注(深棕色)
ax.clabel(
contour,
contour.levels, # 标注所有层级
inline=True,
fontsize=9,
fmt="%d",
colors='#8B4513', # 深棕色
zorder=3
)
# 10. 目标点标注(红色星形)
ax.plot(
target_lon, target_lat,
'r*', markersize=12,
transform=projection,
zorder=4,
label=f'目标点 ({target_lon:.2f}°E, {target_lat:.2f}°N)'
)
ax.legend(
loc='upper right',
fontsize=10,
frameon=True,
framealpha=0.8,
edgecolor='gray'
)
# 在右上角添加指北针
compass_x = lon_max - 0.1 * (lon_max - lon_min)
compass_y = lat_max - 0.1 * (lat_max - lat_min)
ax.text(compass_x, compass_y, 'N',
ha='center', va='center', fontsize=12,
fontweight='bold', transform=projection, zorder=5)
ax.plot([compass_x, compass_x], [compass_y, compass_y + 0.2],
color='black', linewidth=1, transform=projection, zorder=5)
# 13. 调整布局并显示
plt.tight_layout()
plt.subplots_adjust(top=0.92) # 为标题留出空间
plt.show()
finally:
nc.close()
# ------------------- 调用示例 -------------------
if __name__ == "__main__":
wrfout_file = r"C:\1\wrfout_d03_2024-01-21_18" # 替换为实际路径
target_lon = 92.070 # 目标经度
target_lat = 31.455 # 目标纬度
plot_simple_contour(
wrfout_file,
target_lon,
target_lat,
radius_km=25, # 范围大小
contour_interval=20, # 等高线间隔(可调整)
smooth_terrain=True # 平滑地形(让线条更自然)
)
修改代码,经纬度统一保留三位小数
最新发布