Point sampling along a pole-wrapping coastline with PostGIS

Question

I am working on a task that requires me to get sample points every 1000 km along coastlines, and have run into a problem with Antarctica. From what I can tell, it seems to be an issue with the use of geometry in the functions, when really I think geography should be used for this operation.

Using the function from this very similar question, I am able to produce a result that looks like this: .

As you can see, ST_AddMeasure() and ST_LocateAlong() don't seem to treat the geometry spherically, which results in many points that sit on the South Pole. A point was even added on the clip along the date line (left side). Per the documentation of these two functions, only geometry can be used.

The code used to generate the polygon and the points can be found here, but this is the SQL used to generate the points:

CREATE TABLE atest AS WITH line AS 
  (SELECT
      id,
      ST_ExteriorRing((ST_Dump(geom)).geom) AS geom
    FROM line_sample_test),
linemeasure AS
    (SELECT
        ST_AddMeasure(line.geom, 0, (ST_Length(line.geom))::int) AS linem,
    generate_series(0, (ST_Length(line.geom))::int, 10) AS i
FROM line),

geometries AS (
    SELECT
        i,
        ST_LocateAlong(linem, i) AS geom 
    FROM linemeasure)

SELECT
    * from geometries;

How can I generate points at every 1000 Km along this coastline?

Answer 1

As suggested in one of the comments, I would first transform the input geometry to a polar stereographic projection.

Addditionally, you'll want to use ST_Buffer on it (with an amount of 0) to get rid of the resulting cut line.

So this would get your desired result:

-- ST_Transform(geom,3031) reprojects to south polar stereographic,
-- in meters.  ST_Buffer(...) doesn't change the shape, but removes
-- the cut line to the pole (at 180 degrees).
WITH line AS (
    SELECT ST_ExteriorRing(
        ST_Buffer(ST_Transform(geom, 3031), 0)
    ) AS geom
    FROM line_sample_test
),

-- This just generates a table of numbers.  In this case, from 0
-- to the geometry length, counting by 1,000,000 (1000 km).
linemeasure AS (
    SELECT generate_series(0, ST_Length(geom)::int, 1000000) AS i
    FROM line
),

-- Convert those values to a fraction of the overall length (for
-- use as input to ST_LineInterpolatePoint)
linefraction AS (
    SELECT i / ST_Length(geom) AS fraction
    FROM line, linemeasure
),

-- Do the interpolation
geometries AS (
    SELECT ST_LineInterpolatePoint(l.geom, lf.fraction) AS geom
    FROM linefraction lf, line l
),

-- Convert back to EPSG:4326 (i.e. lat/lon coords)
geometries_4326 AS (
    SELECT ST_Transform(geom, 4326) AS geom FROM geometries
)
SELECT * FROM geometries_4326

Note that this query does assume that there's only one row in the line_sample_test table, so tweak as necessary for your real input data.