使用 GNU parallel 來平行運算

通常我們在使用一些簡單的 shell 工具時,都只會用到一個 CPU 核心,但對多核心的機器來說可以使用 GNU parallel 來善用所有的資源,並減少執行時間。以下舉簡單的例子來說明如何使用 GNU parallel:

檔案

very_big_file.csv,tab 分隔的檔案,裡頭是 GBIF 下載物種的資訊,這個文字檔有 34 GB

使用程式

parallel
grep
測試平台: Mac OS X 10.10, 4 cores (8 threads), 16 GB ram

測試結果

Test A

只用 grep 來擷取 "Acanthiza katherina" 的記錄

$ time grep "Acanthiza katherina" very_big_file.csv
366.08s user 16.88s system 99% cpu 6:25.09 total

總共花了 385 秒左右

Test B

$ cat Chordata_small.csv | parallel --pipe -u grep \"Acanthiza katherina\"
0.82s user 30.26s system 8% cpu 6:16.76 total

總共花了 377 秒左右,疑? parallel 似乎也沒有快到哪邊。因為這個例子中瓶頸是輸出入(從硬碟讀取資料然後丟給 grep 擷取),所以試看看每次將 100 M 行的資料轉送給 grep 同時處理:

Test C

$ cat Chordata_small.csv | parallel --pipe --block 100M -u grep \"Acanthiza katherina\"
478.78 s user 135.86 s system 486% cpu 2:06.71 total

改用 block=100M 之後,僅需要 127 秒,節省了 66% 左右的時間。

結論

如果處理的東西不複雜,也許可以用 parallel 搭配 awk/sed/grep 來做平行運算,畢竟寫個 script 只要幾分鐘,相較寫 MPI 或是 Hadoop 這類比較複雜的平行運算系統,parallel 是 CP 比較高的選擇。

Vim 正規表示式(1)

