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?

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);
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.

  • Where do the ids: 131072 & 176148 come from? – Devdatta Tengshe Jul 29 '16 at 9: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 9: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 9:06

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.