3

I am using Shapely to do a basic intersection operation. However, Shapely doesn't return the intersection line I am looking for. I though the reason might be Shapely's intersection doesn't have ArcGIS's tolerance setting.

As shown in this image, I want to use red 'pass' polygon to clip the green rectangle polygon and get yellow polygon which is the polygon intersection between red and green polygon. Then I'd like to get the intersection line between red 'pass' polygon and yellow interesction polygon. enter image description here

However, as shown in this image, only part of the intersection line is returned. (The red line is the intersection line. The upper-right corner line is missed)

enter image description here

Here is my script:

import shapely;import geopandas as gpd;import os

def Shapely2GDF(ShapelyGeom):
    import pandas as pd
    import geopandas as gpd
    SampleDict = {}
    SampleDict[0] = [ShapelyGeom]
    SampleDF  = pd.DataFrame.from_dict(SampleDict, orient='index', columns=['geometry'])
    SampleGDF = gpd.GeoDataFrame(SampleDF, geometry = 'geometry')
    return(SampleGDF)
folder_path = r"YourFolderPath"
CRS = {'proj': 'laea',
     'lat_0': 33,
     'lon_0': 87,
     'x_0': 0,
     'y_0': 0,
     'datum': 'WGS84',
     'units': 'm',
     'no_defs': True}

Pass_wkt = 'POLYGON ((189758.9127545664 -66191.8185786074, 180329.3792152433 -116707.8439576719, 178965.1189195188 -124016.6660195934, 178800.8029019888 -123985.4142005378, 158503.0483487034 -233805.3173365155, 156868.4312602524 -242649.6541701788, 143913.7618549466 -312744.5093643814, 143758.953929076 -312715.355133106, 140197.0832145372 -332169.2871406805, 138084.9767773161 -332169.2871406805, 127821.3433724809 -332169.2871406805, 107552.6633413704 -332169.2871406805, 89774.74872879591 -332169.2871406805, 95330.73292863928 -301751.8496292587, 95488.81857719552 -301781.133725035, 98283.33938273881 -286624.8484239653, 107910.5678397901 -234411.9796224758, 129528.0923211118 -117176.5374084665, 130316.8087657643 -112899.3632158041, 130484.6345863547 -112930.7501034001, 132069.1199091757 -104422.3853548318, 148632.9652992385 -15482.73408142012, 158222.2461373163 -15482.73408142012, 178599.7741489997 -15482.73408142012, 192824.6413851511 -15482.73408142012, 199225.1582247838 -15482.73408142012, 189758.9127545664 -66191.8185786074))'
Polygon_wkt = 'POLYGON ((172956.2479248047 -124143.5767211914, 188165.8128662109 -120930.9952392578, 190789.2002563477 -133351.0950317383, 175579.6354980469 -136563.6763916016, 172956.2479248047 -124143.5767211914))'

Pass_geom = shapely.wkt.loads(Pass_wkt)
Test_polygon_geom = shapely.wkt.loads(Polygon_wkt)

Pass = Shapely2GDF(Pass_geom); Pass.crs = CRS
Test_polygon = Shapely2GDF(Test_polygon_geom); Test_polygon.crs = CRS

intersection_polygon = gpd.overlay(Pass, Test_polygon, how = 'intersection')
OutFileName = 'Polygon_intersection.shp'
OutFilePath = os.path.join(folder_path, OutFileName)
intersection_polygon.to_file(OutFilePath)

Intersection_geom = intersection_polygon.iloc[0].geometry
#Pass_geom         = Pass.iloc[0].geometry 

Contain_result   = Pass_geom.contains(Intersection_geom)
Intersect_result = Pass_geom.intersects(Intersection_geom) 

CurrentGeometry = Pass_geom; PointList = []
for pt in list(CurrentGeometry.exterior.coords): 
    PointList.append(pt)
Pass_boundary = shapely.geometry.LineString(PointList)

intersection_line = Intersection_geom.intersection(Pass_boundary)
LineType = type(intersection_line)

intersection_line_gdf = Shapely2GDF(intersection_line)
intersection_line_gdf.crs = CRS
intersection_line_gdf.to_file('Result_line.shp')
  • Sorry I miss the below information in my question: The raw input data only have the yellow 'intersection' polygon and red 'pass' polygon, the green one does not exist. To better illustrate where does this yellow 'intersection' polygon come from, I mannually clip it from the green polygon. – ToYoungSekiro Nov 21 '19 at 17:19
2

The way you are getting that intersection line is a bit odd and is prone to precision errors.

What happens is that when you get that yellow intersection polygon, two points are created which were not present initially:

from shapely.wkt import loads
a = loads('POLYGON ((189758.9127545664 -66191.8185786074, 180329.3792152433 -116707.8439576719, 178965.1189195188 -124016.6660195934, 178800.8029019888 -123985.4142005378, 158503.0483487034 -233805.3173365155, 156868.4312602524 -242649.6541701788, 143913.7618549466 -312744.5093643814, 143758.953929076 -312715.355133106, 140197.0832145372 -332169.2871406805, 138084.9767773161 -332169.2871406805, 127821.3433724809 -332169.2871406805, 107552.6633413704 -332169.2871406805, 89774.74872879591 -332169.2871406805, 95330.73292863928 -301751.8496292587, 95488.81857719552 -301781.133725035, 98283.33938273881 -286624.8484239653, 107910.5678397901 -234411.9796224758, 129528.0923211118 -117176.5374084665, 130316.8087657643 -112899.3632158041, 130484.6345863547 -112930.7501034001, 132069.1199091757 -104422.3853548318, 148632.9652992385 -15482.73408142012, 158222.2461373163 -15482.73408142012, 178599.7741489997 -15482.73408142012, 192824.6413851511 -15482.73408142012, 199225.1582247838 -15482.73408142012, 189758.9127545664 -66191.8185786074))')
b = loads('POLYGON ((172956.2479248047 -124143.5767211914, 188165.8128662109 -120930.9952392578, 190789.2002563477 -133351.0950317383, 175579.6354980469 -136563.6763916016, 172956.2479248047 -124143.5767211914))')
c = a.intersection(b)

enter image description here
And, most probably, the higher point appears to be lying not really on the boundary of the red polygon, and, therefore, you are losing a corresponding segment.

line = a.boundary.intersection(c)

enter image description here

Instead, you can get the line by simply calculating the intersection of the boundary of a red polygon with the green polygon:

line = a.boundary.intersection(b)

enter image description here


And if, for some reason, you have only the red polygon (a), and its intersection (c) with the green polygon, but no information about the green polygon itself (b), then you can take one of the boundaries, dilate it using buffer method, and get the intersection line with another boundary:

line = a.boundary.buffer(1e-9).intersection(c.boundary)
# or
line = c.boundary.buffer(1e-9).intersection(a.boundary)

This will produce the same result as in the image above.

  • Thank you for your suggestion! – ToYoungSekiro Nov 21 '19 at 17:07
  • Sorry I miss the below information in my question: The raw input data only have the yellow 'intersection' polygon, the green one does not exist. For the purpose of sharing this issue, it is more convinient to let shapely clip the green one first. Otherwise, if I only provide yellow and red polygons, people may doubt my interesction of clip method. In other word, what should I do if I only have your polygon-a and polygon-c? Assume polygon-b does not exist. – ToYoungSekiro Nov 21 '19 at 17:15
  • @ToYoungSekiro I updated the answer. Please, check. – Georgy Nov 22 '19 at 9:41
  • Thank you so much Georgy! – ToYoungSekiro Nov 22 '19 at 18:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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