[GIS] Convert N-sided polygon to nearest fitting rectangle

gdalogropenstreetmappolygonpython

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:

enter image description here

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.

Best Answer

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.

Related Question