在 PostGIS 中创建一个空的平铺栅格表

Question

我是 postgres/postgis 的新手，正在尝试弄清楚如何使用 postgis 制作一个空的平铺栅格。

我想生成一个包含 5000 X 2000 个像元的空栅格，我想稍后查询它以在位置 x/y 处查找特定像元或在位置 x/y 处添加像元值。将来我基本上会将单个单元格值添加到空栅格中（例如，随着时间的推移报告单元格中的动物目击事件）。

我发现通过首先创建一个栅格表可以实现这一点：

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

然后添加一个空栅格：

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

（参考：http://suite.opengeo.org/docs/latest/dataadmin/pgGettingStarted/raster2pgsql.html）

我还添加了一个波段并为其中一个栅格单元添加了一个值：

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)

（参考：https://gis.stackexchange.com/questions/14960/postgis-raster-value-of-a-lat-lon-point）

我现在可以查询我设置的栅格单元的值了：

// 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

我发现很多帖子说，对于较大的栅格，查询性能至关重要，栅格的平铺范围约为每个图块 100X100 个单元格（参考：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/）

现在我不知道：

1. 如何在 PostGIS 中从我的大栅格中创建切片，或者如何使用 PostGIS 从一开始就创建切片栅格？换句话说：我需要使用什么查询来创建一个栅格切片集合，该集合覆盖 5000X2000 栅格，其中切片的切片范围为 100X100 像元，包含多个波段？
2. 创建切片后，我需要创建单独的空间索引还是自动创建？
3. 最后：如何在平铺栅格后查询特定位置的像元值？

任何帮助将不胜感激！

Answer 1

经过反复试验，我似乎找到了解决方案。

要在 PostGIS 中创建一个包含空切片的新表，这些表加起来可以扩展为 5000 X 2000 个单元格，我使用了以下查询：

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

这会使用 ST_MakeEmptyRaster 创建一个范围为 5000 X 2000 的临时栅格。然后，它使用 ST_Tile 将临时栅格平铺成 10X10 单元格平铺，并将这些平铺添加到表格中。

然后我可以添加我的乐队：

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

最后，我可以添加和检索值：

// 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

如您所见，这已将查询时间减少到原始时间的约 2/5。

如果我另外添加一个索引：

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

然后再次执行查询，我得到：

// 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

如您所见，查询时间再次缩短了一半多一点，导致查询时间不到原始查询的 1/5。 我还尝试了不同的磁贴尺寸，发现 10 X 10 的磁贴尺寸会产生最佳性能，但仅比 100 X 100 的磁贴尺寸略好。

如果有人有进一步的优化，欢迎补充。

编辑 (08.07.2016)

我已经写了一篇关于这个的小博客文章。有兴趣的可以看这里：http://www.geonet.ch/postgres-postgis-of-rasters-and-geojsons/