ST_SetZ seems to hang indefinitely if a vertex is not in raster

[robe]

Note this example is very close to what is in the manual but with an added vertex not in raster

WITH test_raster AS (
SELECT
ST_SetValues(
  ST_AddBand(
    ST_MakeEmptyRaster(width => 2, height => 2,
      upperleftx => 0, upperlefty => 2,
      scalex => 1.0, scaley => -1.0,
      skewx => 0, skewy => 0, srid => 4326),
    index => 1, pixeltype => '16BSI',
    initialvalue => 0,
    nodataval => -999),
  1,1,1,
  newvalueset =>ARRAY[ARRAY[10.0::float8, 50.0::float8], ARRAY[40.0::float8, 20.0::float8]]) AS rast
)
SELECT
ST_AsText(
  ST_SetZ(
    rast,
    band => 1,
    geom => 'SRID=4326;LINESTRING(1.0 1.9, 1.0 0.2, 10 100)'::geometry,
    resample => 'bilinear'
))
FROM test_raster

After loosing my patience waiting forf an answer, I canceled, and saw tons of these notices:

NOTICE:  Attempting to get pixel value with out of range raster coordinates: (10, -98)
NOTICE:  Attempting to get pixel value with out of range raster coordinates: (10, -98)
NOTICE:  Attempting to get pixel value with out of range raster coordinates: (10, -98)
NOTICE:  Attempting to get pixel value with out of range raster coordinates: (10, -98)

[robe]

I was thinking about what should be done in the event as geometry is not fully covered by a tile.

1) Allow user to have a fill in value for the vertex

or

2) Clip the geometry by the raster and only return that portion of the geometry that is covered by the raster.

The #2 is more useful to me.

I'm thinking about 2 in particular cause the use case I want to use this with is tiled-rasters so doing ST_Clip/ST_Union dance to use this function is actually slower than the old way of just use ST_Value, ST_DumpPoints like so:

Reference data: ​http://gis.ess.washington.edu/data/raster/tenmeter/hawaii/kauai.zip

raster2pgsql -s 26904 -Y -I -C -M kauai/*.bil -t 200x200 ch12.kauai | psql -d postgis_in_action

e.g. old way 10 ms

SELECT 
    ST_AsText(                              
        ST_MakeLine( -- <1>
            ST_Translate( -- <2>
                ST_Force3D((gd).geom),  -- <3>
                0,0,
                COALESCE(ST_Value(rast,(gd).geom),0) -- <4>
            )
        )
    ) As line_3dwkt
FROM 
    (
        SELECT ST_DumpPoints(
            ST_GeomFromText(
                'LINESTRING(
                    444210 2438785,434125 2448785,
                    466666 2449780,466670 2449781
                )',
                26904
            )
        ) As gd
    ) As t -- <5>
    LEFT JOIN 
    ch12.kauai 
    ON ST_Intersects(rast,(t.gd).geom); -- <6>

— 256ms using ST_SetZ

SELECT ST_AsText(
      ST_SetZ(k.rast, geom) -- <1>
  )
FROM (SELECT ST_GeomFromText(
                'LINESTRING(
                    444210 2438785,434125 2448785,
                    466666 2449780,466670 2449781
                )',
                26904
        ) AS geom ) AS t -- <2>
        CROSS JOIN LATERAL 
        (
          SELECT ST_Union( -- <3>
            ST_Clip(rast, 
              ST_Expand(t.geom,10)  -- <4>
            )  
          ) AS rast
             FROM ch12.kauai AS k
            WHERE ST_Intersects(k.rast, t.geom)
      ) AS k;

I haven't tested the timing of clipping the geometry and then unioning the geometry slivers to see if that is faster. I suspect it would be

[robe]

okay that did not work out well at all, not sure why this didn't work:

Same hanging issue

NOTICE:  Attempting to get pixel value with out of range raster coordinates: (200, 63)

SELECT  ST_AsText(
      ST_SetZ(ST_SetBandNoDataValue(k.rast,0),ST_Intersection(t.geom, k.rast::geometry)  ) -- <1>
    )
FROM (SELECT ST_GeomFromText(
                'LINESTRING(
                    444210 2438785,434125 2448785,
                    466666 2449780,466670 2449781
                )',
                26904
        ) AS geom ) AS t 
       INNER JOIN ch12.kauai AS k
           ON ST_Intersects(k.rast, t.geom);

[pramsey]

I think in-filling a default is what's going to happen, I don't think clipping the input is a great idea. If people want clipped, then they can pre-clip.

[komzpa]

Documentation needs a nice example how to quickly add z to a geometry that spans multiple tiles imported by raster2pgsql -t auto.

[Paul Ramsey <pramsey@…>]

In 2851cd19/git:

Avoid doom loop when sampling a vertex outside the range of a geometry in ST_SetZ(). Closes #4898