python绘制降水色斑图

前言

本来想绘制如下的色斑图，但一开始不知到它叫这个名字，从等值线开始查起，发现等值线绘制是个比较大且难的问题，会出现等值点计算、等值点追踪、等值线裁剪等一些列的部分所组成，预想到最后还可能会出现效率问题，然，现今已有很多成熟和软件已集成（实现）了该功能，故在查找方法（方式）过程中小伙伴发现了它原来的真实名字，好了，废话到此结束。

一、色斑图绘制

1. 加载数据

做加载的数据包括，经度、维度、降水值三类数据，我把它们存在了一个csv文件中了，由于数据提前已进行了插值处理，因此这里不再进行插值计算。

def get_origin_data(data_path, num):"""获取经纬度数据， 雨强数据:param data_path: 数据路径:param num: 格网数:return:"""df = pd.read_csv(data_path)     # 读取数据df.columns = ['lon', 'lat', 'rain']    # 列命名olon = np.array(df['lon']).reshape(num, num)    # 经度olat = np.array(df['lat']).reshape(num, num)    # 纬度rain = np.array(df['rain']).reshape(num, num)   # 雨强return olon, olat, rain

若需要插值可使用下面的代码，需要自行选择经纬度的范围和分辨率

olon = np.linspace (125,131,120) # 经纬坐标，0.05°分辨率
olat = np.linspace (44,47,60)         #  纬度坐标，0.05°分辨率
olon,olat = np.meshgrid(olon,olat)  # 生成坐标网格 meshgrid网格化
func = Rbf(lon,lat,rain,function='linear') #插值函数 调用Rbf插值函数中的 cubic 插值法linear
rain_data_new = func(olon,olat) #插值
rain_data_new[rain_data_new <0 ] = 0

2、等值线绘制和填充

def plot_contourf(olon, olat, rain, clevs, cdict):"""画等值线或色斑图:param olon: 经度:param olat: 纬度:param rain: 降水值:param clevs: 颜色级别:param cdict: 颜色列表:return:"""# 画布及绘图声明fig = plt.figure(figsize=(16, 9.6), facecolor='#666666', edgecolor='Blue', frameon=False)  # 画布ax = fig.add_subplot(111, projection=ccrs.PlateCarree())  # 绘图区my_cmap = colors.ListedColormap(cdict)  # 自定义颜色映射 color-mapnorm = mpl.colors.BoundaryNorm(clevs, my_cmap.N)  # 基于离散区间生成颜色映射索引#  绘制等值线、等值线填色cf = ax.contourf(olon, olat, rain, clevs, transform=ccrs.PlateCarree(), cmap=my_cmap, norm=norm)ct = ax.contour(olon, olat, rain, clevs)  # 绘制等值线ax.clabel(ct, fmt='%i')position = fig.add_axes([0.82, 0.2, 0.05, 0.2])  # 位置[左，下，宽。高]plt.colorbar(cf, cax=position)  # 颜色参照表position.set_yticklabels((0, 10, 25, 50, 100, 250, 500, 2000))ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))ax.yaxis.set_major_formatter(LatitudeFormatter())ax.set_xticks(np.arange(125,131, 2), crs=ccrs.PlateCarree())  # x轴ax.set_yticks(np.arange(44, 47, 2), crs=ccrs.PlateCarree())  # y轴ax.gridlines()  # 显示背景线plt.show()

我这里用的是这篇文章的数据，在此表示感谢，经过等值线的绘制和填充出来的图片如下，色级定义为：
clevs = [0.1, 10., 25., 50., 100., 250., 500] # 自定义色级，颜色列表
cdict = [‘#A9F090’, ‘#40B73F’, ‘#63B7FF’, ‘#0000FE’, ‘#FF00FC’, ‘#850042’, ‘#FF8C00’]

3、等值线和区域裁剪

上面的图片我们只留下感兴趣的区域，对于此外的区域直接裁剪掉。

 clip = maskout.shp2clip(cf, ax, './shijie.shp', 'hrb')

maskout函数的内容如下：

import shapefile
from matplotlib.path import Path
from matplotlib.patches import PathPatchdef getPathFromShp(shpfile, region):try:sf = shapefile.Reader(shpfile)vertices = [] codes = [] paths = []for shape_rec in sf.shapeRecords():# if shape_rec.record[3] == region: if shape_rec.record[2] == region:  pts = shape_rec.shape.pointsprt = list(shape_rec.shape.parts) + [len(pts)]for i in range(len(prt) - 1):for j in range(prt[i], prt[i + 1]):vertices.append((pts[j][0], pts[j][1]))codes += [Path.MOVETO]codes += [Path.LINETO] * (prt[i + 1] - prt[i] - 2)codes += [Path.CLOSEPOLY]path = Path(vertices, codes)paths.append(path)if paths:path = Path.make_compound_path(*paths)else:path = Nonereturn pathexcept Exception as err:print(err)return Nonedef shp2clip(originfig, ax, shpfile, region):path = getPathFromShp(shpfile=shpfile, region=region)patch = Noneif path:patch = PathPatch(path, transform=ax.transData, facecolor='none', edgecolor='black')for contour in originfig.collections:contour.set_clip_path(patch)return path, patch

其中code是用来记录等值点的标志，1表示开始，2表示中间，97表示结束。
裁剪后的图如下所示：

shp边界数据来源：
链接：https://pan.baidu.com/s/14t_TTgQLChy9_NJS1zCnmg
提取码：113d

END

参考文献：

1、 https://blog.csdn.net/qq_39425958/article/details/104790575
2、 https://gitee.com/cnxd365/qixiangshujukeshihua/blob/master/maskout.py