Pobieranie próbek punktowych wzdłuż owijającej linii brzegowej za pomocą PostGIS

11

Pracuję nad zadaniem, które wymaga ode mnie zdobywania punktów próbnych co 1000 km wzdłuż wybrzeży i napotkałem problem z Antarktydą. Z tego, co mogę powiedzieć, wydaje się, że jest to problem z wykorzystaniem geometrii w funkcjach, kiedy naprawdę uważam, że do tej operacji należy użyć geografii .

Korzystanie z funkcji z tym bardzo podobne pytanie , jestem w stanie produkować wynik, który wygląda tak: zły wynik.

Jak widać ST_AddMeasure()i ST_LocateAlong()wydaje się , że nie traktuje geometrii sferycznie, co powoduje, że wiele punktów znajduje się na biegunie południowym. Punkt został nawet dodany do klipu wzdłuż linii daty (po lewej stronie). W dokumentacji tych dwóch funkcji można zastosować tylko geometrię .

Kod użyty do wygenerowania wielokąta i punktów można znaleźć tutaj , ale jest to SQL użyty do wygenerowania punktów:

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;

Jak mogę generować punkty na każde 1000 km wzdłuż tej linii brzegowej?

jczaplew
źródło
Czy próbowałeś ST_Segmentize? Może także działać tylko na geometrii, ale przynajmniej wydaje się szybszym sposobem na generowanie punktów. W każdym razie, dlaczego nie usunąć punktów na słupie? Wygląda bardziej jak efekt uboczny zastosowanej projekcji niż błąd.
lynxlynxlynx
5
Z twojego zdjęcia wygląda na to, że masz geometrię w EPSG: 4326. Antarktyda lepiej nadaje się do polarnej projekcji stereograficznej, takiej jak EPSG: 3031. Nawet wtedy wygląda na to, że będziesz musiał poradzić sobie z linią podziału na słupie iz powrotem, wzdłuż linii daty.
Toby Speight

Odpowiedzi:

3

Jak zasugerowano w jednym z komentarzy, najpierw przekształciłbym geometrię wejściową w polarną projekcję stereograficzną.

Dodatkowo będziesz chciał użyć ST_Buffergo (z pewną ilością 0), aby pozbyć się powstałej linii cięcia.

Dzięki temu uzyskasz pożądany rezultat:

-- 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

Zauważ, że to zapytanie zakłada, że ​​w line_sample_testtabeli jest tylko jeden wiersz , więc dostosuj w razie potrzeby rzeczywiste dane wejściowe.

csd
źródło
Nie znałem ST_Buffer(geom, 0)sztuczki, aby wyeliminować linię podziału - to przydatne!
Toby Speight