【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("/mnt/c/Users/Downloads/下载/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. from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter # 经纬度的Formatter
   
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. lon_formatter = LongitudeFormatter(zero_direction_label=False)
    
13. lat_formatter = LatitudeFormatter()
    
14. 
    
15. ax1.set_title('Temperature of Loess Plateau (Ordinary Kriging)',loc='center')
    
16. ax1.set_xticks(x_extent, crs=ccrs.PlateCarree()) #添加边框的经纬度标记
    
17. ax1.set_yticks(y_extent, crs=ccrs.PlateCarree())
    
18. ax1.xaxis.set_major_formatter(lon_formatter)
    
19. ax1.yaxis.set_major_formatter(lat_formatter)
    
20. ax1.set_extent(extents, crs)
    
21. # 绘制海洋、河流、海岸线，需连接网络！！！！！
    
22. ax1.add_feature(cfeat.OCEAN.with_scale('110m'), linewidth=0.5, zorder=1)
    
23. ax1.add_feature(cfeat.RIVERS.with_scale('10m'), linewidth=1,color='skyblue')
    
24. ax1.coastlines()
    
25. # 绘制插值结果（imshow）和站点实际值（scatter，散点颜色、边框、坐标系、与其他图的位置）
    
26. s1 = ax1.imshow(temp1,origin='lower', cmap='RdYlGn_r', transform=crs,extent=[df['Longitude (°)'].min(), df['Longitude (°)'].max(), df['Latitude (°)'].min(), df['Latitude (°)'].max()])
    
27. ax1.scatter(df['Longitude (°)'], df['Latitude (°)'], c=df['Temperature (°C)'], cmap='RdYlGn_r', edgecolors='black', transform=ccrs.PlateCarree(), zorder=2)
    
28. cb1 = fig.colorbar(s1,ax=ax1,orientation='horizontal') # 水平色带图例
    
29. cb1.set_label('°C')
    
30. plt.savefig('/mnt/c/Users/Downloads/下载/25263637/OK6.png',bbox_inches='tight',dpi=300)
    
31. 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')
    
28. plt.savefig('/mnt/c/Users/Downloads/下载/25263637/Temperature IDW.png',bbox_inches='tight')

复制代码

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. from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
   
9. proj =ccrs.PlateCarree()
   
10. coords_df = [list(xy) for xy in zip(df['Longitude (°)'],df['Latitude (°)'])]
    
11. df_pts = gpd.GeoDataFrame(df['Precipitation (mm)'],geometry = [Point(x, y) for x, y in coords_df], crs = proj)
    
12. min_lon_voro, min_lat_voro, max_lon_voro, max_lat_voro = df_pts.buffer(2).total_bounds
    
13. 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('/mnt/c/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'})
    
16. plt.savefig('/mnt/c/Users/Downloads/下载/25263637/Voronoi.png',dpi=300,bbox_inches='tight',)

复制代码

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