OGR – Get Bounding Box of Geometry Crossing the Antimeridian

antimeridianextentsgdalgeometryogr

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

Best Answer

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.

Related Solutions

Python – Minimum Bounding Circle of Geometry Crossing the 180th Meridian

This is how I solved the problem:

import psycopg2
import re

def replace(matchObj):
    # Group 0 is all groups, group 1 is the first match contained in
    # parentheses. Group 2 is the whole float value.
    value = float(matchObj.group(2))
    value = value - 360
    return matchObj.group(1) + str(value) + " "

def main():
    # connect to postgresql
    conn = psycopg2.connect("conn info")
    cur = conn.cursor()
    cur.execute("select st_astext(the_geom) from tl_2010_us_state10 \
                WHERE name10='Alaska'")

    # Get query results
    rows = cur.fetchall()
    alaskaText = rows[0][0]

    # The regular expression finds a positive float between 170 and 179.9, 
    # inclusive. Matches any character not a negative sign to ensure that the
    # longitude value is positive, then it matches the rest of the number.  
    # To ensure that the number matched is longitude, we check that it is
    # followed by a space (the format being "longitude latitude,".
    matchText = re.compile('([^-])(1[67][0-9]\.[0-9]*) ')

    # Substitute longitude values
    alaskaText = matchText.sub(replace, alaskaText)

    # insert into dabase and commit transaction
    cur.execute("insert into tl_2010_us_state10 (the_geom, name10) \
                VALUES (st_geomfromtext('%s', 4269), 'Alaska2')" %(alaskaText))
    conn.commit()
    conn.close()

if __name__ == "__main__":
    main()

[GIS] Using ST_MakeEnvelope over the antimeridian

For an ordinary grid -180/180 use the geography, just check the command maybe you will need to switch to ST_makebox2D.

Otherwise if dealing with regions around the -180/180 it is better to switch to a 0/360° or consider using another projection.

Best Answer

Related Solutions

Python – Minimum Bounding Circle of Geometry Crossing the 180th Meridian

[GIS] Using ST_MakeEnvelope over the antimeridian

Related Question