【Python教程】气象数据空间插值与绘图

2184

Python气象数据空间插值与绘图

本教程：简述利用Python进行空间插值、气象等值线绘制的基本工具和方法。泛泛之谈难免有遗漏，诚邀各位指正、交流讨论。

在气象领域，所能获取的气温、降水、风速等数据常是以离散点形式存在的气象站点值，特别是一些位置偏僻、难以到达的地区存在气象数据缺测，即在时空上的数据空白。这对流域、区域乃至全球的连续场分析造成一定阻碍。

通过本贴，可实现的气象数据处理结果如下：

图1 基于Python的气象数据空间插值与等值线效果图

0. 预备知识：空间插值——由点推面的过程

空间插值，提供了一种途径，利用有限个采样点（位置、属性值），并通过空间自相关性等原理，推测空间随机场的分布与恢复，简言之就是以点推面、复原空间空间连续分布（空间样本值+坐标 →插值 →连续场）。这需要找出一个合适的函数，其中自变量是已知点的值及其空间关系，因变量是最终需要估算的位置的值。

图2 空间插值

*常见插值方法*：

（1）基于权重：泰森多边形、反距离权重法、自然邻域

泰森多边形（Thiessen Polygons或），又称Voronoi图，是非常经典的用于计算区域降雨量的方法，由美国气候学家Alfred H. Thiessen提出。将区域划分为数个多边形，用位于多边形中或附近的样本点赋值，实际上是用的最邻近方法而不是插值，最邻近样本点代替区域内未知点、忽略未知点信息。因此，这导致预测的精确性有问题，多边形边界内无变化、边界处突变，插值结果构造的不是一个曲面、不是空间要素的真实分布。

图3 泰森多边形

反距离权重法（IDW）对周围样点按距离赋予不同的权重，近的权重大，远小。

图4 反距离权重

（2）克里金插值

克里金插值由南非地质学家、工程师克里金（Krige）在20世纪50年代提出。该方法最初用于采矿勘探，后推广到气象等多领域。细分为简单克里金、普通克里金、泛克里金、回归克里金等多个小类。详见代码部分。

图5 克里金插值

（3）回归模型和多项式

例如，全局回归模型：最小二乘法等。利用一个变量和其他一个（一元回归）/多个（多元回归）变量的相关性（线性/非线性互相关性），建立回归关系进行计算。整个研究区只有一个方程，系数适用于整个研究区，但不一定可推广到其他区域。地理加权回归：局部空间回归方法，基于空间异质性（地理学第二定律），把研究区划分为小区域，每个小区域用一个回归方程拟合。

趋势面：利用一个数学函数（多项式）定义的平滑表面与所有样本点进行建模。径向基函数（RBF）：非线性核函数，常用于多项式不可分或复杂非线性问题。样条函数。

图6 多项式回归

（4）其他新方法：机器学习和专业模型

*软件与工具准备*：

（1）Python-Anaconda

请参考论坛教程版已有的介绍帖子。

（2）工具包安装

Pandas、numpy、matplotlib请参考“论坛教程”版已有的帖子。

cartopy：由英国气象局开发的Python 包，专用于地理空间数据处理，即生成地图和其他地理空间数据分析。安装代码任选其一，在电脑终端（Terminal）命令行输入指令，如果是在Jupyter notebook的.ipynb文件则指令前加“!”，以下其他工具包同理。（https://scitools.org.uk/cartopy/docs/latest/reference/index.html）

1. pip install cartopy
   
2. pip install git+https://github.com/SciTools/cartopy.git

复制代码

pykrige：克里金插值。安装代码任选其一：如果安装了anaconda则推荐前两个代码，其他使用pip。（https://geostat-framework.readth ... atest/contents.html）

1. conda install pykrige
   
2. conda install -c conda-forge pykrige
   
3. pip install pykrige

复制代码

metpy：多用途气象包，读写、计算、可视化气象数据，包含插值工具。https://unidata.github.io/MetPy/latest/index.html

