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?

Craeft
  • 197
  • 7

1 Answers1

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

.

gene
  • 54,868
  • 3
  • 110
  • 187