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.

Answer 1

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.

Answer 2

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

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