用到的主要有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，两个库都可以直接指定经纬度坐标，做省市的地图，但是很多参数没有看懂，比如需要指定坐标投影方法(不能省略)，我没搞懂如何利用自己手上已有的坐标(当地北东坐标)来进行绘图，所以依然坚持……不使用这俩。
利用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
图形边界文件涉及到地图信息隐私，所以我截取了一个角，能看出来已经满足了我的需求。
后续还要接着搞纵坐标格式、等值线图图形内坐标的标签。
转播
