Create an empty tiled raster table in PostGIS

Question

I am new to postgres/postgis and am trying to figure out how to make an empty tiled raster with postgis.

I would like to generate an empty raster containing 5000 X 2000 cells, which I would like to query later to find a specific cell at location x/y or add a cell value at location x/y. I will basically add individual cell values to the empty raster in the future (e.g. report animal sightings in cells over time).

I found out that most of this this is possible by first creating a raster table:

CREATE TABLE myRasterTable(rid serial primary key, rast raster);

and then adding an empty raster:

INSERT INTO myRasterTable(rid,rast)
VALUES(1, ST_MakeEmptyRaster( 5000, 2000, 2485869.5728, 1299941.7864, 100, 100, 0, 0, 2056) );

(Reference: http://suite.opengeo.org/docs/latest/dataadmin/pgGettingStarted/raster2pgsql.html)

I have also added a band and added a value to one of the raster cells:

UPDATE myRasterTable SET rast = ST_AddBand(rast, 1, '32BF'::text, 0);
UPDATE myRasterTable SET rast = ST_SetValue(rast, 1,ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056),987.654321)

(Reference: https://gis.stackexchange.com/questions/14960/postgis-raster-value-of-a-lat-lon-point)

I can now query the value of the raster cell I have set:

// Location with value
SELECT rid, ST_Value(rast, ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056),false) val FROM myRasterTable
// Return = 987.654296875 in 0.5224609375 Seconds

// Location without value:
SELECT rid, ST_Value(rast, ST_Transform(ST_SetSRID(ST_MakePoint(7.0,48.5),4326),2056),false) val FROM myRasterTable
// Return =  0 in 0.51311993598938 Seconds

I have found numerous posts saying, that for larger rasters, it is essential in regards to query performance, that rasters are tiled with extents of about 100X100 cells per tile (references: http://postgis.17.x6.nabble.com/raster-loading-and-ST-Value-performance-td4999924.html AND https://duncanjg.wordpress.com/2013/09/21/effect-of-tile-size-and-data-storage-on-postgis-raster-query-times/ )

Now I do not know:

1. How to create tiles out of my large raster in PostGIS or how I could create a tiled raster from the start with PostGIS? In other words: What query do I need to use to create a raster tile collection covering a 5000X2000 raster with tiles of an extent of 100X100 cells containing multiple bands?
2. Once the tiles are created, do I need to create a seperate spatial index or is that done automatically?
3. And lastly: how do I then perform a query of what cell value a specific location has after the raster has been tiled?

Any help would be appreciated!

Answer 1

After quite some trial and error, I seem to have found the solution.

To create a new table in PostGIS containing empty tiles which add up to an extend of 5000 X 2000 cells, I used this query:

CREATE TABLE myRasterTable (rid integer, rast raster);
INSERT INTO myRasterTable(rast) VALUES(ST_Tile(ST_MakeEmptyRaster( 5000, 2000, 2485869.5728, 1299941.7864, 100, 100, 0, 0, 2056), 10,10, TRUE, NULL));

This makes a temporary raster with an extent of 5000 X 2000 using ST_MakeEmptyRaster. It then tiles the temporary raster using ST_Tile into 10X10 cell tiles and adds these tiles to the table.

I can then add my band:

UPDATE myRasterTable SET rast = ST_AddBand(rast, 1, '32BUI'::text, 0, NULL);

And finally, I can add and retrieve values with:

// Add cell value
UPDATE myRasterTable SET rast = ST_SetValue(rast, 1,ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056),987654321) WHERE ST_Intersects(rast, ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056));

// Get cell value
SELECT rid, ST_Value(rast, ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056), false) FROM myRasterTable WHERE ST_Intersects(rast, ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056));
// Return = 987654321 in 0.22717499732971 Seconds

As you can see, this already reduced query time to ~2/5 of the original time.

If I additionally add an index:

CREATE INDEX myRasterTable_rast_gist_idx ON myRasterTable USING GIST (ST_ConvexHull(rast));

And then perform the query again, I get:

// Get cell value
SELECT rid, ST_Value(rast, ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056), false) FROM myRasterTable WHERE ST_Intersects(rast, ST_Transform(ST_SetSRID(ST_MakePoint(7.5,48.5),4326),2056));
// Return = 987654321 in 0.084153890609741 Seconds

As you can see, the querytime was again cut in a bit more then half, resulting in a querytime of under 1/5 of the original query. I also tried different tile sizes and found that a tile size of 10 X 10 resulted in the best performance, but only slightly better than a tile size of 100 X 100.

If anyone has further optimisations, feel free to add.

EDIT (08.07.2016)

I have written a small blog post about this. If interested, you can have a look here: http://www.geonet.ch/postgres-postgis-of-rasters-and-geojsons/