分組運算作法

今天學生問我分組運算的問題,就順便來。舉個最簡單的例子來說,如下表呈現,目標是以山頭(mountain)和坡向(aspect)為一組,來計算thermophilization indicator,計算的方式是把不同的物種依照海拔分布上下限分成不同等級(rank; R),再乘以在該樣區之覆蓋度(cover; C),接下來將所有物種等級覆蓋度乘積加總,再除以該樣區覆蓋度總和,可用下列公式來表示:

\frac{\sum_{i}^{n} C \cdot R}{\sum_{i}^{n} C}

先來簡化問題,若要計算每一座山頭坡向的物種覆蓋度總和,該如何計算。以下用 R 的 data.table 之 group 概念來實作,也就是 DT[, j, by=group_by]:

# 先使用 data.table fread() 來讀取資料
library(data.table)
dummy <- fread('/path/to/the/file.csv', sep = ',')
# 再依照 DT[, j, by=group] 來運算
dummy[, sum(cover), by = .(mountain, aspect)]

就會得到依山頭、方位的覆蓋度總和

mountain aspect V1
1: SEN W11 40.5
2: SEN W13 33.1
3: SUN E11 21.5
4: SUN S13 51.0

因此,我們可以再進一步來計算 thermophilization indicator

Tsi <- dummy[, sum(cover*r)/sum(cover), by = .(mountain,aspect) ]

這個也就是在資料庫結構化查詢語言(Structural Query Language; SQL)的 group by 概念,SQL 寫法如下:

SELECT 
    mountain,
    aspect,
    sum(cover*r)/sum(cover) as tsi
FROM 
    table_name 
GROUP BY 
    mountain,aspect;

範例檔案及 R code:

mountain aspect species cover r
SEN W11 A 30 1
SEN W11 B 10 3
SEN W11 C 0.5 2
SEN W13 A 15 4
SEN W13 D 5 3
SEN W13 E 3 2
SEN W13 F 0.1 1
SEN W13 G 10 4
SUN E11 A 15 2
SUN E11 B 5 3
SUN E11 C 0.5 2
SUN E11 E 1 4
SUN S13 D 35 1
SUN S13 E 10 1
SUN S13 F 5 2
SUN S13 G 1 4
view raw dummy.csv hosted with ❤ by GitHub
library(data.table)
exampleFile <- 'https://gist.githubusercontent.com/mutolisp/4898bb2e25a33487385bf5ede8e553da/raw/f101cc6e126eb2475d2345377bd02bbee3474799/dummy.csv'
dummy <- fread(exampleFile, sep = ',', header = T)
# calculate cover summation by mountain and aspect
dummy[, sum(cover), by = .(mountain, aspect)]
# calculate thermophilization indicator
Tsi <- dummy[, sum(cover*r)/sum(cover), by = .(mountain,aspect) ]
view raw tsi.r hosted with ❤ by GitHub

將 ESRI Shapefile 轉成 postgis 格式匯入 PostgreSQL 資料庫

可以用 QGIS 的 SPIT 外掛,或是 shp2pgsql 指令,使用 SPIT,圖形介面很簡單,就不用多介紹了

Screen Shot 2012-12-25 at 1.05.25 AM

 

 

用 shp2pgsql 最簡單

shp2pgsql -I -D -c -s 3826 filename.shp > filename.sql  (-s 代表座標系統,常用的是 WGS84 (4326) 或 TWD97 TM2 (3826) )

psql -d database -f filename.sql

利用空間資料庫計算多邊形交集以及計算面積

話說「交集」是個常常被掛在嘴邊的語句,像我和一口巾工讀生們往往都沒有交集,他們不了解我的明白,我也不了解他們的明白啊。交集在幾何上可以用下圖來簡單表示,但是在 GIS 計算交集上,除了空間上的交集外,還有的資料表的交集(ex: join),這個我就不多花篇幅談。Image

以下簡單介紹用空間資料庫的概念來實作多邊形的面積以及交集

一、準備材料:

二、步驟

qgis_layout

上圖是八色鳥在的分布範圍(圖層: range_fairy_pitta,橘色虛線淺綠色底的多邊形區域,簡稱 sp),深綠色是世界的陸域地圖(圖層:map_world,資料屬性表中有國家名稱(cntry_name)以及大陸名稱(continent),簡稱 world),淺橘色則是 sp 和 world 交集的區域。當然我知道做這個很簡單,從 ArcMap 中的工具箱選擇交集工具就可以做出來了,或是用 QGIS 的向量外掛(fTools)的 Geoprocessing Tool/Intersect 即可做出來,但是 GUI 用起來就是有不踏實感 XD 所以我們用空間資料庫來實作!

繼續閱讀 “利用空間資料庫計算多邊形交集以及計算面積"

Postgis 匯入大檔案的 raster

