0

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

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']
| improve this 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.

| improve this answer | |

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.