嘿嘿，让我们采取倒叙的方式，先看成品，再看制作过程。
jz.png
- 准备焦作市各区县的面图文件，具体操作参考宝藏文章，一步一步地就做好了。
Arcgis(一) 制作全国行政区shp文件(精确到县级)_皖渝的博客-CSDN博客_arcgis制作shp文件
需要用到全国地理信息库：
成果数据1:100万 (webmap.cn)
需要电脑上有arcgis软件。
安装一些python库
参考python geopandas库安装 - 知乎 (zhihu.com)我的代码
## 导入库，定义函数
import pandas as pd
import numpy as np
from pykrige.ok import OrdinaryKriging
import plotnine
from plotnine import *
import geopandas as gpd
import shapefile
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import cmaps
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from scipy.interpolate import griddata
#from cartopy.io.shapereader import Reader
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
def shp2clip(originfig, ax, shpfile):
sf = shapefile.Reader(shpfile)
vertices = []
codes = []
for shape_rec in sf.shapeRecords():
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]
clip = Path(vertices, codes)
clip = PathPatch(clip, transform=ax.transData)
for contour in originfig.collections:
contour.set_clip_path(clip)
return contour
## 读数据
df = pd.read_excel(r"D:\jz_21sta_MonthMean.xlsx",sheet_name="jz_18")
pos1 = pd.read_excel(r"D:\所有空气点位经纬度登记表.xlsx",sheet_name="sta-18")
pos2 = pd.read_excel(r"D:\所有空气点位经纬度登记表.xlsx",sheet_name="sta-12")
## 循环绘图
j= -1
need_time = ["2019/1/31","2019/2/28","2019/3/31","2019/4/30","2019/5/31","2019/6/30","2019/7/31",
"2019/8/31","2019/9/30","2019/10/31","2019/11/30","2019/12/31","2020/1/31",
"2020/2/29","2020/3/31","2020/4/30","2020/5/31","2020/6/30","2020/7/31",
"2020/8/31","2020/9/30","2020/10/31","2020/11/30","2020/12/31",
"2021/1/31","2021/2/28","2021/3/31","2021/4/30","2021/5/31","2021/6/30",
"2021/7/31","2021/8/31","2021/9/30","2021/10/31","2021/11/30","2021/12/31"]
figname = ["2019-1","2019-2","2019-3","2019-4","2019-5","2019-6",
"2019-7","2019-8","2019-9","2019-10","2019-11","2019-12",
"2020-1","2020-2","2020-3","2020-4","2020-5","2020-6","2020-7",
"2020-8","2020-9","2020-10","2020-11","2020-12","2021-1","2021-2","2021-3","2021-4","2021-5",
"2021-6","2021-7","2021-8","2021-9","2021-10","2021-11","2021-12"]
for i in need_time[:28]:
j = int(j)+1
name = str(j)+"month_mean_"
##2020年5月之后只有12个站的数据
if j>16:
pos = pos2
lons = pos2["经度"]
lats = pos2["纬度"]
else:
pos = pos1
lons = pos1["经度"]
lats = pos1["纬度"]
data = pd.DataFrame(df[df["date"]==i]["O3"])
min_data = data.min().values
max_data = data.max().values
#print(data)
grid_lon = np.linspace(112.4, 113.8,140)
grid_lat = np.linspace(34.7, 35.6,90)
OK = OrdinaryKriging(lons, lats, data, variogram_model='linear',nlags=6)
z1, ss1 = OK.execute('grid', grid_lon, grid_lat)
#z1.shape
xgrid, ygrid = np.meshgrid(grid_lon, grid_lat)
df_grid = pd.DataFrame(dict(long=xgrid.flatten(),lat=ygrid.flatten()))
df_grid["Krig_gaussian"] = z1.flatten()
#df_grid.to_csv(name+".csv",encoding='utf_8_sig')
sh = gpd.read_file(r'D:\arcgis\jiaozuo\final\jiaozuo_result.shp')
#sh
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_extent([112.5, 113.7, 34.75, 35.55], crs=ccrs.PlateCarree())
province = shpreader.Reader(r"D:\arcgis\jiaozuo\final\jiaozuo_result.shp")
ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=1,edgecolor='k',facecolor='none')
#cdict=['white', '#A6F28F','#3DBA3D','#61B8FF','#0000FF','#FA00FA','#800040','#60281E']#自定义颜色列表
cf = ax.contourf(xgrid, ygrid, z1, levels=np.linspace(120,130,21),cmap="Spectral_r",
transform=ccrs.PlateCarree())
cb = plt.colorbar(cf,shrink=0.8,pad=0.01)
cb.set_label('O${_3}$($\mu$g/m${^3}$)',fontsize=15)
shp2clip(cf, ax, r'D:\arcgis\jiaozuo\final\jiaozuo_result.shp')
#ax.set_extent([112., 114., 34, 36])
ax.scatter(lons,lats,c=data["O3"],s=90,ec="k",lw=0.5,cmap="Spectral_r",vmin = 120,vmax =130)
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.set_xticks(np.arange(112.6,113.7,0.2))
ax.set_yticks(np.arange(34.8,35.5,0.1))
ax.tick_params(labelsize=18) #坐标轴刻度字体大小18
ax.set_title(figname[j],fontsize=22)
plt.show()
fig.align_ylabels(ax)
fig.savefig(name+'.png',bbox_inches = 'tight',dpi=300)
为了严谨，coordinates_type="geographic"应该补充在OrdinaryKriging这个函数里。不过对于中纬度而言，加与不加差距很小了。克里金插值的coordinates_type参数 - Heywhale.com
