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-map
    norm = 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 PathPatch


def 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.points
               prt = 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 = None
       return path
   except Exception as err:
       print(err)
       return None
       
def shp2clip(originfig, ax, shpfile, region):
   path = getPathFromShp(shpfile=shpfile, region=region)
   patch = None
   if 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