geopandas doesn't find point in polygon even though it should?

user_7124155 0
Viewed 105

I have some lat/long coordinates and need to confirm if they are with the city of Atlanta, GA. I'm testing it out but it doesn't seem to work.

I got a geojson from here which appears to be legit:

https://gis.atlantaga.gov/?page=OPEN-DATA-HUB

import pandas as pd
import geopandas
atl = geopandas.read_file('Official_City_Boundary.geojson')
atl['geometry']   # this shows the image of Atlanta which appears correct

I plug in a couple of coordinates I got from Google Maps:

x = [33.75865421788594, -84.43974601192079]
y = [33.729117878816, -84.4017757998275]
z = [33.827871937500255, -84.39646813516548]

df = pd.DataFrame({'latitude': [x[0], y[0], z[0]], 'longitude': [x[1], y[1], z[1]]})
geometry = geopandas.points_from_xy(df.longitude, df.latitude)
points = geopandas.GeoDataFrame(geometry=geometry)

points

                     geometry
0  POINT (-84.43975 33.75865)
1  POINT (-84.40178 33.72912)
2  POINT (-84.39647 33.82787)

But when I check if the points are in the boundary, only one is true:

atl['geometry'].contains(points)

0     True
1    False
2    False

Why are they not all true? Am I doing it wrong?

geopandas pandas gis python-3.x
2 Answers
  • I found some geometry similar to what you refer to
  • an alternative approach is to use intersects() to find the contains relationship. NB use of unary_union as the Atlanta geometry I downloaded contains multiple polygons
import pandas as pd
import geopandas
from pathlib import Path
atl = geopandas.read_file(Path.home().joinpath("Downloads").joinpath('Official_City_Council_District_Boundaries.geojson'))
atl['geometry']   # this shows the image of Atlanta which appears correct

x = [33.75865421788594, -84.43974601192079]
y = [33.729117878816, -84.4017757998275]
z = [33.827871937500255, -84.39646813516548]

df = pd.DataFrame({'latitude': [x[0], y[0], z[0]], 'longitude': [x[1], y[1], z[1]]})
geometry = geopandas.points_from_xy(df.longitude, df.latitude)
points = geopandas.GeoDataFrame(geometry=geometry, crs="epsg:4326")

points.intersects(atl.unary_union)

0    True
1    True
2    True
dtype: bool
user_94414040

As it is said in documentation:

It does not check if an element of one GeoSeries contains any element of the other one.

So you should use a loop to check all points.

user_173091080
Related