Python determine the mean and extract maximum value inside a polygon over a grid in an Array

Question

I have 3 NumPy arrays which consist of UTM-X(256) and UTM-Y(256) coordinates, and the accumulated Rainfall(65536) for a Weather Radar 256x256 (km) in UTM.

I also have a Polygon inside the Grid bounds that is a Catchment Boundary in UTM.

I need to determine the Average Rainfall over just the catchment polygon (a clipped sub set of the RADAR data), and the maximum, and the location of the maximum. I have already determined the average over the entire RADAR grid.

So the question is: How do I perform analysis on a subset of a NumPy array that is determined by the Polygon? I would have thought that this would be a very common operation, but have not found any Python scripts to perform this operation.

Here is an illustration of the data set:

Answer 1

Here is an outline of a possible approach.

First find the polygon that bounds the catchment boundary. Presuming you know which of the UTM coordinates of your full set of points form that catchment boundary, say it's like this,

catchment = an np.array of (UTM_X, UTM_Y) point tuples

you could find the boundary of that point set using scipy.spatial.ConvexHull

boundary= scipy.spatial.ConvexHull(catchment)

Next, for your array of rainfall data, you would have to test whether the coordinates fall inside or outside of the boundary of the convex hull.

This previous SO question has some good answers explaining ways to do that coordinate test.

Finally you would gather those rainfall data points that passed the test of being inside the boundary and perform whatever statistical calculations you want to do with appropriate NumPy/SciPy statistical functions.

Answer 2

Assuming the boundary is given as a list of the polygon vertices, you could have matplotlib generate a mask for you over the data coordinates and then use that mask to sum up only the values within the contour.

In other words, when you have a series of coordinates that define the boundary of the polygon that marks the region of interest, then have matplotlib generate a boolean mask indicating all the coordinates that are within this polygon. This mask can then be used to extract only the limited dataset of rainfall within the contour.

The following simple example shows you how this is done:

import numpy as np
from matplotlib.patches import PathPatch
from matplotlib.path import Path
import matplotlib.pyplot as plt

# generate some fake data
xmin, xmax, ymin, ymax = -10, 30, -4, 20 
y,x = np.mgrid[ymin:ymax+1,xmin:xmax+1]
z = (x-(xmin+xmax)/2)**2 + (y-(ymin + ymax)/2)**2
extent = [xmin-.5, xmax+.5, ymin-.5, ymax+.5]
xr, yr = [np.random.random_integers(lo, hi, 3) for lo, hi
         in ((xmin, xmax), (ymin, ymax))] # create a contour

coordlist = np.vstack((xr, yr)).T  # create an Nx2 array of coordinates
coord_map = np.vstack((x.flatten(), y.flatten())).T # create an Mx2 array listing all the coordinates in field

polypath = Path(coordlist)
mask = polypath.contains_points(coord_map).reshape(x.shape) # have mpl figure out which coords are within the contour

f, ax = plt.subplots(1,1)
ax.imshow(z, extent=extent, interpolation='none', origin='lower', cmap='hot')
ax.imshow(mask, interpolation='none', extent=extent, origin='lower', alpha=.5, cmap='gray')
patch = PathPatch(polypath, facecolor='g', alpha=.5)
ax.add_patch(patch)
plt.show(block=False)
print(z[mask].sum())  # prints out the total accumulated

In this example, x and y represent your UTM-X and UTM-Y dataranges. z represents the weather rainfall data, but is in this case a matrix, unlike your single-column view of average rainfall (which is easily remapped onto a grid).

In the last line, I've summed up all the values of z that are within the contour. If you want the mean, just replace sum by mean.