# Determine if point is within an irregular polygon using Python

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!

• A quick search with python+point+in+polygon should give over 400 already existing questions and answers with code examples. For instance this here with shapely gis.stackexchange.com/questions/10033/… Commented Dec 2, 2013 at 15:22
• And you should clarify if you want to test if a point is within a given polygon or if you want to generate a hull where all points are within (in your second sentence you talk about creating a polygon?) Commented Dec 2, 2013 at 15:23
• Hmmm... somehow I picked some lousy keywords for searching, as that does yeild quite a lot more. I have a set of coordinates that I would like to create a polygon from. I would then like to test another set of coordinates to see if they fall within this polygon. Commented Dec 2, 2013 at 15:29
• This post illustrates the confounding problem with the keywords on this heikkitoivonen.net/blog/2009/01/26/point-in-polygon-in-python Commented Dec 9, 2013 at 16:26

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

# 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):
``````
• what is a list of grid lat and lon here? Commented Mar 1, 2021 at 7:02

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))
``````
• Does this package require ArcGIS? Commented Dec 2, 2013 at 15:44
• Yes this would require an arcgis installation Commented Dec 2, 2013 at 15:45
• Why choose a solution for which you have to pay even though pure Python solutions exist, as mentioned by Curlew ?
– gene
Commented Dec 2, 2013 at 16:40
• Perhaps then answer it rather than critiquing another person's answer? ;) Commented Dec 2, 2013 at 19:55