Determine if a polygon can be translated to contain another

Question

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)

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:

• how to fit polygon inside another polygon [closed]: here polygons are allowed to be rotated too

• How to determine if one polygon is contained in another?: translations not allowed. Not in Python.

• How to judge if a polygon is inside another polygon in Python?: translations not allowed

• What algorithms can I use to check if a polygon fits inside of another polygon?: rotation also allowed, unlike in my case.

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.

Answer 1

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()

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:

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