PostGIS Raster 空间查询

说明:本例中准备的Raster为dem高程数据,只有一个波段,多波段查询原理相同.完整文件名D:\dem.tiff.

1 栅格数据检查

本章节在cmd中执行,要求先安装gdalinfo.

1.1 栅格空间投影检查

在本例中坐标系统一采用EPSG:4326,因此拿到数据后先检查空间参考是否正确.

gdalinfo -stats "D:/dem.tif"
#如内容太长可输出至文本
#gdalinfo -stats "D:/dem.tif">d:\result.txt

如果Coordinate System is:输出内容和下面相同表示空间参考一至无需再转换了.

Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (100.957279207132430,21.849704436532608)
Pixel Size = (0.000280315066780,-0.000280315066780)

1.2 栅格空间投影坐标系变换

gdalwarp  -t_srs EPSG:4326 -r near -of GTiff "D:/dem.tif" "D:/dem1.tif"
gdalinfo -stats "D:/dem1.tif"
del "D:/dem.tif"
rename "D:/dem1.tif" "dem.tif"
rename "D:/dem1.tif.aux.xml" "dem.tif.aux.xml"

2 栅格数据入库

2.1 入库命令

本章节在cmd中执行,要求先安装psql客户端和postgis客户端.
raster2pgsql使用帮助

raster2pgsql -s 4326 -b 1 -t 8x8  -P -d -C -T "nuts" -X "idxnuts" -Y "D:/dem.tif "  public.dem | psql -h localhost -p 5432 -U postgres -d 数据库名称

重点

• -b:如果是多波段的栅格要全部入库则不要设置这个参数;
• -t:入库后瓦片的大小,这个参数非常非常重要,直接关系到您的查询响应时间,稍候详细讲解;
• -P:根据-t参数自动补全瓦片的值,自动补的值为0x7fff(32767),使用时要注意处理.详细请参看raster2pgsql使用帮助;
• 全用-Y参数可以显著增加入库速度,特别在栅格比较大的情况下;
• 入库的数据暂时不创建空间索引,但主键使用-X指定的索引表空间.

2.2 -t参数详解

• 栅格不要直接存到数据库,否则只有一行数据,加载查询都非常慢.要使用-t参数切片后入库
• -t 设置的值是入库后瓦片的宽高(以栅格像素宽高为单位,对应的函数是ST_Width和ST_Height).
• -t 设置的值一般为2的倍数(2*2,4*4,8*8,32*32,64*64,128*128,256*256,512*512),具体多大需要根据栅格的范围确定,如果范围非常大就设置为512*512或1024*1024先入库,然后根据需求使用相同的空间查询做对比响应时间,以确定-t使用的值.在本例中的dem.tif因为范围比较小,经过反复查询比较,最终选择了8*8.

2.3 确定-t参数的值

本章节在psql中执行.

2.3.1 手动创建空间索引

drop index gidx_dem_geom;
create index gidx_dem_geom on dem using gist(st_convexhull(rast))  with (fillfactor=100) tablespace idxnuts;
vacuum  freeze verbose  analyze dem;

2.3.2 执行sql

在本例中的需求为输入坐标点然后输出此位置高程,因此查询sql如下:

explain (analyze,verbose,costs,buffers,timing) 
select 
       ST_Value(t1.rast,1,t2.point) 
   from dem as t1
   inner join (select (ST_SetSRID(ST_Point(100.998680,21.814713),4326)) as point) as t2 on st_convexhull(t1.rast)~t2.point

• 在第2.1节中设置-t 8x8,至少执行5次sql,分别记录下sql的执行时间
• 重新从第2.1开始,设置-t 2x2,-t 4x4,-t 16x16,-t 32x32,-t 64x64,记录下这个sql的执行时间.
• 最终根据根据sql的执行时间选择一个比较理想的-t值

3 处理-P自动补全的数据

ST_Value函数有多个重载版本,注意根据自己的需求选择,这里使用的是第二个重载的版本.

double precision ST_Value(raster rast, geometry pt, boolean exclude_nodata_value=true);

double precision ST_Value(raster rast, integer band, geometry pt, boolean exclude_nodata_value=true);

double precision ST_Value(raster rast, integer x, integer y, boolean exclude_nodata_value=true);

double precision ST_Value(raster rast, integer band, integer x, integer y, boolean exclude_nodata_value=true);

注意:第三个和第四个重载版本使用的像素单位,像素的值不能超过ST_Width和ST_Height的范围.

with cte(val) as(
    select 
        ST_Value(t1.rast,1,t2.point) 
    from dem as t1
    inner join (select (ST_SetSRID(ST_Point(100.998680,21.814713),4326)) as point) as t2 on st_convexhull(t1.rast)~t2.point
)select (case val when 32767 then --0x7f77
    null
else
    val
end) from cte;