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!

  • 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 gis.stackexchange.com/questions/10033/… – Curlew Dec 2 '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 2 '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 2 '13 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 – ryanjdillon Dec 9 '13 at 16:26
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                                              
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))
  • 1
    Does this package require ArcGIS? – ryanjdillon Dec 2 '13 at 15:44
  • Yes this would require an arcgis installation – michael macdonald Dec 2 '13 at 15:45
  • 2
    Why choose a solution for which you have to pay even though pure Python solutions exist, as mentioned by Curlew ? – gene Dec 2 '13 at 16:40
  • 2
    Perhaps then answer it rather than critiquing another person's answer? ;) – michael macdonald Dec 2 '13 at 19:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.