1. pip install metpy
   
2. conda install -c conda-forge metpy

复制代码

scipy：用于算法和数据结构领域，包含插值工具。https://docs.scipy.org/doc/scipy/index.html

1. python -m pip install scipy
   
2. conda install scipy
   
3. sudo apt-get install python3-scipy # 用于Ubuntu和Debian
   
4. brew install scipy # 用于macOs

复制代码

1. 空间插值实例与代码（见下文）

2. 空间插值结果（见下文）

2184

1. 空间插值实例与代码

数据准备：

案例数据是根据土壤剖面推导的黄土高原气温、降水（下载地址：https://doi.org/10.6084/m9.figshare.25263637.v1）。

1. import pandas as pd
   
2. import numpy as np
   
3. df = pd.read_excel("/下载/25263637/Paleoclimate data of the 180 modern soils.xlsx",sheet_name='Data')

复制代码

读取结果：

Sample        Longitude (°)        Latitude (°)        Precipitation (mm)        Temperature (°C)        cfd (m3/kg)        Hematite (g/kg)
0        MS-1        108.355767        37.709567        356.361        8.333        1.918159        1.912016
1        MS-2        107.626517        37.418967        331.914        8.286        1.329787        2.068267
2        MS-3        104.430650        35.711850        365.322        6.643        1.756757        2.861642
3        MS-4        104.765900        35.665383        380.359        5.960        4.011461        2.591716
4        MS-5        105.158983        35.570183        397.854        5.653        2.631579        2.836020
...        ...        ...        ...        ...        ...        ...        ...
175        MS-176        105.260210        35.441560        401.010        5.666        4.048583        2.809483
176        MS-177        105.019330        35.715930        387.716        5.856        2.200825        2.888728
177        MS-178        104.951970        35.944250        356.024        6.589        5.738881        3.541073
178        MS-179        103.620770        35.724660        368.705        7.873        0.257069        3.444593
179        MS-180        106.704390        37.571310        295.204        8.147        0.721154        1.745166
180 rows × 7 columns

1.1 Kriging

路径：数据读入→数据预处理→变差函数计算→变差函数拟合→插值→等值线绘制

（0）相关概念

（I）平稳假设

区域化变量的随机函数F(x)在各阶矩相对恒定，一般只要求一阶矩（数学期望）、二阶矩（方差）平稳，即在随机场中不同点的随机变量产生的概率分布相等。

（II）无偏估计。克里金权重之和等于1。

（III）变差函数（variogram或semi-variogram）

含义：空间随机场上两个点形成的随机变量，随机变量方差之大作为函数。两点（i, j）距离h越小、方差越小。数学含义为：方差与两点的协方差（cov(i,j)）相减，等价于两点（x=i, j）在F(x)上两个函数值之差的方差的二分之一。

类型（代码中的variogram_model）：1.球状模型（spherical）。2.指数模型（exponential）。3.高斯模型（gaussian）。

重要参数：①块金（nugget）：随机场粗糙程度。②基台值（sill）：空间随机场起伏程度。③变程（range）：空间相关性存在范围。pykrige还有nlags：对于半方差模型的平均箱数（bins），默认为6。

（1）普通克里金（Ordinary Kriging）

方程组考虑条件极值问题，采用拉格朗日算法。数学期望未知，只要满足二阶平稳假设。

利用pykrige工具的函数OrdinaryKriging()插值黄土高原气温点。可以选择半变异函数（variogram_model）类型。调整对于半方差模型的平均箱数（nlags，默认为6）：是将成对的地点归为一个距离类别的大小。

1. # Ordinary Kriging
   
2. from pykrige.ok import OrdinaryKriging
   
3. # 建立插值网格，分辨率0.25°×0.25°（即~25km）。oklon为经度，oklat为纬度
   
4. oklon = np.arange(df['Longitude (°)'].min(),df['Longitude (°)'].max(),0.25)
   
