11

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: bad result.

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?

  • Have you tried ST_Segmentize? It also can only work on geometries, but at least seems like a quicker way to generate the points. Anyway, why not just remove points at the pole? It looks more like a side effect of the used projection than a bug. – lynxlynxlynx Jan 27 '15 at 19:05
  • 5
    From your picture, it looks like you have your geometry in EPSG:4326. Antarctica is better suited to a polar stereographic projection such as EPSG:3031. Even then, it looks like you might need to deal with a cut-line to the pole and back, along the Date Line. – Toby Speight Jun 3 '15 at 17:15
3

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.

  • I didn't know the ST_Buffer(geom, 0) trick to eliminate the cut-line - that's useful! – Toby Speight Apr 4 '17 at 15:43

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.