Convert N-sided polygon to nearest fitting rectangle?

Question

I'm working with OSM data and many of my objects are simple 4-polygons (parallelograms), which can be expressed by 4 unique corner points. But the full OSM data has the shape represented by more than 4 points, with multiple points per line segment (see image). This means that the shapes are expressed with redundant points, 16+1 in the case of this polygon:

For an arbitrary polygon, how can I find either

1 the minimum bounding (rotated) rectangle?

2 or the best fitting rectangle?

Example GeoJSON:

{"type":"FeatureCollection","features":[{"type":"Feature","properties":{"osm_id":"1269601","type":"multipolygon","leisure":"pitch","sport":"soccer"},"geometry":{"type":"Polygon","coordinates":[[[6.6131123,46.5124914],[6.6129421,46.5125385],[6.6127998,46.5125783],[6.6126291,46.512626],[6.6125308,46.5124593],[6.6127016,46.5124121],[6.6128452,46.5123724],[6.6130153,46.5123244],[6.6131123,46.5124914]]]}}]}

http://bl.ocks.org/d/465e77bb03acfd7085b2f9a1d2dd08d5

What I've tried:

ogr2ogr's simplify flag and the simplify() function in Shapely, which doesn't force my polygons to be 4 sided, and is also smoothening things I don't want (like corners). It's my understanding that simplify won't work if the points aren't near enough to each other.
Shapely's buffer() function
Bounding Box, which gives me the smallest non-rotated rectangle that fits my polygon.

I get the data from Planet.osm (or some piece of it) and then use osmconvert and osmfilter to keep my tags of interest, and then use ogr2ogr to convert to geojson — philshem, Feb 28 '18 at 20:14
Possible duplicate of Free tool to simplify parcel shapefiles? — tinlyx, Mar 01 '18 at 01:51
I also thought it was a duplicate, but simplify doesn't give me any control to force the shape to be 4 sided — philshem, Mar 01 '18 at 08:09
Try to decrease simplifying factor 0.01 of a degree is a few hundred of meters... — Jendrusk, Mar 01 '18 at 08:34
Yeah, that "works", but as I updated in my question, I get "invalid points" error message. — philshem, Mar 01 '18 at 08:43
You could try minimum bounding rectangle, but that would make many buildings larger. — BERA, Mar 01 '18 at 11:37

score 5 · Answer 1 · answered Mar 18 '18 at 19:51

In QGis, the Processing toolbox has a "Oriented minimum bounding box" algorithm, which does exactly what you want (your first choice). Be careful, you need to have the data saved in the correct coordinate system (your example data is saved in EPSG:4326, even though you visualize it in EPSG:3857, so the stored data is not a rectangle and the algorithm will not give the expected result).

Open QGis and load the data you want (I had to convert your example geojson to EPSG:3857 first)
Open the Processing toolbox (menu Processing -> Toolbox), then search for "Oriented minimum bounding box"
In the small wizard window, choose the layer to apply the algorithm to and where to save the results, there are no other settings

thanks, this is useful but I prefer not to require QGIS unless absolutely necessary — philshem, Mar 19 '18 at 08:12

score 5 · Accepted Answer · edited May 03 '20 at 14:43

You can use minimum_bounding_rectangle() function in Finding minimum-area-rectangle for given points?

For your GeoJSON text, to get "minimum area bounding rectangle (MABR)":

import numpy as np

data = {"type":"FeatureCollection",
 "features":[{"type":"Feature",
              "properties":{"osm_id":"1269601",
                            "type":"multipolygon",
                            "leisure":"pitch",
                            "sport":"soccer"},
              "geometry":{
                  "type":"Polygon",
                  "coordinates":[[
                      [6.6131123,46.5124914],
                      [6.6129421,46.5125385],
                      [6.6127998,46.5125783],
                      [6.6126291,46.512626],
                      [6.6125308,46.5124593],
                      [6.6127016,46.5124121],
                      [6.6128452,46.5123724],
                      [6.6130153,46.5123244],
                      [6.6131123,46.5124914]]
                  ]
              }
             }
            ]
      }
# polygons can have holes, so, ["coordinates"][0] gives you boundary of polygon.
# If you have multipolygon, ["coordinates"][0][0] gives you the first polygon boundary.
geom = data["features"][0]["geometry"]["coordinates"][0]

mabr = minimum_bounding_rectangle(np.array(geom))

# OUT:
#array[[  6.6131123 ,  46.5124914 ],
#      [  6.61306213,  46.51231129],
#      [  6.6125308 ,  46.5124593 ],
#      [  6.61258097,  46.51263941]]

data2 = dict(data) # copy data to data2    
data2["features"][0]["geometry"]["coordinates"][0] = mabr.tolist()

Now, data2 is GeoJSON text with MABR of polygon. But it is always 'great equal' than source polygon. So, you can think of scaling down polygon by rate of source_polygon_area/mabr_area.

I noticed later. Shapely v1.6 includes minimum_rotated_rectangle method which returns minimum area bounding rectangle. — Kadir Şahbaz, Mar 22 '18 at 12:06

Convert N-sided polygon to nearest fitting rectangle?

2 Answers2