5. oklat = np.arange(df['Latitude (°)'].min(),df['Latitude (°)'].max(),0.25)
   
6. # 普通克里金插值：当 半变异函数（variogram_model）选择 spherical
   
7. OK = OrdinaryKriging(df['Longitude (°)'],df['Latitude (°)'],df['Temperature (°C)'],
   
8.                      variogram_model="spherical",verbose=False,enable_plotting=False,coordinates_type='geographic',nlags=6)
   
9. temp1, ss1 = OK.execute('grid',oklon,oklat)

复制代码

插值得到的气温面数据temp1和原始散点数据df绘图：

1. # 绘图
   
2. import matplotlib.pyplot as plt
   
3. import cartopy.crs as ccrs # 坐标系
   
4. import cartopy.feature as cfeat
   
5. 
   
6. extents = [103, 120, 34, 41] # 经纬度范围
   
7. crs = ccrs.PlateCarree() # 坐标系
   
8. fig = plt.figure() # 可调整，例如figsize=(10,8),dpi=500
   
9. ax1 = fig.add_subplot(111, projection=crs) # 创建子图，111是把画板分成1行1列的第1张图
   
10. x_extent = [105, 110, 115, 120]
    
11. y_extent = [30, 35, 40, 45]
    
12. 
    
13. ax1.set_title('Temperature of Loess Plateau (Ordinary Kriging)',loc='center')
    
14. ax1.set_xticks(x_extent, crs=ccrs.PlateCarree()) #添加边框的经纬度标记
    
15. ax1.set_yticks(y_extent, crs=ccrs.PlateCarree())
    
16. ax1.set_extent(extents, crs)
    
17. # 绘制海洋、河流、海岸线，需连接网络！！！！！
    
18. ax1.add_feature(cfeat.OCEAN.with_scale('110m'), linewidth=0.5, zorder=1)
    
19. ax1.coastlines()
    
20. # 绘制插值结果（imshow）和站点实际值（scatter，散点颜色、边框、坐标系、与其他图的位置）
    
21. s1 = ax1.imshow(temp1,origin='lower', cmap='RdYlGn_r', transform=crs,extent=[df['Longitude (°)'].min(), df['Longitude (°)'].max(), df['Latitude (°)'].min(), df['Latitude (°)'].max()])
    
22. ax1.scatter(df['Longitude (°)'], df['Latitude (°)'], c=df['Temperature (°C)'], cmap='RdYlGn_r', edgecolors='black', transform=ccrs.PlateCarree(), zorder=2)
    
23. cb1 = fig.colorbar(s1,ax=ax1,orientation='horizontal') # 水平色带图例
    
24. cb1.set_label('°C')
    
25. plt.show()

复制代码

绘图结果：

（2）泛克里金（universal kriging）

用于解决非平稳随机场的方法，在数学期望和协方差不同（协方差不平稳）时使用，可看成在局部范围内平稳。比较不同nlags：

1. from pykrige.uk import UniversalKriging
   
2. import numpy as np
   
3. # 建立插值网格，分辨率0.25°×0.25°（即~25km）
   
4. uklon = np.arange(df['Longitude (°)'].min(),df['Longitude (°)'].max(),0.25)
   
5. uklat = np.arange(df['Latitude (°)'].min(),df['Latitude (°)'].max(),0.25)
   
6. # 克里金插值：当半变异函数（variogram_model）选择 spherical时
   
7. uk = UniversalKriging(df['Longitude (°)'],df['Latitude (°)'],df['Temperature (°C)'],
   
8.                      variogram_model='spherical',verbose=False,enable_plotting=False)# 默认nlags=6
   
9. tempuk, ssuk = uk.execute('grid',uklon,uklat)
   
10. uk1 = UniversalKriging(df['Longitude (°)'],df['Latitude (°)'],df['Temperature (°C)'],
    
11.                      variogram_model='spherical',verbose=False,enable_plotting=False,nlags=1)
    
