在网格数据的 4D numpy 数组中查找不规则区域（纬度/经度）

Question

我有一个大型的 4 维温度数据集 [时间、压力、纬度、经度]。我需要找到由纬度/经度索引定义的区域内的所有网格点，并计算该区域的平均值以留下二维数组。

如果我的区域是矩形（或正方形），我知道该怎么做，但是如何用不规则多边形来做到这一点？

下面的图像显示了我需要一起平均的区域以及数据在阵列中网格化的纬度/经度网格

Answer 1

我相信这应该可以解决您的问题。

下面的代码生成由顶点列表定义的多边形中的所有单元格。它逐行“扫描”多边形，跟踪您（重新）进入或退出多边形的过渡列。

def row(x, transitions):
    """ generator spitting all cells in a row given a list of transition (in/out) columns."""

    i = 1
    in_poly = True
    y = transitions[0]
    while i < len(transitions):
        if in_poly:
            while y < transitions[i]:
                yield (x,y)
                y += 1
            in_poly = False
        else:
            in_poly = True
            y = transitions[i]
        i += 1


def get_same_row_vert(i, vertices):
    """ find all vertex columns in the same row as vertices[i], and return next vertex index as well."""

    vert = []
    x = vertices[i][0]
    while i < len(vertices) and vertices[i][0] == x:
        vert.append(vertices[i][1])
        i += 1
    return vert, i


def update_transitions(old, new):
    """ update old transition columns for a row given new vertices. 

    That is: merge both lists and remove duplicate values (2 transitions at the same column cancel each other)"""

    if old == []:
        return new
    if new == []:
        return old
    o0 = old[0]
    n0 = new[0]
    if o0 == n0:
        return update_transitions(old[1:], new[1:])
    if o0 < n0:
        return [o0] + update_transitions(old[1:], new)
    return [n0] + update_transitions(old, new[1:])


def polygon(vertices):
    """ generator spitting all cells in the polygon defined by given vertices."""

    vertices.sort()
    x = vertices[0][0]
    transitions, i = get_same_row_vert(0, vertices)
    while i < len(vertices):
        while x < vertices[i][0]:            
            for cell in row(x, transitions):
                yield cell
            x += 1
        vert, i = get_same_row_vert(i, vertices)
        transitions = update_transitions(transitions, vert)


# define a "strange" polygon (hook shaped)
vertices = [(0,0),(0,3),(4,3),(4,0),(3,0),(3,2),(1,2),(1,1),(2,1),(2,0)]

for cell in polygon(vertices):
    print cell
    # or do whatever you need to do