4

I have a raster and a vector table in Postgres 12 PostGIS 3.0

The vector polygon(orange) overlaps the slope raster(grid in black, whiteish)

enter image description here

All I want to do is cut out/erase the orange polygon from the black/white raster.

**My question is just focused on this specific raster and vector but in general this going to be a big process of hundreds of thousands of polygons erasing slope rasters

My process so far:

  1. Convert the polygon to a raster
  2. use MapAlgebra and subtract the layers which is essentially acting as a clip

query

drop table if exists corners.rast_test;
create table corners.rast_test as 
with a as(
SELECT ST_AsRaster(m.shape, (SELECT rast FROM corners.slope_clip LIMIT 1), '32BF', -9999) rast,id
FROM (select m.shape,m.id
from corners.mappluto m 
join corners.corners_final c on st_intersects(m.shape,c.n9x9) and m.id=2019) m
)
SELECT ST_MapAlgebra(b.rast,a.rast, '([rast1]-[rast2])') rast,cornerid
from a, corners.slope_clip  b where cornerid=1085547

the result

enter image description here

So this is the piece of the raster I do not want. I also want to keep the original slope value during this whole process

I found the expression I need to use from this question I posted about the ArcGIS Raster Calculator

Con(IsNull(rast2),rast1)

how do I turn this above expression into the correct st_mapalgebra statement

4
  • 1
    Why don't use st_clip function of postgis to clip the raster
    – Arun
    Commented Dec 11, 2019 at 6:51
  • I need the opposite of clip
    – ziggy
    Commented Dec 11, 2019 at 12:23
  • need to get Pierre Racine for this question! gis.stackexchange.com/users/6362/pierre-racine
    – ziggy
    Commented Dec 11, 2019 at 20:46
  • You want to cut the polygon or the raster??
    – Arun
    Commented Dec 12, 2019 at 5:40

1 Answer 1

1
+50

If I get this right, why not ST_Clip(rast) the raster by the difference of the ST_Envelope(rast) and your Polygon:

SELECT  ST_Clip(rast, ST_Difference(ST_Envelpe(rast), <polygon_geom>), false)
FROM    <raster>,
        <poly_table>
;

This is obviously not optimized; but if running this on a large amount of rasters, it might be the more performant solution to prepare a geometry table with the actual polygons to clip by (the blade, i.e. the difference of the ST_Envelope(rast) and any of your polygons) including the (PK) reference to the respective raster rid and run the ST_Clip(rast) between raster and blade with matching PK/rid.

6
  • This is a good idea, trying it out soon!
    – ziggy
    Commented Dec 12, 2019 at 12:13
  • this got me the answer but it says I need to wait 4 hours until I give you the bounty. thanks
    – ziggy
    Commented Dec 12, 2019 at 14:50
  • sweet! yeah just wait until the the end, maybe someone comes up with a better answer.
    – geozelot
    Commented Dec 12, 2019 at 14:56
  • although this worked do you think there is a way to do this with just rasters? like just using the map algebra?
    – ziggy
    Commented Dec 12, 2019 at 15:01
  • 2
    @ziggy pretty sure, in the end you just want to set null values to some pixels...but the PostGIS raster suite is not my expertise! What I took from my time using it is to avoid touching the actual raster data as much as possible, and work with geometries wherever you can.
    – geozelot
    Commented Dec 12, 2019 at 18:11

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.