12. tempuk1, ssuk1 = uk1.execute('grid',uklon,uklat)

复制代码

绘图代码同普通克里金，结果如图：

（3）协同克里金（Co-Kriging）

使用条件是，主变量的信息不足，而次要变量的信息充分（例如主变量为温度，气象台站点稀少，用更多辅助变量（高程）协助插值）；主次变量存在一定相关性（例如温度和高程）；估值时在搜索半径内至少存在一个主变量采样点。

（4）回归克里金（Regression kriging）

1. # RK
   
2. # 具体请参考https://geostat-framework.readthedocs.io/projects/pykrige/en/latest/generated/pykrige.rk.RegressionKriging.html
   
3. from pykrige.rk import RegressionKriging
   
4. RegressionKriging()

复制代码

2184

1.2 IDW

反距离权重法（IDW）对周围样点按距离（d）赋予不同的权重，近的权重大，远小。权重的衰减采用距离的倒数（1/d或1/d^2或1/d^3...）计算：样本到实测点反比占整个距离反比和的比重作为权重，和距离反比成正比，近-反比大-权重大。

方法一：利用metpy的inverse_distance_to_grid()函数。

1. from metpy.interpolate import inverse_distance_to_grid
   
2. idwlon = np.arange(df['Longitude (°)'].min(),df['Longitude (°)'].max(),0.25)
   
3. idwlat = np.arange(df['Latitude (°)'].min(),df['Latitude (°)'].max(),0.25)
   
4. idwlon_mes,idwlat_mes = np.meshgrid(idwlon,idwlat)
   
5. # 可调整的参数：r：从网格中心算起的半径，在此半径内的观测数据将被考虑并加权。min_neighbors：对一个点进行插值所需的最小邻点数，默认值为 3
   
6. IDW = inverse_distance_to_grid(df['Longitude (°)'],df['Latitude (°)'],df['Temperature (°C)'],
   
7.                                idwlon_mes,idwlat_mes, r=1, min_neighbors=1)

复制代码

绘图代码同普通克里金，结果图如下：

方法二：

1. # 定义距离函数
   
2. def distance_matrix(x0, y0, x1, y1):
   
3.     obs = np.vstack((x0, y0)).T
   
4.     interp = np.vstack((x1, y1)).T
   
5.     # 两个观测值之间距离的矩阵
   
6.     d0 = np.subtract.outer(obs[:,0], interp[:,0])
   
7.     d1 = np.subtract.outer(obs[:,1], interp[:,1])
   
8.     return np.hypot(d0, d1)
   
9. # 定义idw插值函数
   
10. def simple_idw(x, y, z, xi, yi):
    
11.     distance = distance_matrix(x,y, xi,yi)
    
12.     weights = 1.0 / distance # 权重
    
13.     weights /= weights.sum(axis=0) # 权重和为1
    
14.     z_idw = np.dot(weights.T, z)
    
15.     return z_idw
    
16. # 数据准备
    
17. lonnumber = int((df['Longitude (°)'].max()-df['Longitude (°)'].min())*10) # 列数
    
18. latnumber = int((df['Latitude (°)'].max()-df['Latitude (°)'].min())*10)
    
19. idwlon = np.linspace(df['Longitude (°)'].min(), df['Longitude (°)'].max(), lonnumber)
    
20. idwlat = np.linspace(df['Latitude (°)'].min(), df['Latitude (°)'].max(), latnumber)
    
21. idwlon, idwlat = np.meshgrid(idwlon, idwlat)
    
22. idwlon, idwlat = idwlon.flatten(), idwlat.flatten()
    
23. # 调用插值函数并绘图
    
24. grid1 = simple_idw(df['Longitude (°)'], df['Latitude (°)'], df['Temperature (°C)'],idwlon,idwlat)
    
25. grid1rsh = grid1.reshape((lonnumber, latnumber))
    
