I'd like to calculate the "Oriented Bounding Box" (OBB) which is the enclosing rectangle aligned to the longest extent of a point set (which is not necessarily aligned with the axes of the co-ordinate system).
The solution could be similar to Finding minimum-area-rectangle for given points? but I'm looking for a robust solution in PostGIS SQL (and PL/Python if necessary).
This seems to be peanuts for JTS. Let's hope that the method is ported into GEOS so that you can enjoy from it with Python and PostGIS.
JTS has a minimum diameter method tsusiatsoftware.net/jts/javadoc/com/vividsolutions/jts/… that can also directly return "the minimum rectangular Polygon which encloses the input geometry". Convert you point set into multipoint and run the function. But based on your screen capture you may have a polygon to start with, and that you can use directly as input.
I just was pointed to an answer by Rémi-C over at PostGIS mailinglist. He wrote:
"Here is a working solution in plpgsql, designed to be fast to code (and slow to execute :-/)"
That function rc_BboxOrientedFromGeom(geom) does the job in PL/pgsql. Thanks to Rémi-C!
As he said and as I suggested above, most probably there exist more efficient solutions.
Currently, there are two great implementations of this for Python and PostGIS.
The shapely Python library provides the oriented_envelope method.
There's the alias named minimum_rotated_rectangle as well – they operate the same way. It can be used as an attribute of the geometry or function applied to it.
from shapely import Polygon, oriented_envelope
geometry = Polygon([(1, 1), (15, 1), (5, 10), (1, 1)])
# Oriented bounding box
obb = oriented_envelope(geometry)
If you're working with PostGIS, since version 2.5.0 there's the ST_OrientedEnvelope method available.
SELECT ST_OrientedEnvelope(geom) FROM ...
Note that there may be more than one possible enclosing rectangle for a given set of features, although the function will always return the first one found.
References: Shapely docs; PostGIS docs; BostonGIS example
There's a cool solution for this: https://stackoverflow.com/questions/32892932/create-the-oriented-bounding-box-obb-with-python-and-numpy
import numpy as np
from shapely.geometry.base import BaseGeometry
from shapely.geometry import Point, Polygon, MultiPolygon, GeometryCollection
def pca_eigenvectors(pts: np.ndarray) -> np.ndarray:
"""
Returns the principal axes of a set of points.
Method is essentially running a PCA on the points.
Parameters
----------
pts : array_like
"""
ca = np.cov(pts, y=None, rowvar=False, bias=True)
val, vect = np.linalg.eig(ca)
return np.transpose(vect)
def oriented_bounding_box(pts: np.ndarray) -> np.ndarray:
"""
Returns the oriented bounding box width set of points.
Based on [Create the Oriented Bounding-box (OBB) with Python and NumPy](https://stackoverflow.com/questions/32892932/create-the-oriented-bounding-box-obb-with-python-and-numpy).
Parameters
----------
pts : array_like
"""
tvect = pca_eigenvectors(pts)
rot_matrix = np.linalg.inv(tvect)
rot_arr = np.dot(pts, rot_matrix)
mina = np.min(rot_arr, axis=0)
maxa = np.max(rot_arr, axis=0)
diff = (maxa - mina) * 0.5
center = mina + diff
half_w, half_h = diff
corners = np.array([
center + [-half_w, -half_h],
center + [half_w, -half_h],
center + [half_w, half_h],
center + [-half_w, half_h],
])
return np.dot(corners, tvect)
def polygon_from_obb(obb: np.ndarray) -> Polygon:
"""
Returns the oriented bounding box width set of points.
Parameters
----------
obb : array_like
"""
obb = np.vstack((obb, obb[0]))
return Polygon(obb)
This essentially performs a mini Principal Component Analysis to pull out the two perpendicular axes best describing the points, which can then be used to draw. It's all implemented with numpy for the most part so performant for most use cases.
I've pulled that example's implementation and added some utilities for working with shapely and geopandas here: https://github.com/raphaellaude/geo-obb/blob/main/geoobb/obb.py
See example usage: https://github.com/raphaellaude/geo-obb/blob/main/examples/parcel_obbs.ipynb
