CODEWITHSUNDEEP
pythonshapefileesrimatplotlib-basemappyshp
Ash
Ash

Reputation: 3550

Given latitude/longitude say if the coordinate is within continental US or not

I want to check if a particular latitude/longitude is within continental US or not. I don't want to use Online APIs and I'm using Python.

I downloaded this shapefile

from shapely.geometry import MultiPoint, Point, Polygon
import shapefile    
sf = shapefile.Reader("cb_2015_us_nation_20m")
shapes = sf.shapes()
fields = sf.fields
records = sf.records()
points = shapes[0].points
poly = Polygon(points)
lon = -112
lat = 48
point = Point(-112, 48)
poly.contains(point)
#should return True because it is in continental US but returns False

The sample lon, lat is within US boundary but poly.contains returns False. I'm not sure what the problem is and how to solve the issue so that I can test if a point is within continental US.

Upvotes: 2

Views: 3107

Answers (2)

adamcircle
adamcircle

Reputation: 724

I ran your code and plotted the polygon. It looked like this:

bad image of the US

If you had run it with this code:

import geopandas as gpd
import matplotlib.pyplot as plt

shapefile = gpd.read_file("path/to/shapes.shp")
shapefile.plot()
plt.show()
# credit to https://stackoverflow.com/a/59688817/1917407

you would have seen this:

enter image description here

So, 1, you're not looking at the CONUS, and 2, your plot is borked. Your code works though and will return True with the geopandas plot.

Upvotes: 0

Ash
Ash

Reputation: 3550

I ended up checking if lat/lon was in every state instead of check in continental U.S., if a point is in one of the states, then it is in continental U.S..

from shapely.geometry import MultiPoint, Point, Polygon
import shapefile
#return a polygon for each state in a dictionary
def get_us_border_polygon():

    sf = shapefile.Reader("./data/states/cb_2015_us_state_20m")
    shapes = sf.shapes()
    #shapes[i].points
    fields = sf.fields
    records = sf.records()
    state_polygons = {}
    for i, record in enumerate(records):
        state = record[5]
        points = shapes[i].points
        poly = Polygon(points)
        state_polygons[state] = poly

    return state_polygons

#us border
state_polygons = get_us_border_polygon()   
#check if in one of the states then True, else False
def in_us(lat, lon):
    p = Point(lon, lat)
    for state, poly in state_polygons.iteritems():
        if poly.contains(p):
            return state
    return None

Upvotes: 1

Related Questions