Extract the symmetric difference from two polygons?-CodePudding

This code finds the intersection of two polygons and create a polygon out of the intersection. My aim is to do the exact opposite. I would like to create a polygone which has a hole in it in others words it excludes the intersection from the big polygone.

The code i managed to construct is:

from shapely.geometry import shape,Polygon,MultiPolygon,mapping
import geopandas as gpd

g1 = geomap
print(geomap)
g2 = geomap_world
print(geomap_world)

data=[]
for index, orig in g1.iterrows():
    for index2, ref in g2.iterrows():
        if ref['geometry'].intersects(orig['geometry']):                
            data.append({'geometry':ref['geometry'].intersection(orig['geometry'])}) 

df = gpd.GeoDataFrame(data,columns=['geometry'])

where df is :

162  POLYGON ((-2.16991 35.16840, -1.79299 34.52792...

But I am looking to construct the opposite result which is a polygone that excludes the intersection. Any hint?

I also tried symmetric_difference instead of intersection as follows:

data=[]
for index, orig in g1.iterrows():
    for index2, ref in g2.iterrows():
        if ref['geometry'].symmetric_difference(orig['geometry']):
               data.append({'geometry':ref['geometry'].symmetric_difference(orig['geometry'])})     

df = gpd.GeoDataFrame(data,columns=['geometry'])

The result is a polygon that is equal to the intersection i would like to exclude. Where I am going wrong?

CodePudding user response：

The geospatial predicate you need to use in the creation of a polygon with a hole is called difference. Try this code. It plots a country with a hole in it.

import geopandas as gpd
#from shapely.geometry import Polygon    

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
poly = world.geometry.values[12]
hole = world.geometry.values[12].representative_point().buffer(0.5)
poly.difference(hole)

CodePudding user response：

This answer worked for me. it gave me a polygone with a hole in it.

import pandas as pd
import numpy as np
import re    
# Geospatial libraries
import geopandas as gpd    
import glob    
from shapely.geometry import shape,Polygon,MultiPolygon,mapping, LineString, Point    

#geomap contains the polygon that i would like to extract (ref)
geomap = import_shapes_list(path_to_data,shapes_folder,crs='EPSG:4326')

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
# world contains the big polygon (orig)
geomap_world=world[world['name']=='Morocco']


g1=gpd.GeoDataFrame({'geometry': geomap.geometry, 'df1':[1]})
g2=gpd.GeoDataFrame({'geometry': geomap_world.geometry, 'df2':[1]})

data=[]
for index, orig in g1.iterrows():
    for index2, ref in g2.iterrows():
        if ref['geometry'].difference(orig['geometry']):       
            data.append({'geometry':ref['geometry'].difference(orig['geometry'])})    


df = gpd.GeoDataFrame(data,columns=['geometry'])
df.to_file(path_to_output 'difference.geojson', driver="GeoJSON")

Thank you for all your comments and guidance