Scaling rotated rectangle in Shapely

Question

I have a Polygon which looks like this

It's a minimum-oriented bounding box for the LineString geometry contained within it, I generated it using QGIS. Actually, there are many more of these, but let's use this one as an example.

This polygon is part of a georeferenced Shapefile (the CRS is EPSG:32632, WGS 84/UTM zone 32N), and as WKT it looks like

POLYGON ((505317.38000976015 4344246.695509137, 505316.6889838772 4344247.481812444, 505316.7475176216 4344247.533253576, 505317.43854350504 4344246.746950269, 505317.38000976015 4344246.695509137))

Points are ordered as such: top-left, top-right, bottom-right, bottom-left, top-left (to close the geometry)

I want to pad it on the sides so that I have something like:

where the blue box is ultimately what I want to obtain (it's not a precise drawing of course but I hope it gives a clear enough idea)

So, I tried using shapely's affine transformations of rotation, scaling and skewing, since I also know the angle of rotation of the original green box. Here's my code.

from shapely import affinity

def scale_obb(geometry, angle, yfactor, xskew, yskew):
    rotated_a = affinity.rotate(geometry, -angle)
    scaled_obb = affinity.scale(rotated_a, yfact=yfactor)
    rotated_b = affinity.rotate(scaled_obb, angle)
    skewed_obb = affinity.skew(rotated_b, xs=xskew, ys=yskew) 

where angle here is exactly 138.69008 and yfactor, xskew and yskew are variables. The idea is to rotate the polygon such that it is horizontal, then scale it on the y-axis, then rotate it back and if needed, correct the skewness. However, for example using a yfactor = 3, xskew = -13, and yskew = 0, I get:

which is kind of close to what I want, but not quite. How can I straighten this blue box to look sort of like the rectangle I want? I have tinkered a bit with yfactor, xskew, and yskew but haven't found the right values so far. Maybe I need different transformations?

Answer 1

This is perhaps to specific to work on all your polygons, but it is working on your one example polygon.

Extract the two long side lines, buffer them at one side, union with your polygon:

import shapely
import geopandas as gpd #For plotting

g = shapely.wkt.loads("POLYGON ((505317.38000976015 4344246.695509137, 505316.6889838772 4344247.481812444, 505316.7475176216 4344247.533253576, 505317.43854350504 4344246.746950269, 505317.38000976015 4344246.695509137))")

def splitatvertices(bound):
    """A function to return the two longest line segments of an input line"""
    global coords, lines
    coords = [c for c in bound.coords]
    lines = sorted([shapely.geometry.LineString([c1,c2]) for c1,c2 in zip(coords, coords[1:])], key=lambda x: x.length, reverse=True)
    return lines[:2] #Return the two long sides

segments = splitatvertices(bound=g.boundary) #A list of the two longer side lines
segments = [s.buffer(distance=0.1, single_sided=True) for s in segments] #Buffer them at one side, luckily the outside is buffered
union = shapely.ops.unary_union(segments+[g]) #Union/dissolve the original polygon and the buffered sides

#PLOT
s1 = gpd.GeoSeries(g)
s2 = gpd.GeoSeries(segments)
s3 = gpd.GeoSeries(union.boundary)

ax = s1.plot(color="blue", alpha=0.5) #The original polygon is blue
s2.plot(ax=ax, color="green", alpha=0.5)  #The buffered long sides are green
s3.plot(ax=ax, hatch=r"| |", edgecolor="black") #The union of them is black hatched