Get a bounding box of a geometry that crosses the antimeridian using OGR

Question

I'm trying to get a bounding box of the geometry that crosses the antimeridian using OGR:

from osgeo import ogr

geom = ogr.CreateGeometryFromJson('{"type":"MultiPolygon","coordinates":[[[[-180.0,-40.084015],[-180.0,-40.092945],[-179.988885,-40.092945],[-179.988885,-40.084015],[-180.0,-40.084015]]],[[[179.991288,-40.092945],[180.0,-40.092945],[180.0,-40.084015],[179.991288,-40.084015],[179.991288,-40.092945]]]]}')
minX, maxX, minY, maxY = geom.GetEnvelope()

The output is (-180.0, 180.0, -40.092945, -40.084015), which is not what I've expected. According to GeoJSON specification the bounding box should be (179.991288, -40.092945, -179.988885, -40.084015)

Is this a bug in OGR or my expectations are incorrect? Which is the best way to retrieve a bounding box in format defined in GeoJSON specification?

I'm writing code in Python using GDAL 2.3.1

Answer 1

I've got an answer to the github issue from Evan Rouault:

GetEnvelope() computes the envelope in a pure non-geo aware way. It has no notion of SRS / antimeridian / etc. So what you get is the intended behaviour (like some famous chip maker that says their broken CPUs work as intented ;-)). I can see your point that said, but I don't have a solution.

So now I get an envelope of each polygon within a multipolygon, and then calculate the bounding box of the multipolygon based on a set of envelopes.

Answer 2

The easiest way I found to get the bounding box from a geometry that crosses the antimeridian using OGR/GDAL is by using the gdal.VectorTranslate() function, which is the Python binding for OGR2OGR. If you pass the -lco WRITE_BBOX=YES option as a parameter, the output file will contain the correct bounding box.

from osgeo import ogr, gdal
import json

outputFile = "out.geojson"
geom = '{"type":"MultiPolygon","coordinates":[[[[-180.0,-40.084015],[-180.0,-40.092945],[-179.988885,-40.092945],[-179.988885,-40.084015],[-180.0,-40.084015]]],[[[179.991288,-40.092945],[180.0,-40.092945],[180.0,-40.084015],[179.991288,-40.084015],[179.991288,-40.092945]]]]}'

gdal.VectorTranslate(outputFile, geom, options=["-wrapdateline", "-lco", "RFC7946=YES", "-lco", "WRITE_BBOX=YES"]
with open(outputFile) as f:
    jf = json.load(f)
    bbox = jf['features'][0]['bbox']