在Spyder Python中基于.mat数据创建Thiessen多边形的技术咨询

我正在开展一个使用.mat文件数据的项目，目标是利用Thiessen多边形分析数据的分布位置及强度，目前已编写代码计算127×127网格点到9个监测点的欧氏距离，恳请相关实现建议。

现有代码：

for x in range(0,127):
    for y in range(0,127):
        EuD1 = np.sqrt((x - 40)**2 + (y - 30)**2)
        #print(EuD1)
        EuD2 = np.sqrt((x - 20)**2 + (y - 60)**2)
        EuD3 = np.sqrt((x - 50)**2 + (y - 10)**2)
        EuD4 = np.sqrt((x - 60)**2 + (y - 60)**2)
        EuD5 = np.sqrt((x - 80)**2 + (y - 20)**2)
        EuD6 = np.sqrt((x - 20)**2 + (y - 20)**2)
        EuD7 = np.sqrt((x - 115)**2 + (y - 110)**2)
        EuD8 = np.sqrt((x - 90)**2 + (y - 60)**2)
        EuD9 = np.sqrt((x - 20)**2 + (y - 115)**2)

代码优化与核心实现步骤 #

1. 优化距离计算逻辑 #

硬编码9个距离计算既冗余又低效，改用numpy向量化运算批量处理，大幅提升效率且易维护：

import numpy as np

# 存储9个监测点坐标，顺序和原代码的EuD1-EuD9对应
stations = np.array([[40,30], [20,60], [50,10], [60,60], [80,20], [20,20], [115,110], [90,60], [20,115]])
# 生成127×127的网格点矩阵
x_grid, y_grid = np.meshgrid(np.arange(0,127), np.arange(0,127))
# 将网格点转为N×2的数组，方便批量计算
grid_points = np.stack([x_grid.flatten(), y_grid.flatten()], axis=1)
# 批量计算每个网格点到所有监测点的欧氏距离，结果形状为(127*127, 9)
distances = np.sqrt(np.sum((grid_points[:, np.newaxis] - stations)**2, axis=2))

2. 划分Thiessen多边形区域 #

对每个网格点的距离数组取最小值的索引，即可得到该点所属的监测点区域：

# 每个网格点对应的监测点索引（0-8）
region_indices = np.argmin(distances, axis=1)
# 转回127×127的网格形状，方便后续处理
region_grid = region_indices.reshape(127, 127)

3. 加载.mat文件的强度数据 #

用scipy.io.loadmat读取.mat中的监测点强度数据，确保数据顺序和监测点坐标对应：

from scipy.io import loadmat

# 替换为你的.mat文件路径
mat_data = loadmat('your_data_file.mat')
# 假设.mat中存储强度的变量名为station_intensities，转为一维数组
intensities = mat_data['station_intensities'].flatten()

4. 生成强度分布网格 #

根据网格点的区域归属，将对应监测点的强度赋值到网格上：

# 得到整个区域的强度分布网格
intensity_grid = intensities[region_grid]

5. 可视化验证 #

用matplotlib直观展示Thiessen区域划分和强度分布：

import matplotlib.pyplot as plt

# 绘制Thiessen区域划分
plt.figure(figsize=(8,8))
plt.imshow(region_grid, cmap='tab10', origin='lower')
plt.scatter(stations[:,0], stations[:,1], c='red', s=50, label='监测点')
plt.legend()
plt.title('Thiessen多边形区域划分')
plt.show()

# 绘制强度分布
plt.figure(figsize=(8,8))
plt.imshow(intensity_grid, cmap='viridis', origin='lower')
plt.colorbar(label='强度')
plt.scatter(stations[:,0], stations[:,1], c='red', s=50)
plt.title('基于Thiessen多边形的强度分布')
plt.show()

进阶方案：矢量边界生成 #

如果需要Thiessen多边形的矢量边界（而非网格形式），可以用scipy.spatial.Voronoi直接生成：

from scipy.spatial import Voronoi, voronoi_plot_2d

vor = Voronoi(stations)
fig = voronoi_plot_2d(vor, show_vertices=False, line_colors='blue', line_width=2, line_alpha=0.6, point_size=50)
plt.title('Thiessen多边形矢量边界')
plt.show()

内容的提问来源于stack exchange，提问作者Tony