Python Coordinates – Handling Multiple Coordinates from Intersection Function

coordinatesintersectionpythonshapely

I am trying to find intersection vertices of two polygons defined below

PG = Polygon([(0.0,0.0), (0.0,4.0), (2.0, 0.0)])
QG = shapely.affinity.rotate(PG, 10.0, origin="centroid")
print(PG.intersection(QG))

It returns a polygon which represents the intersection vertices but I get multiple vertices belonging to polygon PG and QG and not the exact intersection vertices something like this

POLYGON ((0 1.275007557649384, 0 3.233068398870882, 1.176669689726421 1.646660620547158, 1.901423900453083 0.1971521990938349, 0.7833182180345657 0, 0.2248182330207899 0, 0 1.275007557649384))

The graph looks something like this

How to get the exact coordinates of intersection vertices A, B and C as shown above?

Best Answer

As user30184 says, the intersection of two overlapping polygons is a polygon and with your result (light blue polygon)

The intersection of the boundaries of the polygons (linearRings) gives a multipoint in you case

one = LineString(list(PG.exterior.coords))
two = LineString(list(QG.exterior.coords))
print(one.intersection(two).wkt)
MULTIPOINT (0 1.275007557649384, 0 3.233068398870882, 0.2248182330207899 0, 0.7833182180345657 0, 1.176669689726421 1.646660620547158, 1.901423900453083 0.1971521990938349)

These points are also the vertices of the resulting polygon

points = [Point(pt) for pt in list(PG.intersection(QG).exterior.coords)[:-1]]
for i in points:
     print(i.wkt)
POINT (0 1.275007557649384)
POINT (0 3.233068398870882)
POINT (1.176669689726421 1.646660620547158)
POINT (1.901423900453083 0.1971521990938349)
POINT (0.7833182180345657 0)
POINT (0.2248182330207899 0)

But shapely alone does not allow to select the A,B,C points from the others.

Related Solutions

[GIS] Getting intersection of multiple polygons efficiently in Python

One possible approach could be considering the combination of pairs of polygons, their intersections and finally the union of all the intersections via a cascaded union (like suggested here):

from shapely.geometry import Point
from shapely.ops import cascaded_union
from itertools import combinations

circles = [
    Point(0,0).buffer(1),
    Point(1,0).buffer(1),
    Point(1,1).buffer(1),
]

intersection = cascaded_union(
    [a.intersection(b) for a, b in combinations(circles, 2)]
)
print intersection

A more efficient approach should use a spatial index, like Rtree, in order to deal with a lot of geometries (not the case of the three circles):

from shapely.geometry import Point
from shapely.ops import cascaded_union
from rtree import index

circles = [
    Point(0,0).buffer(1),
    Point(1,0).buffer(1),
    Point(1,1).buffer(1),
]
intersections = []
idx = index.Index()

for pos, circle in enumerate(circles):
    idx.insert(pos, circle.bounds)

for circle in circles:
    merged_circles = cascaded_union([circles[pos] for pos in idx.intersection(circle.bounds) if circles[pos] != circle])
    intersections.append(circle.intersection(merged_circles))

intersection = cascaded_union(intersections)
print intersection

[GIS] Geopandas Line Polygon Intersection

When comparing geodataframes with geometry operations in Geopandas, the geometries are first matched by index. In the case where there is no matching index (because you only have a single polygon for instance) then the result will be False.

If it were to compare each object in the GeoSeries you would instead need to get back a full rectangular dataframe of boolean values, and this would likely be very inefficient.

If you do want to compare all geometries then you have two options. The first (and probably easiest) is to use the geopandas sjoin method:

gpd.sjoin(line_gdf, poly_gdf, op='intersects')

This returns a new GeoDataFrame with the geometries for each object on the left dataframe repeated for each geometry they intersect in the right, with the index of the object in the right, i.e.:

                        geometry  index_right
0  LINESTRING (0.5 0.5, 0.7 0.7)            0
1  LINESTRING (0.9 0.9, 0.2 0.6)            0

The second method is to us the pandas apply method on the GeoSeries to return the rectangular dataframe:

line_gdf.geometry.apply(lambda g: poly_gdf.intersects(g))

Which in turn returns (with increasing inefficiency as the dataframes grow):

index_right     0
index_left
0            True
1            True

In general, unless you needed the square matrix, my advice would be to stick to the sjoin method.

Best Answer

Related Solutions

[GIS] Getting intersection of multiple polygons efficiently in Python

[GIS] Geopandas Line Polygon Intersection

Related Question