2

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?

Devdatta Tengshe
  • 41,311
  • 35
  • 139
  • 263

2 Answers2

2

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);
nronnei
  • 568
  • 1
  • 4
  • 17
Devdatta Tengshe
  • 41,311
  • 35
  • 139
  • 263
1

I have exported a GeoTiff from a clipped raster in Postgis with the following workflow and the use of large objects:

1.Creating a large object:

WITH buffer AS (
SELECT
    id,
    ST_Buffer(geom,50) AS geom
FROM polygon_table),

clip AS (
SELECT 
    partnumber,
    ST_Union(ST_Clip(a.rast, 1,b.geom, TRUE)) AS rast
FROM aspect a
    JOIN buffer b
    ON ST_Intersects(a.rast, b.geom)
        GROUP BY partnumber)

SELECT 
    oid, 
    lowrite(lo_open(oid, 131072), png) As num_bytes
    FROM (
    VALUES (lo_create(0),
           ST_AsGDALRaster(((SELECT ST_Union(rast) FROM clip)),'GTiff')
       )) As v(oid,png)

2.Exporting it to a GTiff file (e.g. with oid=176148)

SELECT lo_export(176148, '/tmp/demo_rast.tif')

(3. Deleting the large object: SELECT lo_unlink(176148))

The code is from a question I have asked.

Here are some more information.

Community
  • 1
Stefan
  • 4,414
  • 1
  • 33
  • 66
  • Where do the ids: 131072 & 176148 come from? – Devdatta Tengshe Jul 29 '16 at 09:00
  • This is a good approach. You can also do it with a pythonu function that saves the bytes directly. – John Powell Jul 29 '16 at 09:05
  • The 131072 is a fix id to export raster files. The 176148 is the oid of the last query (see part 1). It is changing every time you call a query. This oid will shown in the result set of the last part of the query (see part 1) – Stefan Jul 29 '16 at 09:06