在 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和https://duncanjg.wordpress.com/2013/09/21/effect-of-tile-size-and-data-storage-on-postgis-光栅查询次数/ )

现在我不知道：

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 的图块尺寸稍好一些。

如果大家有进一步优化的地方，欢迎补充。

编辑（2016年7月8日）

我写了一篇关于此的小博文。如果有兴趣，您可以在这里查看：http ://www.geonet.ch/postgres-postgis-of-rasters-and-geojsons/