程序由本公众号交流群友友提供，欢迎扫描文末二维码进群！

用途&效果

该程序对珠海市的臭氧多年季节平均值进行计算绘图，实现了 指北针添加，shpfile添加+mask白化，subplots多个子图制作排版，IDW插值算法等功能 ，效果如下，代码附后。

实现代码

"""
@author: licaoming
@contact: [email protected]
@software: jupyter notebook
@file: draw_zhuhai_ato_station_O3_season
@time: 2024/1/6
@function:提升20231226版本画图功能，画季节组图
@purpose:画站点数据图，插值方法IDW
"""

import pandas as pd
import numpy as np
import matplotlib.patches as mpatches
from math import radians, cos, sin, asin, sqrt 
import geopandas as gpd
from shapely.geometry import Point, Polygon
import matplotlib.pyplot as plt
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cartopy.feature as cfeat
import shapefile
from matplotlib.path import Path
#from cartopy.mpl.clip_path
from matplotlib.patches import PathPatch
import cartopy.crs as ccrs
import gma
from matplotlib.pyplot import savefig
#import cartopy as 
plt.rcParams['font.sans-serif']=['SimHei']

#读取站点数据
station_data=pd.read_excel("F:/课题/珠海/20231223/station_data_2.xlsx",sheet_name="臭氧季")#季节平均
#station_data=station_data.loc[station_data["要素"]=="臭氧"]#二氧化氮，臭氧
station_data

def haversine(lon1, lat1, lon2, lat2):
    # IDW插值方法
    R =  6372.8
    dLon = radians(lon2 - lon1)
    dLat = radians(lat2 - lat1)
    lat1 = radians(lat1)
    lat2 = radians(lat2)
    a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2
    c = 2*asin(sqrt(a))
    d = R * c
    return d
def IDW(x, y, z, xi, yi):
    lstxyzi = []
    for p in range(len(xi)):
        lstdist = []
        for s in range(len(x)):
            d = (haversine(x[s], y[s], xi[p], yi[p]))
            lstdist.append(d)
        sumsup = list((1 / np.power(lstdist, 2)))
        suminf = np.sum(sumsup)
        sumsup = np.sum(np.array(sumsup) * np.array(z))
        u = sumsup / suminf
        xyzi = [xi[p], yi[p], u]
        lstxyzi.append(xyzi)
    return(lstxyzi)
    
def add_north(ax, labelsize=18, loc_x=0.1, loc_y=0.93, width=0.05, height=0.08, pad=0.1):
    """
    画一个比例尺带'N'文字注释
    主要参数如下
    :param ax: 要画的坐标区域 Axes实例 plt.gca()获取即可
    :param labelsize: 显示'N'文字的大小
    :param loc_x: 以文字下部为中心的占整个ax横向比例
    :param loc_y: 以文字下部为中心的占整个ax纵向比例
    :param width: 指南针占ax比例宽度
    :param height: 指南针占ax比例高度
    :param pad: 文字符号占ax比例间隙
    :return: None
    """
    minx, maxx = ax.get_xlim()
    miny, maxy = ax.get_ylim()
    ylen = maxy - miny
    xlen = maxx - minx
    left = [minx + xlen*(loc_x - width*.5), miny + ylen*(loc_y - pad)]
    right = [minx + xlen*(loc_x + width*.5), miny + ylen*(loc_y - pad)]
    top = [minx + xlen*loc_x, miny + ylen*(loc_y - pad + height)]
    center = [minx + xlen*loc_x, left[1] + (top[1] - left[1])*.4]
    triangle = mpatches.Polygon([left, top, right, center], color='k')
    ax.text(s='N',
            x=minx + xlen*loc_x,
            y=miny + ylen*(loc_y - pad + height),
            fontsize=labelsize,
            horizontalalignment='center',
            verticalalignment='bottom')
    ax.add_patch(triangle)
    
def maskout_areas(originfig, ax, shp_path, proj=None):
    # mask白化
    sf = shapefile.Reader(shp_path,encoding='gbk')
    vertices = []
    codes = []
    for shape_rec in sf.shapeRecords():
        if shape_rec.record[2] == '斗门区' or shape_rec.record[2] == '金湾区' or shape_rec.record[2] == '香洲区':
            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]):
                    if proj:
                        vertices.append(proj.transform_point(pts[j][0], pts[j][1], ccrs.Geodetic()))
                    else:
                        vertices.append((pts[j][0], pts[j][1]))
                codes += [Path.MOVETO]
                codes += [Path.LINETO] * (prt[i + 1] - prt[i] - 2)
                codes += [Path.CLOSEPOLY]
            clip = Path(vertices, codes)
            clip = PathPatch(clip, transform=ax.transData)
    for contour in originfig.collections:
        contour.set_clip_path(clip)
        
#创建经纬格点
entent = [113, 114.4,21.8, 22.6]
grid = 20
Lon1 = entent[0]
Lon2 = entent[1]
Lat1 = entent[2]
Lat2 = entent[3]
grid_lon_list=np.linspace(Lon1,Lon2,grid)
grid_lat_list=np.linspace(Lat1,Lat2,grid)
# 转换成网格
newlon, newlat = np.meshgrid(grid_lon_list, grid_lat_list)
#数据扁平化
grid_lon_list = newlon.flatten().tolist()
grid_lat_list = newlat.flatten().tolist()

#组图
#画季节平均值图
datayear =["春季","夏季","秋季","冬季"]
latitude =[]
longitude=[]
idw_data=[]
for dy in datayear:
    #对站点数据进行IDW插值
    lon=station_data["经度"].tolist()
    lat=station_data["纬度"].tolist()
    data = station_data[dy]
    sta_data=np.array(IDW(lon,lat,data,grid_lon_list,grid_lat_list))
    sta_lon=sta_data[:,0].reshape(20,20)[0:17,0:10]
    sta_lat=sta_data[:,1].reshape(20,20)[0:17,0:10]
    sta_idw=sta_data[:,2].reshape(20,20)[0:17,0:10]
    latitude.append(sta_lat)
    longitude.append(sta_lon)
    idw_data.append({dy:sta_idw})
season_data=[i for i in idw_data[1].values()][0]

#画图
num=0
entent = [113,113.7, 21.8, 22.5]
shp_path = r"F:/shp/zhuhaimap/zhuhaimap/zhuhai_new.shp"
#建立画布
crs=ccrs.PlateCarree() 
fig,axes =plt.subplots(nrows=2,ncols=2,figsize=(30,30),sharex=True,sharey=True,subplot_kw={'projection': crs},constrained_layout=True,dpi=100)
for x, row in enumerate(axes):
    for y,col in enumerate(row):
        #数据
        season_lon=longitude[num]
        season_lat=latitude[num]
        season_data=[i for i in idw_data[num].values()][0]
        season=[i for i in idw_data[num].keys()][0]
        #设置图片范围
        axes[x,y].set_extent(entent, crs=crs)
        #绘制数据图
        clevs = np.arange(66,128,2)
        cf = axes[x,y].contourf(season_lon,season_lat, season_data,clevs,cmap="turbo",tranform=crs)#"nipy_spectral"
        #裁剪
        maskout_areas(originfig=cf, ax=axes[x,y], shp_path=shp_path,proj=None)
        #添加矢量地图
        shp = gpd.read_file(shp_path)
        zh = cfeat.ShapelyFeature(shp.geometry,ccrs.PlateCarree(), edgecolor='k',facecolor='none'