Python气象绘图教程特刊（二）等值线
本节提要：站点资料的等值线、再分析资料的等值线的入门教程
一、站点资料的等值线
前几天有小伙伴在后台问站点资料和插值的问题，所以做了这一期特刊。
在气象研究领域，限制于世界的地貌和人文地理，大部分的气象原始资料是站点分布的。气象站的分布的特点是北多南少（有闲钱建设气象站的国家基本在北半球，陆地基本集中于北半球，世界人口集中于北半球），陆多海少（陆地易于永久和半永久观测站建设，海上的漂浮测站和轮船的观测不稳定）。中国的气象站密度基本与人口密度的漠河-腾冲县线吻合，表现在东多西少，中间多南北少（河北县级气象局的密度比长江以南任何一个省都高，中原地区又高于其他地区，这些牵扯到历史自然地理和人文地理）。
我前面讲这么多废话是干什么呢？主要是说明站点资料的一些使用困难，比如站点代表性不足（青藏高原西部一个站需要代表十几万平方公里）。下面才是具体的站点绘制等值线教程。
首先，需要大家回忆一下天气学分析基础和典型天气过程分析这两门课，在绘制500百帕高空场和地面形势图的位势线和等压线时，你是怎么连接两个站点的。
下面举个栗子：
你将会如何连接这两个线呢，相信不会直接连接两个1000百帕，而是会连接一根曲线，曲线向圆点方向凹陷。那么，你这用到了哪些方法而不自知呢？
首先，你参考了1002.5百帕这个点，这导致了你决定连线的凹凸性。下面，通过一个图理解什么是站点数据的网格化和插值。
这就是你在进行天气分析时心里默认的连线方式了。现在每0.5单位分隔，然后假定的在其中某些位置存在1000百帕的点，然后将这些点连接起来。
其实，在你的脑子里，你已经进行了网格化和插值的操作，这些红色的数字就是你人脑插值法进行的插值，人脑插值法是你脑子里一种感性的插值法，与我们的使用的径向基函数插值、克里金插值等方法相比，误差较大比较随意。比如上面这幅图，单纯从连线角度来说，以直线连接是最简洁的，但是你一般不会这样做，等压线很少有笔直的情况。你会画一个弧线出来，但是究竟曲率多大，就因人而异了，这是人脑插值法的缺点。
使用径向基函数插值（黑线）和我人脑插值（红线）的比较：
实例一：恩施州累年一月平均气温
from scipy.interpolate import Rbf#引入径向基函数
filename=r'C:\Users\lenovo\Desktop\ex1.xlsx'#数据文件地址
df=pd.read_excel(filename)#读取文件
lon=df['lon']#读取站点经度
lat=df['lat']#读取站点纬度
tem=df['temp']#读取站点气温
olon=np.linspace(108,111,30)#设置网格经度
olat=np.linspace(29,32,30)#设置网格纬度
olon,olat=np.meshgrid(olon,olat)#网格化
func=Rbf(lon,lat,tem,function='linear')#定义径向基函数插值
tem_new=func(olon,olat)#获得插值后的网格气温
cs= ax.contourf(olon,olat,tem_new,levels=np.arange(6,9,0.3),extend='both')#画图
clip=maskout.shp2clip(cs, ax,r'E:\dijishi\cn_province.shp' ,420000) #白化上面我们给了两个变量olon,olat进行网格化操作。经度从108到111，分为30格，每格0.1°；纬度从29到32，分为30，每格0.1°，所以我们获得的新网格精度为0.1°*0.1°。
实例二：恩施州去年累计降水量
from scipy.interpolate import Rbf#引入径向基函数
filename=r'C:\Users\lenovo\Desktop\累年降水数据.xlsx'#数据文件地址
df=pd.read_excel(filename)#读取文件
lon=df['lon']#读取站点经度
lat=df['lat']#读取站点纬度
rain=df['precipitation']#读取站点累计年降水量
olon=np.linspace(108,111,30)#设置网格经度
olat=np.linspace(29,32,30)#设置网格纬度
olon,olat=np.meshgrid(olon,olat)#网格化
func=Rbf(lon,lat,rain,function='linear')#定义插值函数
rain_new=func(olon,olat)#获得插值后的网格累计降水量
cs= ax.contourf(olon,olat,rain_new,levels=np.arange(900,2000,100),cmap='GnBu',extend='both')#画图
clip=maskout.shp2clip(cs, ax,r'E:\dijishi\cn_province.shp' ,420000) #白化二、再分析资料绘制等值线
这里提供了两个再分析资料，一个是sst，一个是地表温度,都是nooa的nc资料。再分析资料其实就是站点资料已经网格化、插值后的资料了，所以只需要直接读取然后画图就可以了。
实例三：再分析地表温度
from netCDF4 import Dataset
nc_file= Dataset(r'E:\\aaaa\air.mon.mean.nc')
lon=nc_file.variables['lon'][:]#读取经度
lat=nc_file.variables['lat'][:]#读取纬度
temp=nc_file.variables['air'][11][:][:]-273.15#读取11时刻气温，-273.15是转换气温为摄氏度
ax.contourf(lon, lat, temp,levels=np.arange(-40,40),cmap='RdBu_r')实例四：再分析海温场
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from netCDF4 import Dataset
from cartopy.util import add_cyclic_point
nc_file=Dataset('E:\\aaaa\sst.mon.mean.nc')
lon=nc_file.variables['lon'][:]#读取经度
lat=nc_file.variables['lat'][:]#读取纬度
sst=nc_file.variables['sst'][0][:][:]#读取0时刻海温
fig = plt.figure(figsize=(8,6),dpi=500)
ax=fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
ax.coastlines()#添加海岸线
ax.add_feature(cfeat.LAND)#添加陆地
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.7, color='k', alpha=1, linestyle='--')
gl.xlabels_top = False # 关闭顶端的经纬度标签
gl.ylabels_right = False # 关闭右侧的经纬度标签
gl.xformatter = LONGITUDE_FORMATTER # x轴设为经度的格式
gl.yformatter = LATITUDE_FORMATTER # y轴设为纬度的格式
gl.xlocator = mticker.FixedLocator(np.arange(-180, 180.1, 30))
gl.ylocator = mticker.FixedLocator(np.arange(-90,90.1, 30))
gl.xlines=False#关闭内部经纬线
gl.ylines=False#关闭内部经纬线
#################以下步骤添加数据循环，防止白条##################
cycle_sst, cycle_lon = add_cyclic_point(sst, coord=lon)
cycle_LON, cycle_LAT = np.meshgrid(cycle_lon, lat)
ax.contourf(cycle_LON, cycle_LAT, cycle_sst,levels=np.arange(-10,35),cmap='RdBu_r')
plt.savefig('海温场',bbox_inches='tight')
plt.show()在绘制海温场时，出现了网格化命令np.meshgrid，请注意，这是为了添加数据循环防止白条，而不是前面对恩施州站点资料进行网格化。
再分析资料链接：
链接：https://pan.baidu.com/s/1wCj95tJSSrF9_G3j_vd6hA
提取码：cyar
下期可能继续讲legend和colorbar。
腾讯云开发者