26. plt.scatter(idwlon, idwlat, c=grid1rsh, cmap='RdYlGn_r')
    
27. plt.title('Temperature IDW')

复制代码

1.3 Metpy的其他算法

Metpy还提供多种插值工具，例如自然邻域（natural neighbor）、interpolate_to_grid函数可选算法。

自然邻域插值（natural_neighbor_to_grid()）：根据每个样本点假定的“影响区域”为每个样本点创建权重。这些区域由每个样本点周围生成的Voronoi多边形确定。原则上，创建的每个网格交点都将位于其中一个多边形中，并且可以分配创建多边形的点的值。这将导致斑块的阶梯状表面。（(2010). A stable and fast implementation of natural neighbor interpolation.）

1. from metpy.interpolate import natural_neighbor_to_grid
   
2. natural_neighbor_to_grid()

复制代码

interpolate_to_grid函数：interp_type参数提供Metpy和Scipy的多个插值算法：“linear”、“nearest”、“cubic”、“natural_neighbor”、 “barnes”、“cressman”。

1. from metpy.interpolate import interpolate_to_grid
   
2. # 参数自行调整。search_radius 插值的搜索半径，则默认观测值平均间距的5倍。hres 生成网格的水平分辨率，单位与xp和yp参数相同，默认50000。
   
3. itg_ln_lon,itg_ln_lat,itg_ln = interpolate_to_grid(xp,yp,df1['Temperature (°C)'].values,interp_type='natural_neighbor', minimum_neighbors=1,search_radius=400, hres=40000)

复制代码

2184

1.4 Griddata

利用工具scipy，参数method可选linear、nearest、cubic，试做比较：

1. # griddata
   
2. from scipy.interpolate import griddata
   
3. gridlon = np.arange(df['Longitude (°)'].min(),df['Longitude (°)'].max(),0.25)
   
4. gridlat = np.arange(df['Latitude (°)'].min(),df['Latitude (°)'].max(),0.25)
   
5. gridlon_mes,gridlat_mes = np.meshgrid(gridlon,gridlat)
   
6. GD = griddata((df['Longitude (°)'], df['Latitude (°)']),df['Temperature (°C)'],(gridlon_mes,gridlat_mes),method='cubic')

复制代码

三个method绘图结果，绘图代码同普通克里金：

1.5 Scipy的其他插值算法

径向基函数Rbf：

1. from scipy.interpolate import Rbf
   
2. rbflon = np.arange(lon.min(),lon.max(),0.25)
   
3. rbflat = np.arange(lat.min(),lat.max(),0.25)
   
4. rbflon,rbflat = np.meshgrid(rbflon,rbflat)
   
5. rbf_ = Rbf(df['Longitude (°)'], df['Latitude (°)'],df['Temperature (°C)'],function='multiquadric')
   
6. rbf_values = rbf_(rbflon,rbflat)

复制代码

泰森多边形Voronoi：

1. from scipy.spatial import Voronoi,voronoi_plot_2d
   
2. Voronoi(point)

复制代码

2184

2.空间插值结果比较与绘制等值线图

2.1 几种插值结果的比较

（1）不同模型比较：现有条件下，克里金方法的结果相对较优（如下图）。①克里金精度较好，较好地满足估计值的数学期望和真实值数学期望相同的“无偏”条件，且估计值和真实值之间的方差最小。②相比于idw和griddata，kriging能更好地实现对采样点分布稀疏区域或研究区边缘的插值外推。③当采样点局部聚集，采样的有偏导致估计有偏时，克里金可以比idw更好地处理clustery现象，不使估计值倾向于样本点聚集处。④根据方向选择协方差、变差函数函数作为系数，兼顾各向异性，还可显示估计方差。

（2）模型内部算法比较：

