[GIS] Postgis: Is there a way to do an intersection query based on lines while taking into account the dateline

antimeridianpostgis

This is an exact copy of a question I asked on Stackoverflow a few weeks back and still haven't found a good answer for it. The original question is from here: original question.

I've been searching about this for a while and can't seem to find anything useful.

In my database, I have an outline of the world coastline saved as lines. I want to take another line and determine where on the coastline that line intersects. My current query is like this:

SELECT st_intersection(makeline(p1.latlong, p2.latlong)::geography, the_geom::geography)
FROM worldcoastline, points p1, points p2 
WHERE p1.id = 1 AND p2.id = 2

All is well and good, and everything seems to line up when I test it in KML and such, except when the dateline is involved. Instead of going across the dateline, it will go all the way around the world, and show me all these collisions it couldn't have had (across Russia, Europe, America). The most interesting thing I find though is if I throw in a st_intersects on the makeline as a geography, it returns false. I'm considering adding it to the where clause, but the problem is lines that cross the dateline and intersect land will return really whacky results anyways!

Is there anything I can do in PostGIS to fix this problem?

EDIT: I just found this in the reference:

ST_Intersection — (T) Returns a geometry that represents the shared
portion of geomA and geomB. The geography implementation does a
transform to geometry to do the intersection and then transform back
to WGS84.

So I guess that's why it's having dateline problems. Any alternatives or suggestions?

EDIT: I just wanted to point out that I did add the st_instersects WHERE clause, but of course if a line does intersect the dateline as well as a coastline I will end up with a messed up ST_Intersection because it will still try to find the intersections going the long way around the world instead of across the dateline

Best Answer

This can be done but it’s not pretty.

The query needs to be broken into 2 parts. One for those lines that do not cross +- 180, another for those that do. Note that you also need to exclude those segments of the coastline that also cross the 180 since they will exhibit the same behavior you see with p1/p2. Also since you are dealing with points and lines you can do everything with geometry objects.

The first part of the query would look similar to what you have now but with the appropriate where clause to exclude all entities (including coastline) that cross 180. Something like this:

SELECT st_intersection(st_makeline(p1.latlong, p2.latlong), the_geom)
FROM worldcoastline, points p1, points p2 
  WHERE p1.id = 1 AND p2.id = 2 and st_length(st_makeline(p1.latlong, p2.latlong) <= 180
AND st_length(the_geom) <= 180

The second part would look like this:

SELECT
CASE 
WHEN st_X(x) > 180 THEN st_translate(x, -360, 0)
            ELSE x
       END
 from
(SELECT st_intersection(st_shift_longitude(st_makeline(p1.latlong, p2.latlong))), st_shift_longitude(the_geom)) as x
FROM worldcoastline, points p1, points p2 
WHERE p1.id = 1 AND p2.id = 2 and st_length(st_makeline(p1.latlong, p2.latlong) > 180
    AND st_length(st_shift_longitude(the_geom)) <= 180)) AS foo`

The second query selects those p1/p2 pairs that cross 180, shifts them by 180 so you know have 0-360 coordinate system. The coastline is also shifted 180. The where clause also excludes those parts of the coastline that cross 0/360 as they now exhibit the wrap around the globe behavior. The CASE clause corrects the longitude of the shifted coordinates back to the standard coordinate system.

UNION the 2 selects and you should have what you want.

Having said that, I generally believe it is asking for trouble doing queries like this in PostGIS. You might be better off re-structuring your data so that simple queries would work. Also doing anything over long distances with geometries usually leads to unpleasant surprises. One thing you could do is turn your p1/p2 pairs into multi lines that do not cross 180. You certainly should do that to your coastline data. Not knowing your data or application it’s hard to have concrete suggestions, but in my experience queries like this are indicative of a data design shortfall.

Related Question