GeoPandas – Unable to Execute Spatial Join After CRS Transformation

coordinate systemgeodataframegeopandaspythonspatial-join

I'm trying to join a points GeoPandas GeoDataFrame with a polygon GeoDataFrame, but I'm getting an empty results (inner join).

The points GeoDataFrame

I'm loading it by using latitude and longitude. epsg=4326. If I execute a points.explore() it is shown correctly in the map.

Puerta del Sol A    POINT (40.41721 -3.70183)
Puerta del Sol B    POINT (40.41731 -3.70160)
Miguel Moya     POINT (40.42059 -3.70584)

enter image description here

The polygon GeoDataFrame (districts)

This GeoDataFrame is loaded by reading a .shp file. automatically epsg=25830 is loaded.

When printing this GeoDataFrame after changing CRS, it is inverted

districts_4326 = districts_25830.to_crs(df.crs)

CENTRO  POLYGON ((-3.69316 40.40735, -3.69320 40.40720...
LATINA  POLYGON ((-3.72261 40.41138, -3.72182 40.40851...
CARABANCHEL     POLYGON ((-3.71932 40.39984, -3.71901 40.39983...
USERA   POLYGON ((-3.68315 40.36474, -3.68314 40.36449...

However, even that the points of the polygons are inverted, if I execute a districts_4326.explore() regions are shown correctly.

enter image description here

The join

The problem comes when I try a spatial join, I always get empty result:

points.sjoin(districts_4326, how="inner", predicate = 'overlaps')

I have tried with all kinds of predicate, but still empty results.

I suppose that the problem is because the transformation of districts is returning inverted coordinates, but I don't understand why it's shown correctly in the map.

Any ideas?

Best Answer

Your Points are flipped.

Folium.Map uses (Lat, Long) as arguments, but Shapely Point object uses (Long, Lat). That is the issue here.

From the doc string of Folium.Map

location: tuple or list, default None
    Latitude and Longitude of Map (Northing, Easting).

From the code you shared, flip the "coordinates" as longitude, latitude and things should work fine. Also, use intersection predicate here.

stations = [{
          "type": "Feature",
          "geometry": { "type": "Point", 
                       "coordinates": [float(station["latitude"]), float(station["longitude"])] }, # <---- here
          "properties": station
       } for station in raw["stations"]]
Related Question