我是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.htmlhttps://duncanjg.wordpress.com/2013/09/21/effect-of-tile-size-and-data-storage-on-postgis-raster-query-times/)

现在我不知道:
  • 如何从PostGIS的大型栅格中创建切片,或者如何从PostGIS开始创建切片的栅格?换句话说:我需要使用什么查询来创建一个覆盖5000X2000栅格的栅格图块集合,该栅格图块具有100X100像元范围且包含多个波段的图块?
  • 创建图块后,我需要创建一个单独的空间索引还是自动完成?
  • 最后:在对栅格进行平铺后,如何对特定位置具有的像元值进行查询?

  • 任何帮助,将不胜感激!

    最佳答案

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

    为了在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/

    关于postgresql - 在PostGIS中创建一个空的平铺栅格表,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/37910493/

    10-15 18:35