[GIS] Using ST_MakeEnvelope over the antimeridian

antimeridianpostgis

I'm using ST_MakeEnvelope in PostGIS to try and make geometries that cross the international dateline (antimeridian), but they are making a shape that wraps the other way around the globe.

The geometry I want is one which covers the eastern half of Australia and the dateline.

It looks something like:

ST_MakeEnvelope(133, -49, -170, 0, 4283)

But:

ST_Intersects(ST_SetSRID(ST_MakePoint(150, -25), 4283), ST_MakeEnvelope(133, -49, -170, 0, 4283));

returns false. The shape goes the wrong way around the world!

If I just take:

ST_Intersects(ST_SetSRID(ST_MakePoint(150, -25), 4283), ST_MakeEnvelope(133, -49, 180,v0, 4283));

I get true, as expected.

How can I set a direction to my envelope? Or should I create a polygon some other way?

(If you're wondering: I'm trying to select geometries in a table for a map view, and when the view includes the dateline, nothing is shown on the map)

[ETA]

If it helps, the box actually arrives in Spherical Mercator, looking something like:

14855319.859171, -6294627.0589708, 21149946.918142, 0

But calling ST_Transform:

st_transform(st_makeenvelope(14855319.859171, -6294627.0589708, 21149946.918142, 0, 900913), 4283)

Results in geometries with the negative-valued longitudes:

SRID=4283;POLYGON((133.447608799776 -49.1163523261112,133.447608799776 0,-170.006794251165 0,-170.006794251165 -49.1163523261112,133.447608799776 -49.1163523261112))

Best Answer

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.

Select the polygons that cross the anti-meridian, cast to geom and shift to 0-360, then split them. Select the split polys with any longitude > 180, and shift back 360 deg

Create table beams_fixed_dump4 as Select ST_Translate(foo2.thegeom, -360, 0) as newgeom from
(Select (ST_DUMP(foo.geom)).geom as thegeom from 
(Select ST_Split(shift_lon.geom, line2) as geom 
    from ST_SetSRID(ST_MakeLine(ST_MakePoint(180, -90),ST_MakePoint(180, 90)), 4326) As line2, 
    (Select ST_ShiftLongitude(geog::geometry) as geom from 
        (SELECT ST_Intersects(poly_coverage.geog, line) as inter, poly_coverage.* from poly_coverage, ST_GeogFromText('LINESTRING(180 -90, 180 0, 180 90)') As line) as tmp where inter = TRUE) as shift_lon) as foo) as foo2 where st_xmax(foo2.thegeom) > 180
Union all

Select the polygon that cross the anti-meridian, cast to geom and shift to 0-360, then split them. Select the split polys with all longitude <= 180.

Select newgeom from 
(Select (ST_DUMP(foo.geom)).geom as newgeom from 
(Select ST_Split(shift_lon.geom, line2) as geom from 
    ST_SetSRID(ST_MakeLine(ST_MakePoint(180, -90),ST_MakePoint(180, 90)), 4326) As line2, 
    (Select ST_ShiftLongitude(geog::geometry) as geom from 
        (SELECT ST_Intersects(poly_coverage.geog, line) as inter, poly_coverage.* 
            from poly_coverage, ST_GeogFromText('LINESTRING(180 -90, 180 0, 180 90)') As line) as tmp where inter = TRUE) as shift_lon) as foo)as foo2 where st_xmax(newgeom) <= 180
Union all

Select the polygons that don't cross the anti-meridian and cast to geom.

SELECT foo.geog::geometry from 
(Select ST_Intersects(poly_coverage.geog, line) as inter, poly_coverage.* from 
poly_coverage, ST_GeogFromText('LINESTRING(180 -90, 180 0, 180 90)') As line) as foo where inter = FALSE;

[GIS] Best-practices for databases and APIs with geographic data spanning the antimeridian

Speaking purely from a data storage and analysis perspective, the geography type for PostGIS was designed with the antimeridian in mind (among several design goals). There are several functions specifically designed for the geography type.

For instance, consider a LineString across Taveuni, Fiji (mapped with Great Circle Mapper), which straddles the antimeridian:

SELECT ST_Length('LINESTRING(179.9 -17, -179.8 -16.7)'::geography);

The length of this geodesic is about 46 km. Similarly, ST_Area would work correctly on a Polygon of the island, even with longitude coordinates jumping between +179 and -179.

Casting a EPSG:4326 geometry to a geography type also normalizes the coordinates, for example the longitude of the last coordinate is >180:

SELECT ST_AsGeoJSON('LINESTRING (179.9 -17, 180.2 -16.7)'::geography);

NOTICE:  Coordinate values were coerced into range [-180 -90, 180 90] for GEOGRAPHY
                           st_asgeojson
------------------------------------------------------------------
 {"type":"LineString","coordinates":[[179.9,-17],[-179.8,-16.7]]}

is converted back to the exact same geography type in the first example, but now with GeoJSON output. You can choose to ignore the NOTICE (or e.g. SET client_min_messages TO WARNING;), and convert all sorts of funny-looking geometries as geography types.

Showing geography types on maps outside PostGIS is a different story, and hopefully better answers will touch on this aspect.

Best Answer

Related Solutions

PostGIS – How to Convert Geography to Geometry and Split Polygons Crossing Antimeridian

Select the polygons that cross the anti-meridian, cast to geom and shift to 0-360, then split them. Select the split polys with any longitude > 180, and shift back 360 deg

Select the polygon that cross the anti-meridian, cast to geom and shift to 0-360, then split them. Select the split polys with all longitude <= 180.

Select the polygons that don't cross the anti-meridian and cast to geom.

[GIS] Best-practices for databases and APIs with geographic data spanning the antimeridian

Related Question