3

I've numpy meshgrid of lat and lon points with shape 1750 X 1750 and its corresponding data (rainfall) of same shape. I need to find the average of data points that falls inside a polygon,read from a shape file.

If my logic is right, I've to find out the index of lat and lon points that falls inside the polygon and then use that index to filter data and then do the average.

I've successfully implemented it in Matlab using Inpolygon function, but I want rewrite the code in Python. I've used matplotlib Path.contains_points function but it didn't work on numpy ND-array.

Could any one suggest a suitable method ? Your help will be highly appreciated.

pkv
  • 107
  • 1
  • 11
  • 1
    Does this help http://stackoverflow.com/questions/24414230/slicing-numpy-array-with-closed-path – paddyg May 25 '15 at 09:43
  • Similar: http://stackoverflow.com/questions/19413259/efficient-way-to-calculate-distance-matrix-given-latitude-and-longitude-data-in – Christophe Roussy Jul 08 '16 at 09:15

1 Answers1

4

To check if a polygon contains some points, you can simply use matplotlib and more precisely, Path.contains_points implemented in matplotlib.path is the solution. It does accept ND arrays, you just have to flatten them beforehand,

import numpy as np
from matplotlib.path import Path

X, Y = np.meshgrid(x, y)  # X, Y are 2D ndarrays
XY = np.dstack((X, Y))
XY_flat = XY.reshape((-1, 2))

mpath = Path( vertices ) # the vertices of the polygon
mask_flat = mpath.contains_points(XY_flat)
mask = mask_flat.reshape(X.shape)  

Alternatively, you could have a look at GeoPandas module, that has more general implementations for working with geospatial data.

vvvvv
  • 25,404
  • 19
  • 49
  • 81
rth
  • 10,680
  • 7
  • 53
  • 77