说明:本例中准备的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;