1

I want postgis st_clip to output clipped result in extent of initial clipping geometry. Currently, it's crop option can be set to output either in initial raster extent or in extent of intersection between initial raster and clipping geometry.

Here is a test case:

WITH blade AS (
  SELECT
    st_tileenvelope(11,1238, 638) geom
),

butter AS (
SELECT
  st_addband(
  st_makeemptyraster(
    256,
    256,
    st_xmin(t),
    st_ymax(t),
    (st_xmax(t) - st_xmin(t)) / 256,
    - (st_ymax(t) - st_ymin(t)) / 256,
    0,
    0,
    3857),
    '8BUI'::text,
    1) AS rast
FROM
  blade,
  st_translate(geom, 1000, 1000) t
)

SELECT
  st_width(clip),
  st_height(clip),
  st_envelope(butter.rast) AS raster_envelope,
  blade.geom AS blade,
  st_envelope(clip) AS clipped
FROM
  butter,
  blade,
  st_clip(butter.rast, blade.geom, TRUE) clip;

I'm using a random tile as a clipping geometry, then making a test raster by just shifting it by 1000 meters north and east, then clipping it. Resulting raster size is 243 x 243 and is an intersection between clipping geometry and initial raster.

Image below shows the operation: Red square is a clipping geometry, green square is initial raster extent, and hashed area is clipping result extent.

clip result

I'd like a result to be 256 x 256 pixels and it's extent to look like this: wanted clip result

I've tried to recompile postgis after changing this code in rtpg_mapalgebra.c file a bit, but with no luck yet.

  • Which value do you want the geom area not intersecting with the raster be filled with? – Pierre Racine Jun 23 at 0:51
  • Any value that NODATA is set to. – nextstopsun Jun 25 at 17:06
1

Try this:

WITH blade AS (
  SELECT ST_MakeEnvelope(0, 0, 4, 4) geom
), butter AS (
  SELECT ST_AddBand(ST_MakeEmptyRaster(4, 4, 1, 5, 1), '8BUI'::text, 4) rast
), blade_as_raster AS (
  SELECT ST_AsRaster(geom, 4, 4, '8BUI'::text, 2) rast
  FROM blade
)
SELECT ST_MapAlgebra(a.rast, b.rast, '[rast2]', '8BUI'::text, 'FIRST', NULL::text, '[rast1]') rast
FROM blade_as_raster a, butter b;
| improve this answer | |
  • Does this work for multichannel rasters? – nextstopsun Jun 25 at 17:05
  • ST_MapAlgebra() works band per band... – Pierre Racine Jun 26 at 18:46
  • So you really need to call ST_MapAlgebra times number of bands and then union results to a multiband raster. That's quite inefficient. Another issue is that you get pixel alignment error if you call ST_AsRaster on arbitrary geometry (e.g. tile boundaries). – nextstopsun Jun 29 at 8:11
  • You can control the alignment with some variants of ST_AsRaster(). – Pierre Racine Jun 30 at 12:33

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.