【问题标题】:Find irregular region in 4D numpy array of gridded data (lat/lon)在网格数据(纬度/经度)的 4D numpy 数组中查找不规则区域
【发布时间】:2015-12-08 02:28:51
【问题描述】:

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

如果我的区域是矩形(或正方形),我知道该怎么做,但是不规则多边形怎么做?

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

【问题讨论】:

  • 你的多边形是如何定义的?
  • 目前多边形只在arcgis中根据地形和观测数据的相关性来定义。覆盖在图像顶部的网格是根据模型 Lat 和 Lon 值创建的,然后偏移,因此模型网格点是上面正方形的中心。
  • 我现在需要在 38x5000x30x20(给或取)的 4D 数组中找到那些网格品脱,如果我将使用 numpy.where 查找东西方和南北之间的索引的区域然后用这些范围对数组进行切片,但我不知道如何最好地找到这些不规则形状。
  • 你能举一个这样的多边形定义的例子(甚至是简化的)吗?它是坐标的大条件函数吗?顶点列表?属于多边形的所有单元格的列表?...它在您的 python 代码中是如何定义的?
  • 例1:绿色小方块 temp[:,:,(lat>=37)&(lat=-122.375) ] 示例2:小的 6 边紫色多边形。我想要顶点 lat=-122.875 和 lat=37.75,lon=37.75,lon>=-122.625 和lat>=37.875,lon>=-122.625 和 lat>=37.875,lon>=-122.875

标签: python arrays algorithm numpy multidimensional-array


【解决方案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

【讨论】:

    【解决方案2】:

    一般类问题称为"Point in Polygon",其中(相当)标准算法基于通过所考虑的点绘制一条测试线并计算它穿过多边形边界的次数(这真的很酷/很奇怪它的工作原理很简单,我认为)。 This is a really good overview which includes implementation information.

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

    • 对于每个区域,形成定义它的所有(纬度/经度)正方形的列表。 根据您的区域是如何定义的,这可能是微不足道的,也可能是烦人的......

    • 1234563 /em> numpy.digitize.
    • 测试该点所在的正方形是否在其中一个区域内。

    如果您仍然遇到问题,请提供有关您的问题的更多详细信息(特别是您的区域是如何定义的)——这样会更容易提供建议。

    【讨论】:

    • 请看上面的评论。谢谢你或你的帮助。
    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 2021-12-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2015-01-25
    • 2017-12-31
    • 2022-10-18
    相关资源
    最近更新 更多