I have a raster table in PostGIS. When I run ST_SetValue(rast, geom , value) on it to change the pixel value at a point, it works very slowly. Do you know how can I create spatial index for a raster column? (I could not create it using PgAdmin interface with gist method)
-
1Did you tile your raster? It makes a big difference to pixel access speeds and spatial search in general.– John PowellJul 12, 2019 at 16:36
-
1Usually, for rasters, you tile the bounding box of each tile, even though mutliple tiles constitute a complete coverage. Obviously, indexing each pixel is pointless and indexing the bounding box of the raster is pointless, but indexing each sub-tile of the raster will make this a lot quicker, as you have to open much fewer points before you can run ST_SetValue. See this answer, for example.– John PowellJul 15, 2019 at 11:21
-
or this one: gis.stackexchange.com/questions/43053/…– John PowellJul 15, 2019 at 11:25
Add a comment
|
2 Answers
After working more with rasters I want to share my experience with you. If you need to create an spatial index for a raster column in PostGIS, you can simply use the following query:
CREATE INDEX my_index_name_idx ON my_table USING GIST(ST_Envelope(rast));
However, personally, I prefer to add two more geometry columns to the table that has a raster column. Then, add index on those geometry columns.
The first added column (tile_extent) will contain a geometry that represents the extent of the tile which is calculated by ST_Envelope(rast), while the second added column (raster_valued_extent) will contain a geometry that represents the extent of non NoData pixels which is calculated by ST_Polygon(raster rast, integer band_num=1).
SELECT AddGeometryColumn ('my_schema','my_raster_table','tile_extent',4326(or any other SRID),'POLYGON',2);
SELECT AddGeometryColumn ('my_schema','my_raster_table','raster_valued_extent',4326(or any other SRID),'MULTIPOLYGON',2);
UPDATE my_raster_table
SET tile_extent = ST_Envelope(rast);
UPDATE my_raster_table
SET raster_valued_extent = ST_Polygon(rast, 1);
CREATE INDEX my_raster_table_tile_extent_idx ON my_raster_table USING GIST(ST_Envelope(tile_extent));
CREATE INDEX my_raster_table_raster_valued_extent_idx ON my_raster_table USING GIST(ST_Envelope(raster_valued_extent));
Finally, you can create a GIST(geom) index for each added geom column and use those ones in the queries instead of rast column itself. This way, you have two indices and you can choose the right one based on the application.
Then, for example, If you want to crop a raster table with a polygon, 1- if you want to have a cropped raster that also contains those parts of the first raster that had NoData value, and are within the polygon, you use the following query:
SELECT ST_clip(rast, geom, ***)
FROM my_raster_table rt, my_polygon_table pt
WHERE ST_Intersects(rt.tile_extent, pt.geom)
and 2- if you want to have a cropped raster that only contains parts of the first raster that has value other than NoData, you can use the following where clause:
SELECT ST_clip(rast, geom, ***)
FROM my_raster_table rt, my_polygon_table pt
WHERE ST_Intersects(rt.raster_valued_extent, pt.geom)
You may say that we can handle NoData value within ST_Clip() function! I would say yes, but for some applications it is necessary that ST_Intersects() return true only if the valued part of the raster has intersection with other geometry. In these cases, you can use the index that you created on second added column.
So when you are using raster2pgsql raster_options_go_here raster_file someschema.sometable > out.sql, use -I tag, that will create a GiST index on the raster column.
For more go here here
