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;