Find irregular region in 4D numpy array of gridded data (lat/lon)

Question

I have a large 4-dimensional dataset of Temperatures [time,pressure,lat,lon]. I need to find all grid points within a region defined by lat/lon indices and calculate an average over the region to leave me with a 2-dimensional array.

I know how to do this if my region is a rectangle (or square) but how can this be done with an irregular polygon?

Below is an image showing the regions I need to average together and the lat/lon grid the data is gridded to in the array

Answer 1

I believe this should solve your problem.

The code below generates all cells in a polygon defined by a list of vertices. It "scans" the polygon row by row keeping track of the transition columns where you (re)-enter or exit the polygon.

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