在网格数据的 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

Answer 2

一般类别的问题被称为“多边形中的点”，其中（相当）标准的算法基于通过所考虑的点绘制一条测试线并计算它穿过多边形边界的次数（这真的很酷/奇怪的是它的工作原理如此简单，我认为）。 这是一个非常好的概述，其中包括实施信息。

特别是对于“你的问题”，因为你的每个区域都是基于少量的方形单元格定义的 - 我认为更暴力的方法可能会更好。也许是这样的：

对于每个区域，形成定义该区域的所有（纬度/经度）正方形的列表。
• 根据您的区域的定义方式，这可能是微不足道的，也可能是烦人的......

对于您正在检查的每个点，找出它位于哪个方块中。
• 由于方块的行为非常好，您可以使用每个方块的对角手动执行此操作，或使用类似

  numpy.digitize

  的方法。 

测试该点所在的方格是否位于其中一个区域中。
如果您仍然遇到问题，请提供有关您的问题的更多详细信息（具体来说，您的区域是如何定义的）——这样可以更轻松地提供建议。