3

I would like to determine whether a polygon could or could not fit inside another polygon. This latter polygon can be translated arbitrarily. A Python example:

from shapely.geometry import Polygon
import geopandas as gpd

p = Polygon([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]])
q0 = Polygon([[1.5, 1.6], [2, 1.6], [2, 2],[1.5, 1.6]])
q1 = Polygon([[1.4, -0.45], [1.6, -0.45], [1.95, 1.4], [1.4, 1.4], [1.4, -0.2]])
q2 = Polygon([[1.1, -0.05], [1.2, -0.05], [1.2,1.05],[1.1,1.05],[1.1,-0.05]])

gpd.GeoSeries([p,q0,q1,q2]).plot(color=['r','g','b','y'],alpha=0.5)

enter image description here

Clearly, the red square (p0) can be translated to contain the green triangle (q0). The blue trapezoid (q1) and the yellow rectangle (q2) cannot be contained, no matter how we translate the square. (Rotation is not allowed in this example.) Given GeoSeries:

qs = gpd.GeoSeries([q0,q1,q2])

and polygon p, I am looking for function does_it_fit, which, when applied to qs:

qs.geometry.apply(lambda row: does_it_fit(row,p))

returns:

0    True
1    False
2    False
dtype: bool

ie:

def does_it_fit(q,p):
    ### what comes here? ###
    return result

Related but different questions:

This problem was inspired by the interesting task of finding parcels on a map which can contain a building with certain size & orientation. Notebook to play with the above example here.

2
  • Do you want to try each combination e.g. (p,q0), (p,q1), (p,q2), (q0,q1), ... or each q against single p?
    – Taras
    Oct 2, 2022 at 8:51
  • Each q against single p. I provided more qs to explain requirements better (green: fits, blue: too big, yellow: only fits if rotation is allowed, so does not fit in this case).
    – zabop
    Oct 2, 2022 at 9:35

1 Answer 1

2

If rotations are strictly forbidden, what you basically want to do is to translate all your polygons towards the origin. More specifically, by attaching their most lowest-left point to the origin.

After that, it will be trivial to see for each one, if it does overlap another one.

Let's dive into the code:


import numpy as np
from shapely.geometry import Point, Polygon
from shapely.affinity import affine_transform
from shapely.ops import nearest_points
import geopandas as gpd
import matplotlib.pyplot as plt

ref = Point([-10,-10])
p = Polygon([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]])
q0 = Polygon([[1.5, 1.6], [2, 1.6], [2, 2],[1.5, 1.6]])
q1 = Polygon([[1.4, -0.45], [1.6, -0.45], [1.95, 1.4], [1.4, 1.4], [1.4, -0.2]])
q2 = Polygon([[1.1, -0.05], [1.2, -0.05], [1.2,1.05],[1.1,1.05],[1.1,-0.05]])

polys = [p, q0, q1, q2]

gpd.GeoSeries([p,q0,q1,q2]).plot(color=['r','g','b','y'],alpha=0.5)

Here, I simply add some imports and defined a reference point ref which will serve us to find the lowest left point for each polygon. Of course, if you have other polygon, you may want to displace this point elsewhere but setting [-10,-10] is perfectly fine with the data your provided.

Then we define convenience functions:

def find_nearest_point(poly,ref):
    """Find the closest point of a polygon to a reference point."""
    p0, _ = nearest_points(poly, ref)
    p0 = tuple(p0.coords)[0]

    return p0

def transform_affine(poly):
    """Translate the polygon so that it's lowest left point
    is at the origin (0,0)."""
    x,y = find_nearest_point(poly,ref)
    m = [1,0,0,1,-x,-y]
    poly_new = affine_transform(poly, m)

    return poly_new

We can now check that everything went fine by plotting them:

new_polys = [transform_affine(poly) for poly in polys]
g = gpd.GeoSeries(new_polys)
fig, ax = plt.subplots()
g.plot(color=['r','g','b','y'], alpha=0.24, ax=ax)
ax.grid()
ax.set_aspect('equal', 'box')
plt.show()

enter image description here

Things become more obvious now...

But we have to verify it mathematically:

n = len(g)
fig, axs = plt.subplots(n,n, figsize=(15, 15))
for i, gi in enumerate(g):
    gs = gpd.GeoSeries(np.roll(g, shift=i+1))
    for j, gsi in enumerate(gs):
        if gsi.contains(gi) and not gsi.equals_exact(gi,0):
            print(f'iter {i}:\n', f'{gsi.wkt} contains {gi.wkt}')
            xa,ya = gi.exterior.xy
            xb,yb = gsi.exterior.xy
            axs[i,j].plot(xa,ya,'-r', lw=4)
            axs[i,j].plot(xb,yb,'-b', lw=2)
            axs[i,j].grid()

This code is purposefully expanded in order to make the logic clear enough.

Results is as follow:

iter 1:
 POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0)) contains POLYGON ((0 0, 0.5 0, 0.5 0.4, 0 0))
iter 3:
 POLYGON ((0 0, 0.200 0, 0.55 1.85, 0 1.85, 0 0.25, 0 0)) contains POLYGON ((0 0, 0.1 0, 0.1 1.1, 0 1.1, 0 0))

This states which polygon contain which other. I removed self containing polygons with the test and not gsi.equals_exact(gi,0).

And the figure shows every occurrence when it happens, otherwise it's a white plot. You can figure out how to discard those by plotting each individual figure separately when conditions are met. But I like this kind of matrix representation:

enter image description here

This is a pairwise operation. As your problem is apparently a subset of that, I let you figure out the was to discard what's unnecessary to you by, e.g. removing a loop if needed. He who can do more can do less.

The mathematics behind: https://en.wikipedia.org/wiki/Transformation_matrix#Affine_transformations

3
  • 1
    Imagine you wish to check for containment of a small pentagon in a square large enough to contain it. There is no vertex of the polygon that you can translate to any vertex of the square and have the pentagon contained in the square.
    – Llaves
    Oct 3, 2022 at 3:41
  • what is the "lower-left" vertex of a polygon?
    – Llaves
    Oct 3, 2022 at 3:42
  • Cool approach! I agree with @Llaves: I fail to see how this approach covers all possible cases. Or, to give another example: I think above method would give wrong results if the red polygon was Polygon([[0,0],[10,0],[10,2],[8,2],[7,0.1],[0,0]]). Still, nice things to learn from this answer!
    – zabop
    Oct 3, 2022 at 7:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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