如何用Python生成含全部原始点的最粗规则网格？是否值得？

问题概述 #

拥有三组数据：存储经度的x数组、存储纬度的y数组、对应经纬度降雨量的z数组。其中x与y为非均匀间隔，且长度一致（例如x包含34.9912, 35.1568, 35.8881及46, 45, 47这类值）。目标是生成包含所有原始点的最粗规则矩阵，用于测试各类非工具包默认的插值方法，核心需求是尽可能减少插值计算时间。

示例说明 #

假设输入数据：

x = [36,38,31]
y = [12,19,15]
z = [1,2,3]

期望生成矩阵（正值对应原始点位置，其余为0）：

0 0 0 0 0 0 0 2
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
3 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0

现有方案的问题 #

此前找到的scipy.interpolate.interp1d方案存在两个缺陷：

• 无法保证包含所有原始点
• 该函数已被标记为即将弃用，需替换为现代方案

同时，网上推荐的griddata生成的是非均匀网格，不符合规则矩阵要求；直接使用linspace(x[0],x[-1],10000)会生成过密网格，导致全球范围数据处理耗时极高。

Python实现方案 #

方案1：适配示例的基础实现（整数步长场景） #

import numpy as np

# 示例数据
x = np.array([36, 38, 31])
y = np.array([12, 19, 15])
z = np.array([1, 2, 3])

# 1. 确定数据范围
x_min, x_max = x.min(), x.max()
y_min, y_max = y.min(), y.max()

# 2. 设置网格步长（示例为整数步长）
dx = 1
dy = 1

# 3. 生成规则网格节点
x_grid = np.arange(x_min, x_max + dx, dx)
y_grid = np.arange(y_min, y_max + dy, dy)

# 4. 创建全零矩阵
grid_matrix = np.zeros((len(y_grid), len(x_grid)), dtype=np.float64)

# 5. 映射原始点到网格索引（y轴倒序匹配示例矩阵）
x_indices = ((x - x_min) / dx).astype(int)
y_indices = ((y_max - y) / dy).astype(int)

# 6. 填充原始数据
grid_matrix[y_indices, x_indices] = z

# 打印结果
for row in grid_matrix:
    print(' '.join(map(str, row.astype(int))))

方案2：适配高精度真实数据的高效实现 #

针对小数点后四位的真实数据，优先选择覆盖原始数据范围的最小均匀网格，而非全球范围网格，同时可根据插值精度需求调整步长：

import numpy as np

# 模拟真实高精度数据
x = np.array([34.9912, 35.1568, 35.8881, 46.0, 45.0, 47.0])
y = np.array([10.1234, 11.5678, 12.9012, 13.3456, 14.7890, 15.2345])
z = np.array([5.6, 7.8, 9.0, 2.3, 4.5, 6.7])

# 1. 确定数据范围
x_min, x_max = x.min(), x.max()
y_min, y_max = y.min(), y.max()

# 2. 设置网格步长（根据插值精度调整，示例为0.01）
dx = 0.01
dy = 0.01

# 3. 生成规则网格节点
x_grid = np.arange(x_min, x_max + dx, dx)
y_grid = np.arange(y_min, y_max + dy, dy)

# 4. 创建全零矩阵
grid_matrix = np.zeros((len(y_grid), len(x_grid)), dtype=np.float64)

# 5. 映射原始点到最近网格节点（四舍五入+边界裁剪）
x_indices = np.round((x - x_min) / dx).astype(int)
x_indices = np.clip(x_indices, 0, len(x_grid)-1)

y_indices = np.round((y_max - y) / dy).astype(int)
y_indices = np.clip(y_indices, 0, len(y_grid)-1)

# 6. 填充原始数据
grid_matrix[y_indices, x_indices] = z

# 输出网格信息
print(f"网格形状：{grid_matrix.shape}")
print(f"覆盖范围：x[{x_min:.4f}, {x_max:.4f}], y[{y_min:.4f}, {y_max:.4f}]")

方案3：强制原始点精确落在网格节点（自适应步长） #

如果必须保证所有原始点精确对应网格节点，可通过计算坐标差值的最大公约数确定步长：

import numpy as np
import math
from functools import reduce

def compute_gcd(numbers):
    """计算一组数的最大公约数"""
    def gcd(a, b):
        return math.gcd(a, b)
    return reduce(gcd, numbers)

# 模拟真实数据
x = np.array([34.9912, 35.1568, 35.8881, 46.0, 45.0, 47.0])
y = np.array([10.1234, 11.5678, 12.9012, 13.3456, 14.7890, 15.2345])
z = np.array([5.6, 7.8, 9.0, 2.3, 4.5, 6.7])

# 1. 将小数转换为整数（适配小数点后四位精度）
x_int = (x * 10000).astype(int)
y_int = (y * 10000).astype(int)

# 2. 计算坐标差值的最大公约数，转换为步长
x_diff = np.abs(np.diff(np.sort(x_int)))
x_gcd = compute_gcd(x_diff) if len(x_diff) > 0 else 1
dx = x_gcd / 10000

y_diff = np.abs(np.diff(np.sort(y_int)))
y_gcd = compute_gcd(y_diff) if len(y_diff) > 0 else 1
dy = y_gcd / 10000

# 3. 生成规则网格节点
x_min, x_max = x.min(), x.max()
x_grid = np.arange(x_min, x_max + dx, dx)

y_min, y_max = y.min(), y.max()
y_grid = np.arange(y_min, y_max + dy, dy)

# 4. 创建矩阵并填充数据
grid_matrix = np.zeros((len(y_grid), len(x_grid)), dtype=np.float64)
x_indices = ((x - x_min) / dx).astype(int)
y_indices = ((y_max - y) / dy).astype(int)
grid_matrix[y_indices, x_indices] = z

# 输出网格信息
print(f"自适应步长：dx={dx:.4f}, dy={dy:.4f}")
print(f"网格形状：{grid_matrix.shape}")

插值效率优化建议 #

1. 限制网格范围：仅覆盖原始数据的经纬度范围，而非全球范围，可大幅减少网格节点数量
2. 调整步长：根据插值方法的精度要求选择合适步长（如0.01或0.1），避免使用过小步长（如0.0001）
3. 降采样原始数据：对非均匀分布的原始点进行降采样，减少需要处理的数据量

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