postgis 在 2.0 版以後支援的 raster 的格式,可以透過 raster2pgsql 來轉成 postgis format,但如果 raster 檔案很大,出現記憶體不足的問題時,如

$ psql -d vegetation -f twdtm_asterv2_30m-tm2_twd97.sql 

BEGIN
psql:twdtm_asterv2_30m-tm2_twd97.sql:2: NOTICE: CREATE TABLE will create implicit sequence “twdtm_asterv2_30m-tm2_twd97_rid_seq" for serial column “twdtm_asterv2_30m-tm2_twd97.rid"
psql:twdtm_asterv2_30m-tm2_twd97.sql:2: NOTICE: CREATE TABLE / PRIMARY KEY will create implicit index “twdtm_asterv2_30m-tm2_twd97_pkey" for table “twdtm_asterv2_30m-tm2_twd97″
CREATE TABLE
psql:twdtm_asterv2_30m-tm2_twd97.sql:5: ERROR: out of memory
DETAIL: Cannot enlarge string buffer containing 1073741808 bytes by 8191 more bytes.
CONTEXT: COPY twdtm_asterv2_30m-tm2_twd97, line 1
psql:twdtm_asterv2_30m-tm2_twd97.sql:6: ERROR: current transaction is aborted, commands ignored until end of transaction block
psql:twdtm_asterv2_30m-tm2_twd97.sql:7: ERROR: current transaction is aborted, commands ignored until end of transaction block
ROLLBACK

此時可以啟用 tile 選項 (-t 寬x高) 來裁切成更小的 size ,如 100×100 就是切成寬 100 個、高 100 個小尺寸的 raster: (-s 3826 代表 EPSG:3826 座標系統, -I 代表建立空間索引, -Y 代表用 COPY 的方式取代 INSERT 的 SQL 指令匯入資料 

$ raster2pgsql -s 3826 -t 100×100 -I -Y twdtm_asterv2_30m-tm2_twd97.tif > twdtm.sql

參考資料:

http://postgis.refractions.net/docs/using_raster.xml.html

[HOWTO] 匯入 ESRI Shapefile 至 PostgreSQL 空間資料庫(postgis)

一、安裝

安裝有好幾種方法,而且都很簡單,如果是 *BSD 平台,就用 ports / pkg_add 之類的套件管理系統,GNU/Linux 例如 Debian ubuntu 就用 aptitude / apt-get 之類的,Mac 就用 Macports 或homebrew。Windows 就去抓 package 自己安裝,例如 http://www.postgresql.org/download/windows/,就…自己看!

要安裝的東西有兩個,第一個是 DBMS,即 PostgreSQL,第二個則是 PostGIS,則是讓 PostgreSQL 具有支援空間資料庫的功能。建議安裝 PostgreSQL > 9.1 版,以及 PostGIS 2.0 以上的版本(因為之後設定比較方便,建議用 2.0 以上版本)。另外如果你是使用 EnterpriseDB 公司所出品的安裝程式,就可以在安裝完 PostgreSQL 後利用 Stack Builder 來選擇延伸模組。

二、設定

先建立一個新的資料庫,可以用 pgAdmin ,或是用 PostgreSQL 的 shell 來建立一個名為 test 的資料庫,然後擁有者是你(選項 -O),並且這個資料庫是 utf8 編碼:

# createdb -E utf8 -l en_US.UTF-8 -O yourusername test

接下來建立 postgis extension,先進入 PostgreSQL shell

# psql -d test

建立 postgis extension

test=# CREATE EXTENSION postgis ;

test=# CREATE EXTENSION postgis_topology;

這樣就讓 PostgreSQL 建立 postgis 支援了!

三、匯入 shape file

這個步驟可以選擇用 QGIS 的 SPIT 延伸套件來匯入,

或是用指令的 shp2pgsql 來匯入:

# shp2pgsql -g geom -W big5 -I -s 3826 -D shpfile > table_name.sql

如果轉換都沒有問題的話,就匯入 psql

# psql -d test -f table_name.sql

之後就可以用 QGIS 連 postgis 資料庫了,如下圖:

Image

另外, postgis 也支援 raster 格式,同樣可以用 raster2pgsql 來匯入

raster2pgsql -s 3826 -t 100×100 -I -C -M raster.tif  | psql -d test

psql 卡住了…

1. DBMS 環境:

psql (PostgreSQL) 8.3.3
contains support for command-line editing

2. 想法:

Table 1 (其中,~ 為一些欄位,除了 ba 之外,所有欄位都有資料)

mapno tsuga ba
——-+ ~ +—–+—————–
bigint int double precision

Table 2

mapno ba
——-+ ~ +—————–+
bigint double precision

我想做的事情是,如果 Table 1 中 tsuga 的欄位是某個特定值,例如說 1 好了
就把 Table 2 中的 ba 值進行 log 計算後複製到 Table 1 中的 ba 欄位

繼續閱讀 “psql 卡住了…"