①griddata的三个插值算法：nearest即最邻近法，将距离插值目标格点最近的采样点值作为格点值，这使得数据真实性较好、并且可外推不规则插值区的边缘格点值，但插值结果呈现斑块状、空间连续性差。cubic插值面连续性最好、但数值有可能失真，且无法外推格点。linear介于前两方法之间，兼顾数据真实性与空间场的平滑程度，但也无法外推格点。

nearest：，

linear：，

cubic：

②克里金三算法比较：球状模型（spherical）、高斯模型（gaussian）和指数模型（expotential）。球状模型和指数模型插值结果较优，高斯模型结果在区域西北边缘存在异常高值区，这可能说明高斯模型不如前二者适合插值采样点稀缺的边缘区。

球状：，

指数：，

高斯：

（3）模型内部参数比较：内部参数，影响搜索值的范围，使插值结果存在聚集或分散的差异。

①idw：

r=min_neighbors=1，
r=min_neighbors=3

②克里金：nlags=1或6

可见，不同数据特征、不同参数组合将影响插值精度。数据离散或聚集特征、空白区、分布不均匀、空值等特征引起插值结果不同，甚至在某些算法中出错或无法计算。因此需认真了解数据、选择适宜的插值方法与参数。

（未完待续……）

yrch007

个人应用之谈，在Python中遇到无法外推的插值方法，其实也可以通过设立虚拟数据点的方式来实现，只要遵循地理第一定律（越近的事物越相关），可适当延展插值区域到所需的绘图边界外。

2184

yrch007 发表于 2024-6-18 19:31
个人应用之谈，在Python中遇到无法外推的插值方法，其实也可以通过设立虚拟数据点的方式来实现，只要遵循地 ...

赞同。比如，制作泰森多边形就会用到虚拟点

2184

2.2 插值图完善：添加等值线

利用matplotlib的contourf、contour函数绘制气温、降水等值线。需设置等值线间隔参数levels=np.arange(3,15.5,1)，即3~15.5℃区间每1℃一个等值线/间隔。将普通克里金代码的imshow等行替换为以下代码即可：

1. # 绘制降水插值结果等值线（contourf、contour）
   
2. cf1 = ax1.contourf(oklon,oklat,temp1, levels=np.arange(3,15.5,1),origin='lower',cmap='RdYlGn_r',transform=crs,extent=[df['Longitude (°)'].min(),df['Longitude (°)'].max(),df['Latitude (°)'].min(),df['Latitude (°)'].max()])
   
3. ct1 =ax1.contour(oklon,oklat,temp1, colors='k', levels=np.arange(3,15.5,1), linewidths=0.5,transform=crs)
   
4. lb1 = plt.clabel(ct1, fontsize=6, inline=True, fmt='%r') # 等值线注记

复制代码

结果图：

2184

1.6 泰森多边形

泰森多边形可以看作找到最接近的“任何东西”，即在根据采样点A生成的多边形内，距离最近的采样点即为点A。泰森多边形的构建是将采样点连成三角网（即Delaunay三角网（TIN）），步骤为：

①连接采样点形成三角网，并对点和三角形编号，如果某点与前一个点的距离小于给定邻近值，则分析时忽略该点；②作出每个三角形边的中垂线，由这些中垂线构成泰森多边形的边，而中垂线的交点是泰森多边形的顶点。

需注意：①任何一个三角形的外接圆内不能包含任何其它离散点；②相邻两个三角形构成凸四边形，在交换对角线之后，六个内角的最小者不再增大。

首先，让我们来进行数据预处理。在此，建立每个黄土高原降水点的经纬度点对coords_df，利用geopandas将转换为具有空间属性（依托ccrs.PlateCarree()坐标系的经纬度）的矢量数据df_pts。

为了解决前文“研究区边缘无法外推”的问题，需要建立虚拟点，即在研究区外围增加点位：建立每个降水点的缓冲区buffer，选取缓冲区边缘最偏东/西/南/北四方的边缘点，与降水点合并为待插值点集coords_voro。

1. import numpy as np
   
2. from scipy.spatial import Voronoi,voronoi_plot_2d
   
