pythongeospatialgeopandasshapelyfiona

Which multipolygon does the point of longtitude and latitude belong to in Python?


I have the longitude and latitude and the expected result is that whichever multipolygon the point is in, I get the name or ID of the multipolygon.

import geopandas as gpd

world = gpd.read_file('/Control_Areas.shp')
world.plot()

Output

   0    MULTIPOLYGON (((-9837042.000 6137048.000, -983...
   1    MULTIPOLYGON (((-11583146.000 5695095.000, -11...
   2    MULTIPOLYGON (((-8542840.287 4154568.013, -854...
   3    MULTIPOLYGON (((-10822667.912 2996855.452, -10...
   4    MULTIPOLYGON (((-13050304.061 3865631.027, -13.

Previous attempts: I have tried fiona, shapely and geopandas to get that done but I have struggled horribly to make progress on this. The closest I have gotten is the within and contains function, but the area of work I have struggled is the transformation of multipolygon to polygon successfully as well and then utilising the power of within and contains to get the desired output.

The shapefile has been downloaded from here.


Solution

  • world.crs gives {'init': 'epsg:3857'} (Web Mercator projection) so you should first reproject your GeoDataFrame in the WGS84 projection if you want to keep the latitude-longitude coordinate system of your point.

    world = world.to_crs("EPSG:4326")
    

    Then you can use the intersects method of GeoPandas to find the indexes of the Polygons that contain your point.

    For example for the city of New York:

    from shapely.geometry import Point
    NY_pnt = Point(40.712784, -74.005941)
    world[["ID","NAME"]][world.intersects(NY_pnt)]
    

    which results in:

    ID  NAME
    20  13501   NEW YORK INDEPENDENT SYSTEM OPERATOR
    

    you can check the result with shapely within method:

    NY_pnt.within(world["geometry"][20])
    

    If you have multiple points, you can create a GeoDataFrame and use the sjoin method:

    NY_pnt = Point(40.712784, -74.005941)
    LA_pnt = Point(34.052235, -118.243683)
    points_df = gpd.GeoDataFrame({'geometry': [NY_pnt, LA_pnt]}, crs='EPSG:4326')
    results = gpd.sjoin(points_df, world, op='within')
    results[['ID', 'NAME']]
    

    Output:

    ID  NAME
    0   13501   NEW YORK INDEPENDENT SYSTEM OPERATOR
    1   11208   LOS ANGELES DEPARTMENT OF WATER AND POWER