這個教學中,主要使用 vim 來示範簡單的正規表示式(regular expression),如果不懂可以看一下 wikipedia 上面的正規表示式。看不懂?沒關係,簡單的說就是用某些匹配字元(某些符號,例如"*", "|", "[, ]" 等來找到輸入文字的特徵。

舉實際的例子來看好了,假設我們要用 latex 來排版,需要加索引的時候,通常都會在要加索引的字後面加上 \index{keywords},例如:

Nulla ut vestibulum \index{vestibulum} dolor. Aenean aliquet nisl et quam 
suscipit, eget porta libero molestie. Fusce eget nisl ac urna maximus laoreet. 
Donec bibendum iaculis orci id tristique. Quisque nunc quam, ultrices id sapien 
et, hendrerit \index{hendrerit} pharetra tortor. Class aptent taciti sociosqu 
ad litora torquent per conubia nostra, per inceptos himenaeos. Aliquam ut massa 
nec enim pharetra varius. Aenean nec arcu id sem pretium congue id ut orci. 
Curabitur ultricies in augue ut dignissim. Nulla facilisi.

但是如果你是一份文件,懶得手動加入索引時(假設需要加入特徵的字有 65535 筆時,一個一個加,不就做到天荒地老?),像下面這個例子中,我想要把「Lorem ipsum 福祿壽」當成索引:

\name[1. ]{Lorem ipsum 福祿壽}
\name[2. ]{Quisque nunc 過新年}
\name[3. ]{Vestibulum sit 來爬山}
...
\name[65535. ]{Curabitur augue 千里馬}

用肉眼判斷,需要加索引值的字元都是放在大括號內({}),所以這就是特徵,如果我們要用 vim 要怎麼做呢?在介紹完整作法之前,先來介紹一下 vim 的編輯模式和指令模式。vim 在預設進入編輯器時,是指令模式,要進入編輯模式需要按 i (insert),要跳出編輯模式進入指令模式則是按 esc 鍵 + : (冒號)。而正規表示式則是需要在指令模式中下指令。下面一步步解釋如何做:

  1. 了解如何取代(為了方便解釋,在指令下方有 1,2,3, ... 數字來指示字元位置)

    :%s/A/B/g
    1 2 3 4 5 
    

    1 (:) 代表進入指令模式

    2 (%s) 代表是所有行號(%)取代(s)

    3 A 代表某個特定 pattern

    4 B 代表若符合 A pattern 時,取代為 B

    5 g 代表全域(也就是在整行中所有符合的匹配規則(pattern)都會被取代掉),若沒有標示 g ,則只會取代掉第一個

    例如要把下面一段文字(取自Flora of Taiwan)中的 "lanceolate" 取代成 "oblong":

    Culms slender, to 70 cm tall, 1.5 mm in diam. Blade to 15 cm long,
    (2-)5-15 mm wide; ligule chartaceous, minute, less than 1 mm long,
    truncate. Panicle open, to 30 cm long, branches 4-6-nate, verticillate.
    Spikelets 1-flowered, 3 mm long; glumes unequal, lower glume narrowly
    lanceolate, 1.8 mm long, chartaceous, 1-nerved; upper glume lanceolate,
    3-nerved, 2.5 mm long, chartaceous; lemma subcoriaceous, lanceolate,
    3-nerved, acute, 3 mm long; palea subcoriaceous, linear-lanceolate,
    acute, 2-keeled,as long as spikelet.

    vim 的指令如下:

    :%s/lanceolate/oblong/g

    輸入上述指令則會將所有 lanceolate 取代成 oblong

  2. 範圍
    如果要表示範圍,可以用中括號"[ ]"來包住要表示的範圍,通則可以這樣寫

    [n-m]

    其中 n 代表起始的值,m 代表結束的值,這些值可以是英數字,例如


    [1-9][0-9]

    上面這個例字代表符合10, 11, 12, ... , 99 的數字都會被抓出來,範圍也可以是用","(逗號)做分隔,例如

    [1,3-5]

    代表符合 1 或是 3, 4, 5 數字的匹配規則都會被找出來。如果是英文字的話,可以使用```[A-Z]```代表符合 A 到 Z 所有的字元都會被抓出來,大小寫是有區別的。也可以使用其他的正規表示式符號來搭配,例如"*"代表符合 0 到多個符合的所有字元等,但如果是要符合所有的字元包括空白行則是要用".*",通常是使用包在括號、引號的字元匹配規則內,細節請參見 vim regular expression 的文件。
    
    再來看一個例子,如果我要把上面 Flora of Taiwan 的文字中,有數字單位的,例如 10 cm, 1.5 mm 等的匹配規則 抓出來,如下:
    

    Culms slender, to 70 cm tall, 1.5 mm in diam. Blade to 15 cm long,
    (2-)5-15 mm wide; ligule chartaceous, minute, less than 1 mm long,
    truncate. Panicle open, to 30 cm long, branches 4-6-nate, verticillate.
    Spikelets 1-flowered, 3 mm long; glumes unequal, lower glume narrowly
    lanceolate, 1.8 mm long, chartaceous, 1-nerved; upper glume lanceolate,
    3-nerved, 2.5 mm long, chartaceous; lemma subcoriaceous, lanceolate,
    3-nerved, acute, 3 mm long; palea subcoriaceous, linear-lanceolate,
    acute, 2-keeled,as long as spikelet.

    這樣就可以使用下列的正規表示式來找匹配規則:


    [0-9,0-9.]*\ [c|m]m

    其中,第一部分的範圍有[0-9,0-9\.]*代表符合 0-9 的數字或是 0-9 開頭且後面有小數點("."在 vim 中有特別意義,所以保險起見用 \ 跳脫字元),後面接的 metacharacter 則是代表所有 0-9 以及 0-9 開頭後面有小數點的匹配規則都會被找出來。第二部分觀察到的匹配規則則是數字和單位之間有空格,所以我們用" "(空格)接在後面,第三部分則是單位,有 cm 及 mm,寫成範圍的正規表示式則是 [c|m]m (c 或是 m,後面接 m,則為 cm 或 mm)。但這樣寫會有個問題,就是如果文字中數字單位前有加範圍,例如 1-3 cm ,這樣就無法抓到,所以得另外再處理。

  3. 符合匹配規則(pattern)的分組
    有時候我們需要幾個不同匹配規則做不同的取代,那就可以用 grouping 的方式把匹配規則分組,分組是用\(分組的元素\)(一對括號,但前面要使用 backslash "\" 來跳脫字元),而符合前述匹配規則的分組,可以用 \ 加上 1~9 的數字來代表,如果 \0 表示符合所有匹配規則的比如:

    :%s/\(30\) \(cm\)/\1 \2/g
    

    \1 代表了符合 30 的字元

    \2 則代表符合 cm 的字元

所以我們在「Lorem ipsum 福祿壽」這個例子中要怎麼寫 vim 的正規表示式呢?第一,先找出匹配規則,這個匹配規則是可以「精確」描述出來的,這是最重要的步驟;第二,試著用搜尋(指令為 /)來試試看第一步中所找出的匹配規則,如果不對的話,再重複第一步驟修正,直到確定找出為止。

我們希望 Lorem ipsum 是一組詞,用空白分隔,加入 \index 變成\index{Lorem ipsum},而「福祿壽」是第二組詞,加入 \index 變成\index{福祿壽},第一組的 pattern 就是第一個字首字大寫,第二個字的字首小寫。第一組和第二組都被大括號包住,而且大括號前有\name[數字加上"." ],因此用 vim 寫符合上述規則的正規表示式如下:

:\\name\[.*\]{\([A-Z].*\) \([a-z].*\) \(.*\)}

 |--1-|--2--|----3------|4|----5----|6|--7--|

解釋:
第一部分(|---1---|)代表符合有\name的字元,因為"\"已是跳脫字元,所以要多加一個"\"

第二部分則是代表[ ](中括號)包起來的字元

第三部分則是代表Lorem ipsum中的第一個字Lorem,首字大寫,所以使用[A-Z].*,因為我們希望它在 index 時有階層的索引(請參見 latex 詳細說明),所以要特別處理。所以把它用\(\)當成第一個群組

第四、六部分是空格

第五部分則是Lorem ipsum中的第二個字ipsum,因此寫成[a-z].*,是第二個群組

第七部分則是最後「福祿壽」

接下來則是把符合的匹配規則取代成\name[1. ]{Lorem ipsum 福祿壽} \index{Lorem!ipsum} \index{福祿壽},vim 的語法如下:

:%s/\\name\[.*\]{\([A-Z].*\) \([a-z].*\) \(.*\)}/\0 \\index{\1!\2} \\index{\3}/g
   |-----------------------1--------------------|--------------2--------------|

解釋:
取代的部份(|---2---|)中,\0代表符合前述匹配規則的字元,\1 則是「Lorem」,\2則是「ipsum」,而\3代表「福祿壽」。

雖然正規表示式寫起來像火星文,但是只要能夠精確找到匹配規則的話,可以很快速的節省手動取代或修改的時間,多練習就會上手了!有興趣學習 vim 的朋友們可以參考李果正前輩寫的「大家來學 vim」

參考文獻及延伸閱讀:
1. vim regexp
2. 大家來學 vim

PostgreSQL 延伸套件

PostgreSQL 可以透過 pgxn client 來安裝一些第三方的套件,使用起來還滿方便的。

因為 pgxn 是 python 寫的,所以安裝只要透過 pip :

$ pip install pgxnclient

安裝第三方套件:

# 編譯並安裝

$ pgxnclient install _pkgname_
# 將編譯好的套件使用 extension 方式安裝進特定資料庫

$ pgxnclient load -d _database_ _pkgname_

之後就能夠開心的使用第三方套件了

References:
https://github.com/dvarrazzo/pgxnclient
http://pgxnclient.projects.pgfoundry.org/install.html

關於香蕉的空照判釋

[注意!]其實這篇文章不應該是漂流木事件的主軸,但有些細節可以讓大家理解航空照片的判釋是具有高度專業及不確定性,不是隨便腦補就可以做出結論的(茶)


圖零、你能分得出來這些是什麼樹嗎?

關於潘建志醫師臉書及蘋果報導[1]的香蕉空照圖及街景圖來看,會誤以為屏東縣九如鄉榮泉街一號屏26縣道附近的香蕉是呈現圓團狀(參考圖一及圖二)


圖一、屏26縣道九如鄉榮泉街一號附近的農田街景圖(2013)


圖二、屏26縣道九如鄉榮泉街一號附近的農田衛星圖(2014)

根據本鍵盤森林學家的判釋,這塊 2013 年街景圖看起來都是香蕉田,其實是混作(可以仔細觀察一下圖一標註位置)。香蕉通常是每年就會重新種植,而有些農民會間作其他多年生的作物,像是檳榔、柑橘類(金桔、檸檬)等,加上 Google 街景圖和衛星圖的時間差,所以有可能會造成街景圖和衛星圖上的影像對不起來。仔細觀察圖二衛星圖的香蕉混作田,雖然都是圓團狀的柑橘類作物田(只是猜測,我也沒辦法判斷確切的種類),但還是可以看出有留幾棵香蕉在田中(標註處)。因此呼籲鍵盤辦案的鄉民在觀察街景圖和衛星圖時,除了注意拍攝時間外,畢竟可公開取得的衛星或航照影像圖解像力有限,還是要靠現場或其他證據輔助判斷。

另外對於航照判釋是門專業的技術,想起當年修測量學提到航空測量判釋時,都需要使用立體鏡觀察立體相片對,因為只有透過立體觀察可以看出景深並判斷物體高度(有時候影子是輔助判斷物體時很好的特徵),平面的影像比較不容易看出三維。至於是否能判斷樹種就需要經驗和野外確認了,像是高海拔樹種單純、特徵明顯,所以可以從航測照片看出台灣冷杉(Abies kawakamii)、台灣鐵杉(Tsuga chinensis var. formosana)等,但低海拔闊葉樹樹冠特徵相近,比較不容易確定,此時只能依照植群型和地點來判斷大概是哪些物種,再透過實地訪查確認。

最後整理一下低海拔從空照圖看起來為傘狀或放射狀的樹種,大部分是棕櫚科(Arecaceae)植物,也有些是芭蕉科(Musaceae)或樹蕨等:

A. 天然林

  1. 黃藤類(C. quiquesetinervius,包括水藤 Calamus formosanus 這類攀援性木本植物),看起來會在一些森林邊緣或樹冠層上,如下圖三,這類的森林多半都是干擾頻度較高、攀藤類的植物也較多,航照圖上看起來 pattern 不明顯(sorry 懶得找航照圖了)


    圖三、黃藤類植物

  2. 山棕(Arenga tremula)
    後來真的有學生跑去實地訪查[3],拍了一些照片,潘醫師 blog 內(提到的木材)的應該是山棕,跟香蕉比起來,山棕可以超過三米,葉幅也較香蕉長,加上山棕為一回羽狀複葉從基部生長,小葉會下垂,所以在較低解析度的空照俯瞰圖上看起來會較瘦長,所以潘醫生才會誤認為是木材。另外山棕通常分布在溪溝谷地較溼的環境,而像這樣的生態環境也可以做為輔助判斷

  3. 台灣芭蕉(Musa basjoo var. formosana)
    這個也可以把香蕉的部分拿來一起談,兩個外型很接近的。台灣芭蕉多半零星分布在林下,不容易直接從空照圖看到,下圖是天然闊葉林內的台灣芭蕉


    圖四、這張圖剛好有台灣芭蕉、黃藤和山棕 XD

  4. 樹蕨類
    台灣因為潮溼,所以樹蕨類(多半指桫欏科)也很常見,尤其是低海拔開闊地、道路旁或是崩塌的溪溝谷地常常可以看到筆筒樹(Cyathea lepifera),其特徵是葉子為三回羽狀複葉,其葉寬也比芭蕉或是山棕來得寬,加上筆筒樹葉子著生的位置接近輪生,所以航照圖看起來像是一朵朵小花,跟香蕉類或棕櫚科植物比起來葉部份的陰影也比較少,形狀也比較圓潤一些,請參照圖五及圖六


    圖五、陽明山一帶的筆筒樹

    圖六、如果你仔細尋找道路旁或山澗溪溝旁的空照圖,也可以找到筆筒樹。紅色圓圈標記處下方有一條黑色的陰影,這個可能是路徑或溪溝。

B. 作物或栽培植物

椰子類(可可椰子、亞歷山大椰子等)或檳榔,因為屏東種植椰子和檳榔很多,所以我就拿這兩種主要的栽培作物來做比較。


圖七、可可椰子和檳榔比較,能看出傘狀,樹冠幅度較大的為椰子。

本鍵盤森林學家老家就是種可可椰子和檳榔的,所以我對這兩種作物都十分熟悉。可可椰子高度可以達到五、六層樓(~15m)、冠幅約 5–7 公尺,所以從空照圖來看陰影會很深,如果使用立體相片對觀察,圖七椰子和檳榔交界處就可觀察的出高低的層次差異,而檳榔因為樹冠幅度窄,種植間距也比較近,這一點也是可以做判斷的依據。

結論:航照及衛星影像圖只能當成輔助判斷的依據,並不能直接腦補當成證據的,要知道真相還是要實地訪查和拍照 XD

參考文獻及引用地圖

[1] 關鍵木材變「香蕉樹」?霹靂潘挨轟邱毅2.0
[2] 高品質香蕉田間栽培管理手冊,台灣香蕉研究所
[3] PTT: [爆卦] 內湖山老鼠Google拍到了 興善宮也是違建
[4] Google Maps 圖五筆筒樹之原圖
[5] Google Maps 圖六筆筒樹之原圖
[6] Google Maps 圖七可可椰子與檳榔比較之原圖

zfs 單一裝置儲存池改為鏡映

標題有點難下 XD

前言及環境:

假設已經有一台 FreeBSD 預設將 root 安裝至 zfs 檔案系統,但一開始並沒有設定 mirror 或是 raidz ,之後想要將系統換硬碟,將原本的舊硬碟資料轉移出來,並設定 mirror (換句話說就是:一開始裝系統時沒錢買多餘硬碟,只好將就著用,等有經費時,可以買兩顆硬碟和新主機時,要怎麼比較無痛把舊機器舊硬碟上的資料整個放在新的系統上,並設定比較安全的儲存環境)因為網路上有許多是一開始就設定好 zfs mirror 或 raidz 取代壞掉的硬碟教學,所以本篇主要是著重在怎麼從無冗餘硬碟的狀況下,再增加 zfs mirror 裝置。

系統: FreeBSD 10.1
舊主機之硬碟: 代號假設為 A
新主機的兩顆硬碟: M0, M1

方法:

如果照原本 UFS / MBR 的方式,就是先將新的硬碟裝上去,分割好並格式化,掛載成 /mnt/root,用 dd 或 dump | restore 的方式將資料倒回去。等資料都倒完後再將舊硬碟取下,並設定成新硬碟開機,步驟不難但是很麻煩,所以這篇教學試著使用 zfs zpool resilvering(「銀鹽刻版複製」,方式類似 RAID 的重建[rebuild])的概念來置換舊硬碟資料,用新硬碟開機,並將原本無 mirror 的 zpool 增加磁碟鏡映群組,新增硬碟並重建資料。做法的概念圖如下:

實作:

  1. 先將新的硬碟裝上去,並用 gpart 分割磁區(partition),設定好 FreeBSD 開機磁區、SWAP等
    假設這顆新的硬碟 M1 的裝置節點是 /dev/ada1 ,以下指令請記得都用 sudo 或 root 身分執行
    # 建立 GPT 磁區
    
    gpart create -s gpt ada1
    # 新增 freebsd 開機磁區為,大小為 512 KB,標籤(-l)為 gptboot
    
    gpart add -t freebsd-boot -l gptboot  -s 512K ada1
    # 新增 swap 磁區,大小為 8192 MB
    
    gpart add -t freebsd-swap -l gptswap -s 8192M ada1
    # 剩下的全部分割給 freebsd zfs 磁區,標籤為 zfs1
    
    gpart add -t freebsd-zfs -l zfs1 ada1
    # 設定下次開機的硬碟為本顆硬碟的第一個分割區(ada1,參數 i)
    
    gpart bootcode -b /boot/pmbr -p /boot/gptzfsboot -i 1 ada1
    
    此時,可以到 /dev/gpt 目錄下確認是否有 zfs1 ,若有的話就可以繼續次步驟
  2. 將新硬碟 M1 的 zfs 磁區加入 zpool,並設定成 mirror

    zpool attach zroot gpt/zfs0 gpt/zfs1
    

    使用 zpool status 來確認,如下圖,此時 gpt/zfs1 顯示 resilvering,可能需要一些時間處理(端看資料量大小)

  3. 將舊硬碟 (A) 離線
    這個步驟將舊硬碟先離線,等待關機之後(不過理論上 SATA 硬碟可以熱拔插)換上另一顆新硬碟:

    zpool offline zroot gpt/zfs0
    

    zpool status 確認時,原本 gpt/zfs0 標籤會變成一串數字,STATE 則會顯示 OFFLINE

  4. 移除舊硬碟,安裝另外一顆新硬碟(M2),並重複步驟 1 的 gpart

    換硬碟時,記得不要拿錯 XD 下面的步驟中,在重開機後,原本的硬碟 M1 變成 ada0,而新的 M2 則是 ada1,所以在分割磁區設定標籤時,將 freebsd-boot 改成 gptboot0,而 M2 上的新 zfs 磁區標籤則為 zfs0

    # 建立 GPT 磁區
    
    gpart create -s gpt ada1
    # 新增 freebsd 開機磁區為,大小為 512 KB,標籤(-l)為 gptboot
    
    gpart add -t freebsd-boot -l gptboot0  -s 512K ada1
    # 新增 swap 磁區,大小為 8192 MB
    
    gpart add -t freebsd-swap -l gptswap -s 8192M ada1
    # 剩下的全部分割給 freebsd zfs 磁區,標籤為 zfs0
    
    gpart add -t freebsd-zfs -l zfs0 ada1
    
  5. 使用 replace 更新磁碟資料

    離線的磁碟在 zpool 中會顯示一串數字,更新磁碟時記得複製這串數字,如下圖

    zpool replace zroot 9450223813957032699 gpt/zfs0
    


    接下來就是等這顆新磁碟自動 resilvering 資料即可
    psilotum@biodiv ~ $ zpool status
    pool: zroot
    state: ONLINE
    scan: resilvered 98.4G in 0h20m with 0 errors on Thu Feb 5 19:01:57 2015
    config:
    NAME          STATE     READ WRITE CKSUM
    zroot         ONLINE       0     0     0
      mirror-0    ONLINE       0     0     0
        gpt/zfs0  ONLINE       0     0     0
        gpt/zfs1  ONLINE       0     0     0
    

    errors: No known data errors


    參考資料

    https://www.freebsd.org/doc/handbook/zfs-zpool.html
    http://www.freebsdwiki.net/index.php/ZFS,_booting_from

替換 windows ^M,並練習組合一些簡單 unix 指令

有時候拿到一些在 Microsoft Windows 處理的文字檔案,會有一些斷行後有 ^M (Ctrl+M) 的字元,在 unix 上處理這類的資料很惱人,會造成不正確斷行判斷,所以可以用 find, xargs, awk, sort 和 vim 一次取代完(^M 要按 Ctrl 加上 M)。分解做法如下:

  1. 找某目錄底下以 .txt 為名的所有檔案 實作:
    $ find . -type f -name "*.txt"
    
    解釋:
    .(dot) 代表現在目錄底下
    -name 後面接的參數代表檔名,可以用 wildcards 或 regular expression
    -type 後面接的是檔案類型, f 代表一般的檔案,d 代表目錄等
  2. 找完後用 xargs (eXecute ARGuments)

    xargs Manual 這樣寫:

    The xargs utility reads space, tab, newline and end-of-file delimited strings
    from the standard input and executes utility with the
    strings as arguments.

    簡單說就是把上個程序輸出的結果當成參數,餵給後面的程式使用,有點難懂對吧?那他跟 pipe (|) 有什麼差別呢?
    pipe 字面上是水管,只負責「傳導」這件工作,像是:

    $ ls | grep keywords
    

    這個指令就代表把 ls 的輸出導給 grep 執行,而 grep 把輸出中有 keywords 的行抓出來。

    若是下列指令則會用 grepls 所列出所有的檔案抓取有 keywords 的行

    $ ls | xargs grep
    

    有了 xargs 的基本概念後,用例子來解釋,下面用 ls 會列出目錄底下有三個 .txt 檔案,分別為 1.txt, 2.txt, 3.txt。「列出這三個檔案」就是標準的輸出

    $ ls
    1.txt 2.txt 3.txt
    

    接下來若導向給 xargs,則會把標準的輸出當成是參數清單(argument list),
    ls | xargs grep keywords 執行所得到的結果則和下列指令相同:

    grep keywords 1.txt
    grep keywords 2.txt
    grep keywords 3.txt
    
  3. awk 取出檔案名稱
    因為 xargs grep keywords 的輸出類似:

    1.txt: ipsum lorem keywords consectetur adipiscing elit, sed do eiusmod
    1.txt: lorem ipsum quasi voluptate velit esse keywords dolore eu fugiat
    2.txt: voluptate velit esse cillum dolore eu fugiat keywords, lorem

    因此我們希望只列出檔案名稱,觀察上述的 pattern,就是取出 : 前的文字,使用指令為:

    awk -F':' '{ print $1}'
    

    -F 後面用 '(single quote 單引號) 包起來的是分隔符號(separate Fields),後面加數字代表分隔符號切分後第幾個欄位,例如上面的 3. 的例子中, '{print $1}' 則會將 1.txt, 2.txt 等顯示出來

  4. sortunique 整理

    前面 awk 處理完後把含有 ^M 的檔案名稱丟給 sort -u (排序整理,並刪除重複)

  5. 把整個程序寫成迴圈,用 vim 取代 ^M
    shell script 迴圈:

    for _var_ in _condition_ 
    do
        _statements_
    done
    

    shell script 中變數是加上 符號,將步驟一到四所找到有 ^M 的檔案名稱當成輸入的清單,跑個迴圈再加上 vim 取代即可。在 shell 中 vim 後面加上 + 則會執行 vim 內部的指令,我使用 regular expression 取代 ^M,指令為:

    %s/被取代的文字/取代的文字/g
    

    記得最後還要加上寫入離開指令 +wq

最後把程式整理一下,完整程序如下:

#!/usr/bin/env bash

for i in `find . -type f -name "*.txt" | xargs grep "^M"| awk -F':' '{print $1}' | sort -u`
do 
    echo "Substitute ^M in ${i}"
    vim +%s/^M//g +wq ${i}
done

參考資料

http://en.wikipedia.org/wiki/Xargs

升級 PostgreSQL 主要版本

PostgreSQL 的版本從 8.3 之後,就是以 X.Y 當成是主要版本(major release),而每個版本都會有小修正的次要版本 X.Y.Z (minor release)。若要升級主要版本,例如從 9.3 升級至 9.4 ,因為主要版本升級會加入一些新功能,而系統資料表的配置(layout of system tables)則會改變,所以需要使用 dump / restore 來重新構建系統資料表(是的,很麻煩),在 PostgreSQL 升級主要版本前,可以用 pg_upgrade 來升級。升級的流程如下,以 MacOS X 底下 homebrew 安裝為例:

  1. 停止 postgres daemon 傳統用法是透過 service 來控制 /etc/rc.d/ (或 /usr/local/etc/rc.d) 底下的 scripts 來控制 daemons,例如
    $ service stop postgres
    
    但 MacOS X 則是使用 /Library/LaunchAgents/Library/LaunchDaemons 底下的 plist 來控制,所以先將其 unload (同時也會呼叫 pg_ctl 去停止 postgres 運作):
    $ launchctl unload ~/Library/LaunchAgents/homebrew.mxcl.postgresql.plist
    
  2. 將舊資料庫儲存位置更名 homebrew 預設會將 PostgreSQL 的資料庫儲存於 /usr/local/var/postgres 底下,因為我是從 9.3 升級至 9.4,所以複製為 postgres93 (whatever)
    $ mv /usr/local/var/postgres /usr/local/var/postgres93 
    
  3. 升級現有版本之 postgres 及相關函式庫 如果你有使用相關延伸函式庫(例如 postgis),請記得也一起升級
    $ brew upgrade postgresql postgis
    
  4. 備份 雖說可以用 pg_upgrade 更新,但還是把所有資料庫備份出來比較保險
    $ pg_dumpall > postgres93.dbout
    
  5. 建立新的主要版本資料儲存庫 新建儲存庫
    $ mkdir /usr/local/var/postgres
    $ initdb -D /usr/local/var/postgres
    
    使用 launchctl 載入 postgresql (記得不要在 tmux or screen 內執行,否則會有 permission 的問題)
    $ launchctl load ~/Library/LaunchAgents/homebrew.mxcl.postgresql.plist 
    
  6. 使用 pg_upgrade 升級主要版本 因為需要用到舊版執行檔,所以記得不要太快移除舊版本(這也就是使用 homebrew 的好處,每個第三方套件都會裝在 /usr/local/Cellar 底下,再連結到 /usr/local 中相關的目錄位置) 四個參數如下: --old-datadir: 舊版本的資料儲存目錄, --new-datadir: 新版本的資料儲存目錄, --old-bindir: 舊版本的執行檔目錄, --new-bindir: 新版本的執行檔目錄。 確認上述的目錄和版本後,接下來就可以使用 pg_upgrade 來升級
    pg_upgrade
    --old-datadir /usr/local/var/postgres93
    --new-datadir /usr/local/var/postgres 
    --old-bindir /usr/local/Cellar/postgresql/9.3.5_1/bin 
    --new-bindir /usr/local/Cellar/postgresql/9.4.0/bin
    
    升級過程會確認版本,然後 dump 舊有資料庫,並轉成新版本:
    Performing Consistency Checks
    -----------------------------
    Checking cluster versions                                   ok
    Checking database user is a superuser                       ok
    Checking for prepared transactions                          ok
    Checking for reg* system OID user data types                ok
    Checking for contrib/isn with bigint-passing mismatch       ok
    Checking for invalid "line" user columns                    ok
    Creating dump of global objects                             ok
    Creating dump of database schemas
                                                            ok
    Checking for presence of required libraries                 ok
    Checking database user is a superuser                       ok
    Checking for prepared transactions                          ok
    ...(略)
    Optimizer statistics are not transferred by pg_upgrade so,
    once you start the new server, consider running:
    analyze_new_cluster.sh
    ...(略)
    Running this script will delete the old cluster's data files:
    delete_old_cluster.sh
    
  7. 新版本資料庫最佳化及刪除舊版本資料 升級完畢之後,會在你執行 pg_upgrade 的目錄底下自動產生兩個 shell script 檔案, 執行 analyze_new_cluster.sh 可以最佳化資料庫($PATH_TO_POSTGRESQL/bin/vacuumdb" --all --analyze-in-stages),執行 delete_old_cluster.sh 則會把舊版本的目錄刪除(其實就是`rm -rf postgres.old)
    ./analyze_new_cluster.sh && delete_old_cluster.sh
    
  8. 升級 Postgis extension 若有安裝 postgis extension,所以就順便升級
    spatial_db=# ALTER EXTENSION postgis UPDATE TO "2.1.5";
    spatial_db=# ALTER EXTENSION postgis_topology UPDATE TO "2.1.5";
    

References:

pg_upgrade

ZFS 檔案系統筆記(1)

#zfs #filesystems #freebsd

ZFS 是一個注重資料完整性(data integrity)的檔案系統,也支援磁碟陣列、檔案系統快照(snapshot)、支援大容量磁碟等特色。ZFS 最早是在昇陽(Sun Inc.) Solaris 作業系統上所開發,但昇陽被甲骨文(Oracle Inc.)併購後,轉成開源的 Open Solaris 作業系統[1]。而眾多 Unix-like 的作業系統也支援 ZFS,例如 Freebsd, Linux, MacOS 等,但社群 在 2013 年 9 月也將 ZFS 移植並創建 OpenZFS 專案繼續維護 ZFS。FreeBSD 在 9.

ZFS 主要的設計目標

  1. 資料完整性,所有的資料都會用雜湊演算法計算出確認索引值(checksum),因此只要檔案有變動或出現錯誤,確認索引值即會不同,這對於裝置上的資料完整性是很好的保障。
  2. 資料儲池(pool),(這就是大水庫理論啦(誤XD)) 傳統檔案系統都是以資料儲存裝置為一個基本單位,啟始資料儲存裝置需要新建立分割區、格式化等工作,在一般的 UNIX 上面需要寫入檔案系統表(fstab),一但需要增加空間或修改則十分麻煩。所以 ZFS 設計上以資料儲池為一個單位,若要增加儲存裝置為鏡映、RAID或是修改都較容易
  3. 效能,支援先進的儲存裝置快取功能,包括 adaptive replacement cache (ARC)、L2ARC、ZIL

我為什麼要使用 ZFS?

  1. 因為宅宅都用跟人家不一樣的檔案系統,既然都用 FreeBSD 了, 內建 ZFS 好物不用嗎(而且現在 10.1-RELEASE 連 root filesystem 都可以直接支援 ZFS 開機了
  2. 和傳統 unix 檔案系統比較,以 UFS (unix file system) 為例,在檔案發生問題(file corruption)時, 使用磁碟掃描(fsck)修復時需要將磁碟卸載為離線狀態,才能修正,而且一般的 fsck 只能針對檔案系統的 metadata 修正,若檔案真的爛掉,是無法救回的。像是我們每年都會停電檢修(忘記關機),要不然偶爾就來個跳電,雖然 FreeBSD 有 background fsck 功能,但就有時還是得進入 single user mode 開個 fsck 修一下
  3. 對於沒錢買許多硬碟的人來說,一但硬碟壞掉很麻煩,即使有 RAID ,rebuild 也很麻煩
  4. 支援檔案系統快照,這個功能非常好用,如果你有用過虛擬機器(virtual machine),這個功能在各大虛擬機器軟體(virtualbox, vmware, etc.)都有類似的功能
  5. 備份很方便啊,建立檔案快照,儲存在另一個裝置或是壓縮丟遠端都可以
  6. 支援 jail!

名詞解釋

dataset(數據集): zfs 檔案系統內元件的通稱,包含複本(clone)、檔案系統、檔案系統快照等
pool(儲池): 資料置放區稱之為資料儲池(zfs storage pool)
mirror(鏡映): 同時儲存資料於兩個及以上之裝置,相當於傳統磁碟陣列中的 RAID1
resilvering(銀鹽複製): 將資料從一個裝置複製到另一個裝置,在 zfs 中稱"resilvering"
scrub(刷洗): 相對於傳統 unix 檔案系統的 fsck,ZFS 稱之為 scrub,若使用 scrub,ZFS 會去檢查儲存資料的確認索引值並修正有誤的

安裝

目前 FreeBSD 10.1-RELEASE 已支援 zfs 於 root filesystem,在預設安裝系統中可選擇 guided zfs,按照指示進行即可

常用指令

1. 建立資料儲存池(zpool create pool)

用法:zpool create poolname DeviceNode

root@m15 ~ # zpool create wd1t da1

在 /dev/da5 上建立名為 jetflash64 的儲存池,並掛載於 /wd1t (-m mount)

root@m15 ~ # zpool create -m /jetflash64 jetflash64 da5

顯示現在有的 zpool

root@m15 ~ # zpool list
NAME     SIZE  ALLOC   FREE   FRAG  EXPANDSZ    CAP  DEDUP  HEALTH  ALTROOT
jetflash64  58.5G    79K  58.5G     0%         -     0%  1.00x  ONLINE  -
zroot        928G  4.83G   923G     0%         -     0%  1.00x  ONLINE  -

2. 刪除資料儲存池

zpool destroy _poolname_

3. 檔案系統快照

ZFS 一個特點就是具有檔案系統快照的功能,可以在任意時間點建立檔案系統快照

zfs snapshot zpoolname@snapshot

列出現有的檔案系統快照

zfs list -t snapshot

回溯檔案系統快照至 20141226

zfs rollback -r zpoolname/home@20141226

ZFS 說明非常詳盡,有興趣的人可以讀一下,之後我把 manual 裡頭的一些範例翻譯出來再給各位參考(續)

參考資料
[1] Open Solaris https://solaris.java.net
[2] OpenZFS http://open-zfs.org/
[3] Wikipedia 關於 ZFS 說明 http://en.wikipedia.org/wiki/ZFS
[4] FreeBSD 使用者手冊中的 ZFS 說明 https://www.freebsd.org/doc/handbook/zfs.html
[5] ZFS tutorial http://thegeekdiary.com/zfs-tutorials-creating-zfs-snapshot-and-clones/
[6] http://article.denniswave.com/9218
[7] http://article.denniswave.com/9221

使用 R 處理 LiDAR xyz 資料產生表面數值地形圖

將 LiDAR 資料轉成 STM 網格 (1)

概念及做法

LiDAR 掃出來的資料是「點」資料,若要靠這些反射回來的點座標推得一個區域內的數值地形圖(Digital Terrain Map; DTM,雖然 LiDAR 掃出來的應該是「表面」Surface Terrain Model,若要取得真正的數值地形圖則必須要再做轉換),直覺上最快的作法就是將每個相同座標的點找最小值。假設我們透過 LiDAR 掃描取得的 xyz 座標為(單位為公尺):

id x y z
1 -0.1061 -0.1624 1178.482
2 -0.1061 -0.1623 1178.475
3 -0.1060 -0.1623 1178.474
4 -0.1060 -0.1624 1178.481
5 -0.1060 -0.1623 1178.473
6 -0.1059 -0.1623 1178.472

示意圖如下,假設解析度是 1 cm,LiDAR 掃完某區域後,得到七個點,這些點分別投影在平面後,屬於三個網格,接下來只要找每個網格單位內,LiDAR 掃到的點之最低 z 值,就是表面的高度了。將這些 z 值和網格座標彙整,利用 GIS 即可轉換成 raster 格式的 STM。

實際上,我們可以只針對 x, y 值去換算成網格的資料,若解析度為 1cm ,1cm 以下的小數位數則捨去不計。例如某一點 xy 座標為 (1.51, -0.53),那轉換成 1cm 解析度網格的座標則是 (1, 0)。轉換成新的網格座標系統後,找出相同網格內的所有 z 值,求最小值則是該網格的地表高度。

在我們上面的範例資料中,若用 0.1cm 解析度的網格計算,則可得下列表:

id x y z
1 -0.106 -0.162 1178.482
2 -0.106 -0.162 1178.475
3 -0.106 -0.162 1178.474
4 -0.106 -0.162 1178.481
5 -0.106 -0.162 1178.473
6 -0.105 -0.162 1178.472

計算相同網格座標內的最小值,則會得到:

id x y z
5 -0.106 -0.162 1178.473
6 -0.105 -0.162 1178.472

R 實作

0. 安裝套件

必須安裝的套件(package): raster
選擇性安裝的套件:rhdfs

1. 將資料匯入

可以使用 read.csv()。因為載入資料數量很大,下面例子採用 RHadoop 為例,你可以依照需求修改:

Sys.setenv(HADOOP_CMD="/usr/local/bin/hadoop")

library(rhdfs)
library(rmr2)
hdfs.init();

# import data into hdfs via hdfs.file()
HDFS_r = hdfs.file("ST000.xyz", "r",buffersize = 104857600)
HDFS_f = hdfs.read(HDFS_F)
data = read.csv(textConnection(rawToChar(HDFS_f)), sep=" ", header=F)

如果你不知道什麼是 Hadoop ,請跳過,可直接使用 read.csv() 來匯入資料

data = read.csv('~/Documents/Dropbox/projects/LiDAR/lidar_process/sample_data.xyz',
                header=F, sep=" ")
# use head to view top six rows
head(data)

上面的資料中,V1V4 欄位分別代表相對座標 (x, y, z) (相對於 LiDAR 儀器擺放的中心),以及反射強度值。

2. 將 (x, y) 座標取出,用無條件捨去法處理

接下來將 dataset 中的 x, y 座標取出來,並假設解析度為 10 cm。因為原始資料單位是公尺,為了方便計算,將其換算為公分後,再無條件捨去,在這個步驟中,我們會使用到 trunc 這個函式:

trunc_xy = trunc(data[,1:2]*100/10)*10
head(trunc_xy)

處理完之後,把四捨五入的 x 及 y 座標用 column bind (cbind) 合併在一起。

trunc_xyz = cbind(trunc_xy, data[,3])

# 為了讓 dataframe 美觀且容易閱讀,我們加上欄位名稱
colnames(trunc_xyz) = c("x", "y", "z")
head(trunc_xyz)
3. 使用 aggregate 求出相同網格(xy 座標)中的最小值

這個步驟主要是把 (x, y) 合併當成一個群組,相當於 group by 的概念,在一個 group 中求出 z 的最小值。

# aggregate(target, by=_grouping_elements, FUN=_function)
aggr_xyz = stats::aggregate(trunc_xyz["z"], by=trunc_xyz[c("x", "y")], FUN=min)
aggr_xyz

所以上面的 dataframe 就是我們所要求的 10cm 解析度的網格 xyz 座標值。

4. 寫成函式

整理一下上述步驟,將處理的程序寫成函式,方便後續的處理:

輸入之變數: 資料集(dataframe)、網格解析度(res)
使用的功能函式:stats::aggregate(), base::trunc()

proc_stm = function(dataframe,res) {
  r_cm = 100
  trunc_xyz = cbind(trunc(dataframe[,1:2]*r_cm/res)*res, dataframe[,3])
  colnames(trunc_xyz) = c("x", "y", "z")
  output = aggregate(trunc_xyz["z"], by=trunc_xyz[c("x","y")], FUN=min)
  return(output)
}

處理 STM:

output.xyz = proc_stm(data, 10)
output.xyz
5. 將 xyz 座標轉換成 raster 網格

因為 LiDAR 具有距離限制,我們希望以 LiDAR 架設點為原點, 30 公尺的正方形限定範圍:

# limit the target area to 30 by 30 meters
bound=15
output_cm.xyz = output.xyz[(abs(output.xyz[,1]/100) < bound & abs(output.xyz[,2]/100) < bound),]
#換算回公尺
output_m.xyz = cbind(output_cm.xyz[,1:2]/100, output_cm.xyz[,3])
output_m.xyz

接下來使用 rasterFromXYZ() 將 xyz 座標轉成 raster 格式,並使用 writeRaster() 輸出為 GeoTiff raster 格式:

# use raster library to output raster file
library(raster)
rast = rasterFromXYZ(output_m.xyz)
writeRaster(rast, filename = 'output.tif', overwrite=TRUE)

下圖則是輸出的範例(非本例中之資料輸出)

如何有效率的「問」問題

在很久很久以前,Eric Raymond(@esr, 他就是寫了《教堂與市集》的資深駭客,堪稱為自由軟體界的歷史學家)寫了一篇《How to ask question the smart way》,直到今年五月份 esr 仍然還有修訂這份經典之作。如果你沒有讀過這份文件,我強烈建議你讀一次,雖然是英文的,但文件本身非常輕鬆好讀。

簡單講一下一些大原則,在你提出問題時,可以問看看你自己。

搜尋

  • 你真的有好好用過 Google / DuckDuckGO 等搜尋引擎找關鍵字嗎?(關鍵字:stackexchange, stackoverflow)
  • 你真的搜尋不到你要的相關主題或提示嗎?
  • 拼字有拼對嗎?

搜尋主要是要有針對性與特定的問題,才能比較好過濾資訊,這個每個人都會有一些小撇步。依照我的經驗法則,通常可以用「方法」加上「使用軟體/步驟」及「文件類型」等來搭配,常常在搜尋引擎前三到四個頁面就可以找到我想要的答案(但當你什麼都不懂時,請先去學校上基礎課程,或是找一些 open course 來自修)。例如:

我想要找空間資訊中的交集問題,你可以用

gis + intersection

但這樣的結果往往是出現最多人使用/或最常被搜索的關鍵字,所以若你已經知道有哪些軟體或方式可以做,你可以嘗試修改一下關鍵字:

postgis + intersection

但這樣的搜尋結果也許不盡滿意,如果只是想看一下 postgis 做交集的範例說明,可以用下列的關鍵字搜尋:

postgis + intersection + tutorial

若你不是想找文件和教學,而是實際上遇到了問題卡住,例如你今天做交集時,遇到 ERROR: function st_intersect(geometry, raster) does not exist,也許你可以用:

postgis + st_intersect + error 

來找尋相關的解法(上述例子是 st_intersect() 拼錯字,字尾少了一個 s,這個有可能會有搜尋引擎自動幫你更正)

提出問題

我舉幾個最近有人問我的問題來回答(為了相關人的隱私,所以只節錄問題主幹,問候語與署名等皆去除),這樣比較容易讓大家了解(別人問的問題,也許就是你想問的問題)為什麼這樣問,別人無法

範例一

... 我想請問你
未來氣候資料的兩個子情境因子
A,B兩個分別是代表什麼啊??

解析:

這個例子中,很明顯就是事主連搜尋都不願意花功夫,即使用「未來氣候+情境+A+B」這樣的中文關鍵字也能找到答案。
如果你問我這個問題,我回答的意願趨近於零。

範例二

... 今天老師覺得上次請你幫忙的DEM及坡度好像又有問題
你上次給我的檔,之後老師有再轉成網格檔,我再切成各行政區,然後再用ZONAL STAT算出每個地籍的平均坡度。
也許是我轉檔的過程中又有問題,但一時間也想不出來 ...

解析:

這個「問題」的問題出在語焉不詳,雖然已經有提出使用的方法,但是過程仍然沒有交代清楚。以下舉幾個常犯的錯誤(我以前也是有犯過類似的錯誤):

這個看起來怪怪的(哪裡怪?有像怪老子一樣怪嗎?)
看起來有問題(怎麼看?閉著一隻眼睛看?還是倒立看?)
算出平均坡度後有問題(你怎麼算?)

事主一直遇到問題,一直問了「問題」,但是沒有詳細的步驟,也是個很難有人想回答的問題。

如果你這樣問問題,我會樂意幫你找出解答:

我用 gdaldem 把你上次給我的 dem.tiff 檔案轉成坡度檔,轉檔的方式為:

    gdaldem slope dem.tiff slope.tiff

但我用 ArcMap 匯入 slope.tiff 和之前做的 slope2013.tiff 比較後發現數值差很大,附檔是我用 arctoolbox 中的zonal statistics 計算 land.shp 和 slope.tiff 的輸出值.... (後略)

如果是這樣問問題,我就能馬上告訴事主,他使用的 gdaldem 少加了一個參數 -p (gdaldem 預設是計算度,要加 p 才會計算 percent),只需要我看完問題就能馬上告訴事主問題點在哪。

範例三

... 匯入X.Y點位後要輸出成shapefile檔,可是檢查資料夾時卻只看到look檔,是哪裡出問題了呢?

解析:

我會問你的是:

你用什麼軟體做?ArcGIS/QGIS/R? 什麼是 look 檔?是附檔名 .look (當然我猜有可能是 .lock 檔)?還是資料夾中只有 look 這個檔案或目錄?
然後我就覺得很麻煩,你問我問題我還需要擲杯問一下媽祖還是土地公,回答的意願也是趨近於零。

若你這樣問問題:

題目:我用 ArcMap 10.0 版將 xy 座標匯入後,轉成 shape 檔失敗,並只出現附檔名為 .look 的檔案,請問如何解決呢?

使用的軟體及輸入格式:
ArcMap 10.0
xy 座標檔案格式: csv

步驟:
將 csv 檔案匯入 ArcMap 中後,按右鍵 display xy,接下來載入該圖層後,另存新檔至 C:/tmp/ 則出現失敗字樣
error 代碼及訊息如下
...(略)...
重複數次都是得到相同的結果,請問這個問題大概出在哪邊呢?

上面的方式有附上錯誤訊息及相關的操作步驟,這個會更容易找出錯誤在哪。有時候我光是把我的作法寫下來,看一下錯誤訊息,透過搜尋相關的問答,換個方法或調整錯誤的參數,試一下大概就能解決了。

最後請記得,沒有人有義務回答你的任何問題,大家都是花了自己額外的時間和精神互相幫忙,在別人幫你之前,請你自己先幫自己,作業要自己做,不要抄別人的。