[GIS] Determine if point is within an irregular polygon using Python

maskingmatplotlibpolygonpythonshapely

I would like to create a mask of a grid where points are evaluated as either 1's or 0's depending on if they lay outside of an irregular (i.e. not rectangular) polygon.

I have a set of high-latitude coordinates in lat/lon that define the boundary of an area I would like to make a polygon.

From this polygon I would then like to evaluate another set of lat/lon points as to whether they lay within the polygon.

Is there a library best suited for creating such a polygon (given a list of lats/lons that define it)? shapely appears to be geared towards 2D/projected points.

What would be the best method of evaluating other points as to whether they lie within or out of this polygon?

Thanks!

Best Answer

So... my solution was as follows. Thank you Michael for the response and to others who might have other solutions.

import pyproj
from shapely.geometry import Polygon, Point

# WGS84 datum                                                               
wgs84 = pyproj.Proj(init='EPSG:4326')                                       

# Albers Equal Area Conic (aea)                                             
nplaea = pyproj.Proj("+proj=laea +lat_0=90 +lon_0=-40 +x_0=0 +y_0=0 \       
                   +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
# Define Polygon                                                                         
poly_lons = [03.0, 28.0, 28.0, -14.0, 03.0]                                 
poly_lats = [73.0, 73.0, 62.0,  62.0, 68.0] 
grid_lons = [a list of grid longitudes here]
grid_lats = [a list of grid latitudes here]

# Transform polygon and grid coordinates to northern lat AEAC projection    
poly_x, poly_y = pyproj.transform(wgs84, nplaea, poly_lons, poly_lats)      
grid_x, grid_y = pyproj.transform(wgs84, nplaea, grid_lons, grid_lats)

poly_proj = Polygon(zip(poly_x,poly_y))

# Initiate array to store boolean values to                                 
mask = np.zeros_like(grid_lons)                                             

# Step through elements in grid to evaluate if in the study area                                                             
for i in range(x):                                                      
    for j in range(y):                                                  
        grid_point = Point(grid_x[i][j], grid_y[i][j])                  
        if grid_point.within(poly_proj):                                  
            mask[i][j] = 1
Related Question