4

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!

ryanjdillon
  • 433
  • 1
  • 7
  • 19
  • 2
    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 http://gis.stackexchange.com/questions/10033/python-efficiency-need-suggestions-about-how-to-use-ogr-and-shapely-in-more-e – Curlew Dec 02 '13 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?) – Curlew Dec 02 '13 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. – ryanjdillon Dec 02 '13 at 15:29
  • This post illustrates the confounding problem with the keywords on this http://www.heikkitoivonen.net/blog/2009/01/26/point-in-polygon-in-python/ – ryanjdillon Dec 09 '13 at 16:26

2 Answers2

4

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
ryanjdillon
  • 433
  • 1
  • 7
  • 19
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))
michael macdonald
  • 29
  • 1