4

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?

Improve this question
5
  • try plotting both points and polygon in the same figure? cannot figure why coords would flip. Maybe you can share your data?
    – Aman Bagrecha
    Commented Jun 10, 2022 at 19:21
  • Have you tried with the intersects predicate?
    – Mintx
    Commented Jun 10, 2022 at 19:55
  • @AmanBagrecha, you can find oboth code and datasets over here: github.com/miguellobato84/geodata-maps
    – Miguel Ángel Lobato
    Commented Jun 10, 2022 at 23:28
  • @Mintx, yes, I have tried all kind of predicates
    – Miguel Ángel Lobato
    Commented Jun 10, 2022 at 23:28
  • Wouldn't -3.7 40.4 be correct for Madrid? Have you tried to load and project with a different software?
    – M Dollinger
    Commented Jun 11, 2022 at 2:50

1 Answer 1

Reset to default
3

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"]]
Improve this answer
1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.