CODEWITHSUNDEEP
pythonpolygongeojsonshapely
Awais Afridi
Awais Afridi

Reputation: 138

Calculating final polygon after subtracting all intersecting polygons using GeoJson - Python

I have a Polygons such as "Jumeirah Islands Clusters" as shown in fig. I want to subtract lake 1, lake 2 GeoJson files from the Islands Cluster. I tried below solution but it is not working to create a final Polygon.

I tried it twice to get the desired result but the output shown is not desired as the lake 2 is not subtracted in the final step.

from shapely.geometry import Polygon, mapping

with open("file_polygon_1.geojson") as f:
    feature = json.load(f)

if feature["geometry"]["type"] == "Polygon":
    polygon_1 = Polygon([(coor[0], coor[1]) for coor in feature["geometry"]["coordinates"][0]])

with open("file_polygon_2.geojson") as f:
    feature_2 = json.load(f)

if feature_2["geometry"]["type"] == "Polygon":
    polygon_2 = Polygon([(coor[0], coor[1]) for coor in feature_2["geometry"]["coordinates"][0]])

new_geometry = mapping(polygon_2.difference(polygon_1)) 

new_feature = dict(type="Feature",id="",properties=dict(Name=""), geometry=dict(type=new_geometry["type"], coordinates=new_geometry["coordinates"]),
)

outjson = dict(type="FeatureCollection", features=[new_feature])

Upvotes: 0

Views: 1045

Answers (1)

Felipe D.
Felipe D.

Reputation: 1271

I think you might be confusing how the difference method/function works. Since I don't have access to your GeoJSON files, I can't really give you a precise answer, so the best I can do is give you an example using my own shapely features:

import shapely 

lake1 = shapely.geometry.Polygon([(0,0),(0,3),(3,3),(3,0)])
lake2 = shapely.geometry.Polygon([(2,2),(2,5),(5,5),(5,2)])
islands = shapely.geometry.Polygon([(1,1),(1,4),(4,4),(4,1)])

The code above will create three different shapely geometries as illustrated below:

shapely features used in example

If you want to generate a final geometry that results in the islands minus the lakes, you can do it in two different ways.

Using set terminology, we can say:

Islands_Minus_Lakes = Islands - Lake1 - Lake2

which is equivalent to:

Islands_Minus_Lakes = Islands - (Lake1 + Lake2)

In python, you can write both of these as follows:

# Version 1 
islands_minus_lakes = (islands.difference(lake1)).difference(lake2)

# Version 2
islands_minus_lakes = (islands.difference(lake1.union(lake2)))

Both approaches yield the exact same results:

Results after difference

Once you have that islands_minus_lakes object, you can choose to manipulate it however you want.

Note about plotting

In case you're wondering, I used the matplotlib and geopandas libraries to plot my results:

import geopandas as gpd

# Creating GeoDataFrames for the plots
gdf = gpd.GeoDataFrame({'name':['lake1','lake2','islands'],
                        'geometry':[lake1,lake2,islands]})

gdf2 = gpd.GeoDataFrame({'name':['islands_minus_lakes'],
                         'geometry':[islands_minus_lakes]})

# Generates the plot of the inputs
gdf.plot(alpha=0.4)

# Generates the plot of the outputs
gdf2.plot(alpha=0.4)

Upvotes: 1

Related Questions