pykrig克里金插值后绘制等值线图+边界外白化

参考Python遥感可视化 — Basemap将地面观测站点进行空间插值可视化和Python+Matplotlib画contour图这两篇文章
关于空间插值化，也可以直接参考pykrig的github文档

用到的主要有Matplotlib、Pykrig俩包,首先载入需要的数据文件

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from scipy.interpolate import griddata
import cmaps
 
# 用来正常显示中文
plt.rcParams['font.sans-serif'] = ['SimHei']
# 用来正常显示负号
plt.rcParams['axes.unicode_minus'] = False
# 获取污染物分布数据，其中包括维度、经度(或者以其他坐标系表示的x和y值)，和其他需要图形展现的污染物数据或者地形高度z值
df = pd.read_csv(r'E:\\xxxxx your pathway and you filename.dat',encoding='gbk')

# 剔除无效值NAN
df = df.dropna(axis=0, how='any')

# 获取纬度
lat = np.array(df["X"][:])
# 获取经度
lon = np.array(df["Y"][:])
# 获取温度数据，这里也可以是污染物数据或者地形高度z值
Temp = np.array(df["1-2m"][:])

# 创建格网1这种方法我还没看明白
#grid_x, grid_y = np.meshgrid(lat, lon)

#创建网格，我选择的是这种方法。linspace(min, max, 区分为多少间隔)间隔取值越大，计算所得的网格越密集
grid_x = np.linspace(502363, 503129,100)
grid_y = np.linspace(3067296, 3067809,60)

读取数据时可能会遇到一些编码问题，具体更改encoding可以解决，参考

#使用pykrig进行插值

from pykrige.ok import OrdinaryKriging
import pykrige.kriging_tools as kt

OK1 = OrdinaryKriging(lat, lon, Temp, variogram_model='linear')
# Create ordinary kriging object ignoring curvature:
OK2 = OrdinaryKriging(lat, lon, Temp, variogram_model="linear", verbose=False, enable_plotting=False)

#Execute on grid:
z1, ss1 = OK1.execute('grid', grid_x, grid_y)
z2, ss2 = OK2.execute('grid', grid_x, grid_y)

fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.imshow(z1)
ax1.set_title("geo-coordinates")
ax2.imshow(z2)
ax2.set_title("non geo-coordinates")
plt.show() #两种方法插值出来的效果很相似

也可以选用matplotlib来生成等值线图，效果更好看,这里参考Python+Matplotlib画contour图，后续我想再研究研究matplotlib中关于颜色的选择和坐标轴的设置，本土中画出来的坐标轴只有横轴是正确的，纵轴出现了问题。

# 生成等高线图
a = plt.contourf(grid_x, grid_y, z2, 15, cmap=plt.cm.Spectral)
b = plt.contour(grid_x, grid_y, z2,  15, colors='black', linewidths=1, linestyles='solid')
# 添加colorbar，ticks在这里可省略
plt.colorbar(a, ticks=[0, 0.25, 0.5, 0.75, 1])

#添加图内标签
plt.clabel(b, inline=True, fontsize=10)
plt.show()

Contour-example.png

后续：
希望能够将边界外的区域填充为白色。
看了很多文章……basemap功能非常强大，类似的有Cartopy，两个库都可以直接指定经纬度坐标，做省市的地图，但是很多参数没有看懂，比如需要指定坐标投影方法（不能省略），我没搞懂如何利用自己手上已有的坐标（当地北东坐标）来进行绘图，所以依然坚持……不使用这俩。

继续
我找到了 Plotting the geospatial data clipped by coastlines in Python这篇文章
，需求和我完全一致。使用已有地形文件的应该参考这篇文章可以完全解决问题。我没办法形成该文中masking的图形边界，但显然看到思路是在原有图形上加一个填充白色的图层。
我已有个bln边界文件，利用surfer另存为shp格式图形文件后，参考以下几个问题：#Geopandas Polygon to matplotlib patches Polygon conversion、
Plot shapefile with matplotlib、 Plot shapefile with matplotlib。

利用shapefile和descartes两个模块

import shapefile
#导入、形成图形文件
polys  = shapefile.Reader(r'E:\\your shape file.shp')
poly = polys.iterShapes().__next__().__geo_interface__

覆盖外界不需要的部分为白色

from descartes import PolygonPatch
# 生成等高线图
fig, ax = plt.subplots(figsize = (10, 5))
a = ax.contourf(grid_x, grid_y, z2, 15, cmap="inferno",alpha=0.6,)
b = ax.contour(grid_x, grid_y, z2,  15, colors='black', linewidths=1, linestyles='solid')
# 添加colorbar，ticks在这里可省略
plt.colorbar(a, ticks=[0, 0.25, 0.5, 0.75, 1])

#添加图内标签
#plt.clabel(b, inline=True, fontsize=10)

##屏蔽图形边界外侧区域
ax.add_patch(PolygonPatch(poly, fc='white', ec='white', alpha=1, zorder=2 ))
ax.axis('scaled')

Contour-example.png

图形边界文件涉及到地图信息隐私，所以我截取了一个角，能看出来已经满足了我的需求。
后续还要接着搞纵坐标格式、等值线图图形内坐标的标签。