上次以500hpa为例绘制了天气图，这次我们来画一画地面图。
官网查看天气图
官网查看天气图的地址如下：
天气实况_天气图_中国_地面_基本天气分析 (nmc.cn)
如上图所示，地面天气图包含了比湿、温度、露点温度、海平面气压、地面风速等等要素，以下利用代码一一将其可视化。
数据资料与下载
这里我同样使用ECMWF（European Centre for Medium-Range Weather Forecasts）的再分析资料ECMWF Reanalysis v5（ERA5）数据进行可视化。选择了2024年9月21日的24h的1000hpa相对湿度（kg/kg）以及地面的露点温度（K）、温度(K)、10m风（m/s）和海平面气压（Pa）。
ERA5再分析资料的下载地址：ERA5 hourly data on pressure levels from 1940 to present (copernicus.eu)
ERA5 hourly data on single levels from 1940 to present (copernicus.eu)
代码部分
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.ticker as mticker
import xarray as xr
from scipy.ndimage import gaussian_filter,maximum_filter,minimum_filter
import pandas as pd
#-----库包的导入-----
ds1=xr.open_dataset(r"D:\Pangu-Weather-ReadyToGo\Specific-humidity-2024.9.22-24h.nc")
ds2=xr.open_dataset(r"D:\Pangu-Weather-ReadyToGo\T-DT-U-V-msl-2024.9.22-24h.nc")
#-------选取需要的数据和经纬度、时间段--------
sel_time='2024-09-22T00:00:00.000000000'
sel_lev=1000
sq=ds1.q.loc[sel_time,sel_lev,60:0,50:150]
#导入温度、露点
tem=ds2.t2m.loc[sel_time,60:0,50:150]
dtem=ds2.d2m.loc[sel_time,60:0,50:150]
u10=ds2.u10.loc[sel_time,60:0,50:150]
v10=ds2.v10.loc[sel_time,60:0,50:150]
msl=ds2.msl.loc[sel_time,60:0,50:150]
w10=np.sqrt(u10**2+v10**2)
#-----单位换算-----
msl/=100#换算成hpa
sq*=1000#换算成g/kg
tem-=273.15#摄氏度单位换算
lon,lat=tem.longitude,tem.latitude
dtem-=273.15
#导入省级行政区经纬度
df=pd.read_excel(r'全国34个省级行政区.xlsx')
lon_array,lat_array=np.array(df['经度']),np.array(df['纬度'])
#获得各个省级行政区所在经纬度格点的温度和露点
tem_value,dtem_value=[],[]
for i in range(35):
tem_sel=tem.sel(latitude=lat_array[i],longitude=lon_array[i],method='nearest')
dtem_sel=dtem.sel(latitude=lat_array[i],longitude=lon_array[i],method='nearest')
tem_value.append(float(tem_sel))
dtem_value.append(float(dtem_sel))
tem_value,dtem_value=np.array(tem_value),np.array(dtem_value)
#保留一位小数
tem_value = np.round(np.array(tem_value), 1)
dtem_value = np.round(np.array(dtem_value), 1)
#------开始绘图-------
fig =plt.figure(figsize=[10,8],dpi=500)
#设置投影方式
proj=ccrs.LambertConformal(central_latitude=90,central_longitude=104)
# proj=ccrs.PlateCarree()#其他投影方式
ax = plt.axes(projection=proj)
ax.add_feature(cfeature.LAND.with_scale('50m'), facecolor='white',edgecolor='gray') # 白色陆地
ax.add_feature(cfeature.OCEAN.with_scale('50m'), facecolor='#e7ffff') # 浅蓝色海洋
# 添加加粗边框
ax.spines['geo'].set_linewidth(1.5) # 边框加粗
#添加标识
plt.title('UTC '+sel_time[:-13]+' '+'Sea Level Pressure(hPa)',loc='left',c='blue',fontname='SimSun',fontsize=15)
# 在左下角添加标题
plt.text(0.0, -0.03, 'ECMWF Reanalysis v5(ERA5)', transform=ax.transAxes,
fontsize=12, color='blue', ha='left', va='bottom', fontname='SimHei')
#读取shp数据
with open(r"china-geospatial-data-GB2312\CN-border-La.gmt") as src:
context = src.read()
blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
borders = [np.fromstring(block,dtype=float,sep=' ')for block in blocks]
#绘制省界国界
for line in borders:
ax.plot(line[0::2],line[1::2],'-',color='#9b9f91',lw=1,transform=ccrs.Geodetic())
# 设置地图的显示范围（经纬度范围）
ax.set_extent([75, 132, 13, 55.5], crs=ccrs.PlateCarree())
# 添加网格线
gl = ax.gridlines(draw_labels=False, linewidth=0.5, color='k', alpha=0.6, linestyle='--')
# 手动设置经纬度标签
x_ticks = np.arange(80, 130 + 1, 10)
y_ticks = np.arange(10, 50 + 1, 10)
gl.xlocator = mticker.FixedLocator(x_ticks)
gl.ylocator = mticker.FixedLocator(y_ticks)
#绘制经纬度标签
for xtick in x_ticks:
ax.text(xtick, 15 - 3, f'{xtick}°E', horizontalalignment='center',
transform=ccrs.PlateCarree(), color='k', fontsize=10)
for ytick in y_ticks:
ax.text(80 - 2, ytick, f'{ytick}°N', verticalalignment='center',
transform=ccrs.PlateCarree(), color='k', fontsize=10)
msl_smooth = gaussian_filter(msl, sigma=4) # sigma 值控制平滑程度
#绘制等压线
contour_lines = ax.contour(lon,lat, msl_smooth, levels=np.arange(1000,1021,2.5), colors="blue", linewidths=1,
transform=ccrs.PlateCarree())
a1 = ax.contour(lon, lat, msl_smooth, 25, levels=[1005.0,1015.0,1025.0], linewidths=1, colors='blue',
transform=ccrs.PlateCarree())
a2 = ax.contour(lon, lat, msl_smooth, 25, levels=[1000.0,1010.0,1020.0], linewidths=1.5, colors='blue',
transform=ccrs.PlateCarree()) # 画出特定等值线，levels[588]函数画出588线
ax.clabel(a1, fmt="%2.0f", fontsize=10, colors="blue")
ax.clabel(a2, fmt="%2.0f", fontsize=10, colors="blue")
# 找到局部最大值和最小值
def plot_maxmin_points(ax, lon, lat, data, minmax, n_size, symbol, color='k',transform=None):
# outline_effect = [patheffects.withStroke(linewidth=2.5, foreground='w')]
if (minmax == 'max'):
dat = maximum_filter(data, n_size, mode='nearest')
elif (minmax == 'min'):
dat = minimum_filter(data, n_size, mode='nearest')
else:
raise ValueError('No MinMax Input')
m_xy, m_xx = np.where(dat == np.array(data))
lon, lat = np.meshgrid(lon, lat)
for i in range(len(m_xy)):
A = ax.text(lon[m_xy[i], m_xx[i]], lat[m_xy[i], m_xx[i]], symbol,
color=color, size=16,font='Times New Roman',
clip_on=True, clip_box=ax.bbox,
horizontalalignment='center', verticalalignment='center',
transform=transform)
msl_filtered=gaussian_filter(msl,sigma=3)#sigma越大越平滑
trans_crs=ccrs.PlateCarree()
#取数据局部最大最小进行标注，nsize越大越平滑
plot_maxmin_points(ax, lon, lat, msl_filtered, 'max', 35,
symbol='H', color='b', transform=trans_crs)
plot_maxmin_points(ax, lon, lat, msl_filtered, 'min', 35,
symbol='L', color='r', transform=trans_crs)
#----绘制风羽图----
uwnd,vwnd=np.array(u10),np.array(v10)
slice=16#减少数据密度 0.25*16=4°
barb=ax.barbs(lon[::slice],lat[::slice],uwnd[::slice,::slice],vwnd[::slice,::slice],barbcolor='k',
lw=0.5,length=5,barb_increments=dict(half=2,full=4,flag=20)
,transform=ccrs.PlateCarree())#风羽图设置，半个2m，整个4m，20m三角形
sq_filtered=gaussian_filter(sq,sigma=5)#sigma越大越平滑
w10_filtered=gaussian_filter(w10,sigma=5)#sigma越大越平滑
# 画 sq >= 15 的区域为绿色线
sq_contour = ax.contour(lon, lat, sq_filtered, levels=[15], colors='#0fbd0e',
linewidths=2, transform=ccrs.PlateCarree())
# 使用淡紫色栅格表示 w10 >= 10.8 的区域
w10_contourf = ax.contourf(lon, lat, w10_filtered,
levels=[10.8, np.max(w10)], colors='#c3afe4', transform=ccrs.PlateCarree())
plt.scatter(lon_array, lat_array,marker='o',color='k',
transform=ccrs.PlateCarree(),s=6.5)#打点行政区
# 在格点上方和下方添加文本
offset=0.06
for i in range(len(lat_array)):
ax.text(lon_array[i], lat_array[i] + offset, f'{tem_value[i]}',
ha='center', va='bottom', fontsize=8, color='red',transform=ccrs.PlateCarree()) # 上方文本
ax.text(lon_array[i], lat_array[i] - offset, f'{dtem_value[i]}',
ha='center', va='top', fontsize=8, color='blue',transform=ccrs.PlateCarree()) # 下方文本
# 调整边距
plt.subplots_adjust(left=0.0, right=1, top=0.9, bottom=0.1)
plt.tight_layout()
# plt.show()
outputfilename=f'中国天气图-{str(sel_time)[:-16]}-1000hpa.png'
fig.savefig(outputfilename)需要注意的是，这里导入了一个中国34个省级行政区的经纬度数据，以excel的方式存储。以下是将经纬度取出，并将对应经纬度格点的露点温度和温度都取出的代码。
#导入省级行政区经纬度
df=pd.read_excel(r'全国34个省级行政区.xlsx')
lon_array,lat_array=np.array(df['经度']),np.array(df['纬度'])
#获得各个省级行政区所在经纬度格点的温度和露点
tem_value,dtem_value=[],[]
for i in range(35):
tem_sel=tem.sel(latitude=lat_array[i],longitude=lon_array[i],method='nearest')
dtem_sel=dtem.sel(latitude=lat_array[i],longitude=lon_array[i],method='nearest')
tem_value.append(float(tem_sel))
dtem_value.append(float(dtem_sel))
tem_value,dtem_value=np.array(tem_value),np.array(dtem_value)
#保留一位小数
tem_value = np.round(np.array(tem_value), 1)
dtem_value = np.round(np.array(dtem_value), 1)
结果：
结果如下 ：
注意，以上各个参量都与国家气象局绘制的一致。
站点的填图要素包括：风向风速、温度（红色数字）和露点（蓝色数字）；天气系统符号包括：高压中心（H）、低压中心（L）、蓝色等值线为海平面气压场等压线（单位hPa）；紫色填充区为六级以上大风速区；绿色线条为等比湿线，其靠海一侧为比湿达到15g/kg的区域；
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
