4

I've downloaded from the official repository the UTM 5x5 km grid of Spain and uploaded it to a PostgreSQL/PostGIS database (vv. 14 & 3.1.4).

I need to obtain the coordinates of every cell vertices, but my attempts using st_dumppoints() return more than 4 points for each one:

select path, ST_AsText(geom) from (
select (st_dumppoints(geom)).* from referencia.malla_5x5_p) a

{1,1}   "POINT(158179.43709999975 4135320.7804000005)"
{1,2}   "POINT(158147.76989999972 4134822.7887999993)"
{1,3}   "POINT(158144.77589999977 4134822.9791)"
{1,4}   "POINT(158125.79449999984 4134824.1861000005)"
{1,5}   "POINT(157646.78249999974 4134854.6448999997)"
{1,6}   "POINT(157627.35280000046 4134855.8804)"
[...]

From the excellent contribution of @geozelot, I now know that I can't get the "four" vertices of each polygon but how could I get the intersection points of the UTM grid?

1
  • 4
    UTM cells follow the curvature of the given ellipsoid, which is expressed in the data as a varying density of vertices per cell. I.e. these cells are no rectangles, and thus cannot simply be constructed with 4 corners (note that a valid polygon has n + 1 vertices to close the boundary on itself).
    – geozelot
    Commented Apr 28, 2023 at 11:41

1 Answer 1

7

To extract the contextual corner points of each cell, find those vertices that have an angle greater than roughly 90° with respect to the preceding and following vertex:

SELECT
  utm_c.[*],  -- list needed columns explicitly and exclude the geom column
  pts.geom
FROM 
  <utm_cells> AS utm_c
  CROSS JOIN LATERAL (
    SELECT
      dmp.geom,
      ST_Angle(
        LAG(dmp.geom) OVER w,
        dmp.geom,
        COALESCE(LEAD(dmp.geom) OVER w, NTH_VALUE(dmp.geom, 2) OVER w)
      )::INT AS _a
    FROM
      ST_DumpPoints(utm_c.geom) AS dmp
    WINDOW w AS (PARTITION BY id, dmp.path[1:2] ORDER BY dmp.path[3])
  ) AS pts
WHERE
  _a > 3
;

Note that this assumes the cell geometries to be MULTIPOLYGONS; if this is not the case, use:

WINDOW w AS (PARTITION BY id, dmp.path[1] ORDER BY dmp.path[2])

instead!

Here, _a is the angle in radians - rounded down to the nearest integer - between each triplet of vertices, with 3 corresponding to angles greater than ~229.18°.

The result set contains respective corners for every cell, not singular intersection points - you can filter the set with GROUP BY geom, but you loose individual cell information.

2
  • 2
    That's a pretty neat query! One thing, though: the data in question seems to be vectorized pixels, having tiny 'steps' on non-aligned cell borders. Wouldn't this approach fail to exclude them if the 'steps' are at right angles, too?
    – t.ry
    Commented May 10, 2023 at 11:40
  • @t.ry indeed. This would require to apply the same logic iteratively until vertex count is constant...
    – geozelot
    Commented May 10, 2023 at 11:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.