5

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.

  • how do you download the data...? – ziggy Feb 28 '18 at 20:13
  • 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 1 '18 at 1: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 1 '18 at 8:09
  • 1
    Try to decrease simplifying factor 0.01 of a degree is a few hundred of meters... – Jendrusk Mar 1 '18 at 8:34
3
+50

If you prefer to do it with Python code, 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

  • 1
    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
4

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
  • thanks, this is useful but I prefer not to require QGIS unless absolutely necessary – philshem Mar 19 '18 at 8:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.