Determine if point is within an irregular polygon using Python

Question

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!

Answer 1

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                                              

Answer 2

Something like this?

import arcpy

for polygon in arcpy.SearchCursor("path/to/polygon file"):
    feature = polygon.shape # Set the base geom to each polygon
    for point in arcpy.SearchCursor("path/to/point file"):
        feature2 = point.shape # Set the comparison geometry to each point
        print "Point %d is within polygon %d: %s" % (point.FID,polygon.FID,feature.contains(feature2))