1

I am attempting to use the second option in Check if a point falls within a multipolygon with Python to check if a given point is on land or water. The shape file is from Natural Earth: https://www.naturalearthdata.com/downloads/10m-physical-vectors/10m-land/

Here is my code:

#!/usr/bin/python3

import shapefile
from shapely.geometry import Point # Point class
from shapely.geometry import shape # shape() is a function to convert geo objects through the interface

#point = (-80,40) # an x,y tuple
point = (40,-80)
myshp = open('./natural_earth/ne_10m_coastline.shp','rb')
mydbf = open('./natural_earth/ne_10m_coastline.dbf','rb')
shp = shapefile.Reader(shp=myshp, dbf=mydbf) #open the shapefile
all_shapes = shp.shapes() # get all the polygons
all_records = shp.records()
for i in range(len(all_shapes)):
    boundary = all_shapes[i] # get a boundary polygon
    if Point(point).within(shape(boundary)): # see if the point is in the polygon
        print('land')
       myshp.close()
mydbf.close()

I used this code on the USG Tiger data set to tell what state an county a point was in. I have tried reversing the order of the coordinates. This code never finds the point over land. Can someone diagnose this?

2

ne_10m_coastline.shp is a Polyline shapefile and not a Polygon shapefile

If you use a Polygon shapefile

With point = (-80,40)

shp = shapefile.Reader("/natural-earth-vector-master/10m_cultural/ne_10m_admin_0_countries_lakes.shp")
all_shapes = shp.shapes() # get all the polygons
all_records = shp.records()
for i in range(len(all_shapes)):
    boundary = all_shapes[i] 
    if Point(point).within(shape(boundary)): #
          name = all_records[i][3]
          print(name)
United States of America

With point = (40,-80)

for i in range(len(all_shapes)):
    boundary = all_shapes[i] 
    if Point(point).within(shape(boundary)): #
         name = all_records[i][3]
         print(name)
Antarctica

.

| improve this answer | |
  • Thanks! As a follow up question, how could you tell it wasn't a polygon file? – Craeft Mar 26 at 11:36
  • 1
    see in pyshp shp.shapeType – gene Mar 26 at 13:46

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.