用 Postgis 抽取特定 raster 檔案的值

假設有兩個 postgis 的圖層,一個是點圖層,另一個是影像(raster)圖層,我們想要用點圖層上去抽取 raster 的數值,這個在 ArcMap 中可以用 Toolbox/Spatial Analyst 中的 extract values to points 或是 sample 來做到,在 GRASS GIS 則是使用 v.sample 來做,但我想介紹用空間資料庫的作法給大家參考:

  1. 前處理,以及匯入 PostgreSQL [1]
    import_from_raster.sh
    # 先使用 gdal 計算坡度,使用 gdaldem,預設的數值是 degree ()
    gdaldem slope -of HFA twdtm_asterv2_30m.img tw_slope_asterv2.img
    # -Y 用 COPY statements 取代原本的 INSERT statement
    # -t 建立圖磚(tiles),這樣處理大檔案的 raster 速度會比較快
    # -I 建立 GIST 空間索引
    # -s 座標系統, EPSG code, 3826 for TWD1997 TM2, while 4326 for WGS84
    raster2pgsql -t auto -Y -I -s 3826 tw_slope_asterv2.img env.tw_slope_asterv2 | psql -d dbname
    # 加入 constraints
    SELECT AddRasterConstraints('env'::name, 'tw_slope_asterv2'::name, 'rast'::name, 'regular_blocking', 'blocksize');
    
  2. 使用 ST_Intersects() 來計算點圖層和影像圖層的交集,並用 ST_VALUE 取出影像圖層中的數值
    sample_raster.sql
    SELECT 
        pt.field1, -- 原本 sample_table 中的 field1
        pt.field2, -- 原本 sample_table 中的 field2
    -- 使用 ST_Value 去取出 sample_table (alias 成 pt) 中 geometry 欄位
    -- 和 raster 交集的數值
        ST_Value(dtm.rast, pt.geom_twd97) elevation 
        ST_Value(slope.rast, pt.geom_twd97) slope 
    FROM
        env.tw_dtm_asterv2 dtm,
        env.tw_slope_asterv2 slope, 
        hc.sample_table pt 
    WHERE 
        st_intersects(dtm.rast, pt.geom_twd97) AND 
        st_intersects(slope.rast, pt.geom_twd97)
    ORDER BY field1;
    
    抽出來的數值就像下表
field1 field2 elevation slope
a 1 3767 19.0181560516357
b 0 3551 16.7224273681641
c 0 3112 15.4987382888794
a 1 3370 45.6003341674805
a 0 3507 13.5721111297607
c 0 3398 2.47793483734131
d 1 3823 11.9162292480469

很容易吧!

不過我想要和 R 做整合(這樣取出來的數值就可以做一些統計的運算),所以我使用 RPostgreSQL package 來把 postgis 取出來的資料存成 DataFrame:

extract_val_fpoints.R
setwd('/path/to/your/working directory')
# 載入 RPostgreSQL (其實 DBI 也是可以做),
# 如果還沒裝的話記得用 install.packages() 安裝
# install.packages('RPostgreSQL', 
#    contriburl = 'http://cran.r-project.org/src/contrib', type="source")
# 這裡的作法是先建立連線,再送 query,之後把 query 的數值用 fetch 抓回來。
library(RPostgreSQL)
conn = dbConnect(PostgreSQL(), user= "_username_", dbname="_dbname_")
# db_query 就是上面提到的 SQL 
db_query = "
SELECT 
        pt.field1, 
        pt.field2, 
        ST_Value(dtm.rast, pt.geom_twd97) elevation 
        ST_Value(slope.rast, pt.geom_twd97) slope 
FROM
        env.tw_dtm_asterv2 dtm,
        env.tw_slope_asterv2 slope, 
        hc.sample_table pt 
WHERE 
        st_intersects(dtm.rast, pt.geom_twd97) AND 
        st_intersects(slope.rast, pt.geom_twd97)
ORDER BY field1;
"
# 使用 dbSendQuery(connection, SQL statements) 查詢資料表,並用 fetch()
# 抓回資料
rs = dbSendQuery(conn, db_query)
rs_table = fetch(rs)

之後所取得的資料,你就可以開始分析或是畫圖了!這樣就不用再開 Desktop GIS 來做,而且速度很快(如果你有建立 postgis spatial index 以及 tiles 的話,你會發現對於愈大的資料運算速度快很多)!

[1] 關於建立 PostgreSQL 的 postgis 空間資料庫支援及匯入 shapefile ,請參考匯入 ESRI Shapefile 至 PostgreSQL 空間資料庫(postgis)
References:

Comments

comments powered by Disqus