[GIS] How to Clip a PostGIS raster, and get a GeoTIFF out

clippostgisraster

My end goal is to load the DEM of the entire country in the PostGIS Database, and create a service, which will take in a GeoJSON polygon, and give you a tiff file of that area.

For this, my Research is currently stuck in clipping the PostGIS raster table.

This is what I have done so far:

  1. I loaded a small DEM tiff into the PostGIS database, using raster2pgsql. This has created a Table with 16 rows, and it can be successful seen as a Raster in QGIS.

  2. I checked if the Given polygon intersects with the Table, using the following command:

SELECT ST_Intersects(
ST_SETSRID(ST_GeomFromGeoJSON('{ "type": "Polygon", "coordinates": [ [ [ 73.53, 18.57], [ 73.50, 18.45], [ 73.76, 18.48], [ 73.72, 18.58], [ 73.53, 18.57
] ] ] }'), 4326),
rast) from carto_dm1;

This showed me that 2 rows of the table intersect with the given geometry, & 14 rows do not.

If I use ST_Clip, I get multiple records back. I'm using the following command:

Select ST_Clip(rast,
ST_SETSRID(ST_GeomFromGeoJSON('{ "type": "Polygon", "coordinates": [ [ [ 73.53, 18.57], [ 73.50, 18.45], [ 73.76, 18.48], [ 73.72, 18.58], [ 73.53, 18.57] ] ] }'), 4326)
) from carto_dm1
where
ST_Intersects(
ST_SETSRID(ST_GeomFromGeoJSON('{ "type": "Polygon", "coordinates": [ [ [ 73.53, 18.57], [ 73.50, 18.45], [ 73.76, 18.48], [ 73.72, 18.58], [ 73.53, 18.57] ] ] }'), 4326),
rast);

Now how do I get one clipped Raster out, which can be passed to ST_AsTIFF, so that I can write a GeoTiff from it?

Best Answer

One way of doing this, is to use the ST_Union function. It combines all rasters into a single raster. This function is useful if your raster is tiled, or if you loaded several files into a table using raster2pgsql.

This can be the complete Query:

Select ST_AsTIFF(ST_UNION 
                (ST_Clip(rast, 
                        ST_SETSRID(ST_GeomFromGeoJSON('{ "type": "Polygon", "coordinates": [ [ [ 73.534496527777776, 18.573055555555555 ], [ 73.505329861111122, 18.452743055555555 ], [ 73.76782986111111, 18.480086805555555 ], [ 73.729548611111113, 18.585815972222221 ], [ 73.534496527777776, 18.573055555555555 ] ] ] }'), 4326))))
            from carto_dm1
            where 
            ST_Intersects(
                    ST_SETSRID(ST_GeomFromGeoJSON('{ "type": "Polygon", "coordinates": [ [ [ 73.534496527777776, 18.573055555555555 ], [ 73.505329861111122, 18.452743055555555 ], [ 73.76782986111111, 18.480086805555555 ], [ 73.729548611111113, 18.585815972222221 ], [ 73.534496527777776, 18.573055555555555 ] ] ] }'), 4326)
                ,rast);
Related Question