3. import geopandas as gpd
   
4. from shapely.geometry import Polygon, Point
   
5. import matplotlib.pyplot as plt
   
6. import cartopy.crs as ccrs
   
7. import cartopy.feature as cfeat
   
8. proj =ccrs.PlateCarree()
   
9. coords_df = [list(xy) for xy in zip(df['Longitude (°)'],df['Latitude (°)'])]
   
10. df_pts = gpd.GeoDataFrame(df['Precipitation (mm)'],geometry = [Point(x, y) for x, y in coords_df], crs = proj)
    
11. min_lon_voro, min_lat_voro, max_lon_voro, max_lat_voro = df_pts.buffer(2).total_bounds
    
12. coords_voro = coords_df + [[min_lon_voro,min_lat_voro],[max_lon_voro,min_lat_voro],[max_lon_voro,max_lat_voro],[min_lon_voro,max_lat_voro]]

复制代码

缓冲区图

2184

用Voronoi(coords_voro)构建泰森多边形voro。基于顶点voro.vertices创建多边形，gpd.GeoDataFrame()创建空间矢量图层voro_polys，最后把多边形内的降雨点值赋值给多边形。

1. voro = Voronoi(coords_voro) # vpts,coords_df,coords_voro
   
2. voro_poly_list=[]
   
3. for region in voro.regions:
   
4.     if -1 in region:
   
5.         continue
   
6.     else:
   
7.         pass
   
8.     # 选择非0值区域
   
9.     if len(region) != 0:
   
10.         voro_poly_region = Polygon(list(voro.vertices[region]))
    
11.         voro_poly_list.append(voro_poly_region)
    
12.     else:
    
13.         continue
    
14. # 构建空间矢量数据
    
15. voro_polys = gpd.GeoDataFrame(voro_poly_list, columns = ['geometry'], crs = proj)
    
16. # 将泰森多边形内降雨点的数值赋给多边形
    
17. voro_polys_values = gpd.sjoin(df_pts , voro_polys, how = "right", op = 'within')
    
18. display(voro_polys_values.head())

复制代码

多边形属性表：

index_left        Precipitation (mm)        geometry
0        121        595.031        POLYGON ((112.902 31.420, 111.763 29.659, 112....
1        23        412.586        POLYGON ((112.722 39.365, 112.768 39.510, 114....
2        2        365.322        POLYGON ((104.493 34.931, 104.635 35.954, 104....
3        178        368.705        POLYGON ((103.999 34.013, 103.804 33.284, 96.8...
4        117        603.988        POLYGON ((113.148 34.681, 112.972 34.601, 112....

最后绘制降水点及其泰森多边形插值效果图，省界图层为附件：

1. # 设定图幅
   
2. fig, ax = plt.subplots(1, 1, figsize = (8, 5))
   
3. # 图形属性与显示范围
   
4. plt.style.use('bmh')
   
5. plt.xlim(103,115)
   
6. plt.ylim(33,41)
   
7. plt.xlabel('longitude')
   
8. plt.ylabel('latitude')
   
9. # 依次叠加泰森多边形、降雨点、省界三图层
   
10. voro_polys_values.plot(column='Precipitation (mm)',ax = ax, cmap = 'RdBu', edgecolor = 'white', linewidth = 0.5,legend=True,)
    
11. df_pts.plot(ax = ax, marker = 'o', color = 'limegreen', markersize = 15)
    
12. province =gpd.read_file('/Users/Downloads/下载/25263637/中华人民共和国.shp')
    
13. province.plot(ax = ax,facecolor = 'none', edgecolor = 'dimgrey', linewidth = 0.5,legend=True,)
    
14. # 设置标题与保存图片
    
15. ax.set_title('Rainfall Points & Thiessen Polygons in Loess Plateau', fontdict = {'fontsize': '12', 'fontweight' : '10'})

复制代码

降水点及其泰森多边形插值效果图