0

我正在尝试使用 Scipy interpolation.griddata 执行简单的 2D 线性插值，但它的行为方式很奇怪：它永远运行并且计算不能中断（100% CPU，RAM 不动，我必须杀死过程）。应该插值的数据+坐标位于多边形内。起初，我尝试使用仅包含位于该多边形内的点的网格网格，但它并没有解决问题。

我还尝试减少到​​一个非常小的插值网格（10 * 10），但它一直冻结。

然后我尝试逐个插入网格点，发现只有其中一些导致了问题，我不知道为什么。

这是代码示例：

# Sets negative levels to zero
x_a = x_a.clip(min = 0)
coordinate = np.array(coordinate)

# Experiment zone polygon.
if city == "sf":
    zone = [[-122.397752, 37.799137],
            [-122.423329, 37.795746],
            [-122.418523, 37.772886],
            [-122.39243, 37.793982]]
    zone.append(zone[0]) # to close the polygon.

# Background map extent lat/long coordinates.
x_min_map = min([xy[0] for xy in zone])
x_max_map = max([xy[0] for xy in zone])
y_min_map = min([xy[1] for xy in zone])
y_max_map = max([xy[1] for xy in zone])

def IsPointInsidePolygon(x, y, poly):
    n = len(poly)
    inside = False

    p1x, p1y = poly[0]
    for i in range(n + 1):
        p2x, p2y = poly[i % n]
        if y > min(p1y, p2y):
            if y <= max(p1y, p2y):
                if x <= max(p1x, p2x):
                    if p1y != p2y:
                        xinters = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x
                    if p1x == p2x or x <= xinters:
                        inside = not inside
        p1x, p1y = p2x, p2y

    return inside

# Creates a mesh grid adapted to the sf zone, on which to interpolate.
N = 10
xi = np.linspace(x_min_map, x_max_map, N)
yi = np.linspace(y_min_map, y_max_map, N)
xi, yi = np.meshgrid(xi, yi)
xi, yi = xi.flatten(), yi.flatten()
all_points = np.vstack((xi, yi)).T
delete_index = []
for i, point in enumerate(all_points):
    if not IsPointInsidePolygon(point[0], point[1], zone):
        delete_index.append(i)

grid_points = np.delete(all_points, delete_index, axis = 0)

# Interpolation
zi = interpolate.griddata(coordinate, x_a, grid_points, method = "linear")

“坐标”变量是一个数组，其中包含位于“区域”多边形内的大约 1e4 个浮点坐标。

x_min_map 等是多边形的最大和最小 x、y 值，因此原始网格网格是包含整个多边形的最小矩形。

任何建议将被认真考虑 ！

4

0 回答 0