A Big Fish in a Little Pond (もきログ): 2013

2013年7月25日木曜日

シカ再び

随分時間がかかったけど、ようやくシカ論文の修正原稿を書き上げた。
共著者に確認してもらっているので、8月中には投稿したい。

...が、ひとつ大きな問題が....
お金がないので英文校閲が受けられないかも....

2013年6月30日日曜日

メッシュコードから地図（GeoTiff画像）を作成するプログラム（Mesh2Map）を作成しました。

日本の統計データは大抵メッシュコード（１次、2次、3次などなど）で管理・公開されています。
メッシュコードから直接地図を作成できると便利だなと思って（2年くらい前に）作りました。

作ったと言っても、実態はGDALツールのgdal_rasterizeのラッパーなので、大したことはありません。
単に、メッシュコードから緯度経度座標を計算して、それを元にベクター・ラスター変換（gdal_rasterizeを実行）しているだけです。

手前味噌ですが、ちょっと作ってみたデータを簡単に表示させたい場合なんかにとても重宝します。
あくまでもデータ確認用なので、CRSはWGS84に決め打ちしています。

GUIや着色機能も追加しようと思ったのですが、あんまり必要性を感じられなかったので、そのまま公開することにしました（気が向いたらGUI作ります）。

特に問題ないと思いますが、使用の際は自己責任でお願いします。

使い方

1. プログラムのダウンロード&好きな場所に解凍

ダウンロードはこちら

ちなみにWindowsでの動きます。

2. データの準備

テキストエディタで、以下のような書式でデータを作成＆保存。
エクセルで、[csv形式で保存]でも良い。

データは先ほど解凍したフォルダ（プログラムがある場所）に保存する。

Data.txt
**************************************************************
MeshCode Data1 Data2 Data3
54391210 1 2 4
54391211 2   3    5
54391212 3     8     2
54391213 4     5     1
54391214 5    1   3
. . . .
. . . .
. . . .
**************************************************************
※1行目には列名を入力する
※列名はアルファベット推奨, 日本語可, 数字のみ不可
※2行目以降はデータ
※区切り文字は空白かまたはカンマ、タブ
※1列目はメッシュコード
※メッシュコードは1次, 2次, 3次, 5km, 100mメッシュに対応
※2列目以降は属性データ
※拡張子は何でも良い（.csv, .txt, .dat）

2. 解凍したフォルダ内の[usage.bat]を右クリックして[編集]から、以下の書式で実行コマンドを入力
ちなみに、実行したくない行は文頭に"::"を入力（usage.batを参照の下さい）。

sample_sjis_win.exe <入力ファイル名> <出力csvファイル名> <地図化する列名>

例

Mesh2Map_sjis_win.exe Climate.txt Result.csv WI

3. usage.batをダブルクリックするとプログラムが実行される。
上の例では、Climate.txtというデータを入力すると、WIという列名のデータを使って、WI.tifという地図画像（GeoTiffファイル）が出力される。
出力した画像はQGISで確認できる。

また、Result.csvに緯度経度座標がくっついたデータが出力される。
このファイルも、QGISのデリミテッドテキストで開くとそのまま使うことができる。

Result.csv
**************************************************************
Lon,Lat,MeshCode,Data1,Data2,Data3
139.25625,36.09583333,54391210,1,2,4
139.26875,36.09583333,54391211,2,3,5
139.28125,36.09583333,54391212,3,8,2
139.29375,36.09583333,54391213,4,5,1
139.30625,36.09583333,54391214,5,1,3
. . . .
. . . .
. . . .
**************************************************************

注意事項

※色づけ, 投影変換などなどGdal toolで可能な動作は実装可能
※CRSは緯度経度座標系のWGS84 です。
※要望があれば対応します（異なるCRSへの対応、色付けなど）
※Linux版もあります。

2013年6月29日土曜日

生物多様性地域戦略を策定するためのタウンミーティングに行ってきた

今日は高知県の生物多様性地域戦略を策定するためのタウンミーティング（PDF）＠物部川流域に行ってきた。

参加者は30人ぐらい。
個人的には、会の趣旨が読み取れず消化不良でした。

せっかくなので気の利いたことを発言しようとしたけどあえなく撃沈....
やっぱり話術のスキルが足りません....
発言の練習不足です。

MABでも思うのですが、こういう話は生き物屋さんだけでは限界がある。
そればかりを考えていた会でした。

2013年6月26日水曜日

X200

ちかごろメインマシンのLenovo X200の調子が悪い...
今日はやたらと無線LANの接続が切れる...
ファンエラーで起動しないこともあるが、これはファンを掃除すれば解決する。

2009年の4月頃に購入したので丸4年
あと3年ぐらいは使い込みたいのだが

ジャンクを購入して自前でパーツ交換すればいけるのだろうか、今はその余裕すらない

持ってくれることを願いつつ、手持ちの無線LANアダプタ（PLANEX GW-USValue-EZ）を装着して急場をしのぐ
これは以前、アクセスポイントを作成しようとして購入したのだが、結局使わずにしまわれていたもの。

USBに差し込んで再起動するとubuntuでも問題なく動いてくれる。

ありがたい

2013年6月18日火曜日

Rのコンソールの出力をカラー表示にするパッケージ

パソコンのシステム入れ替えにともなって、いつものようにR用のエディタ（VIM）も再設定。

その際、Rのコンソールの出力をカラー表示にするパッケージを発見。

なにやらCRANのポリシーに反するそうで、インストールは以下のようにR上で直接ダウンロード＆インストール

download.file("http://www.lepem.ufc.br/jaa/colorout_1.0-1.tar.gz", destfile = "colorout_1.0-1.tar.gz")
install.packages("colorout_1.0-1.tar.gz", type = "source", repos = NULL)

パッケージをロードすると、自動でカラー出力になる。

まあ.....ちょっと気持ち悪い気もするが.....いや気持ち悪い

カラーハイライトはエディタだけでいいかな

2013年6月16日日曜日

二日間消耗してようやくCUDA5 & pyCUDA on Debian wheezy

研究費でnVidiaのGPUを購入出来たので、前々から挑戦してみたかったGPGPUプログラミングに挑戦して見ることに。

ついでにSSDとメモリも購入したいので、OSごと再インストール。

メインのOSはこれまでububtuを愛用していたのだけれど、ubuntuを使い続けることに不安を感じてきたので、ここらで他方面に乗り換えてみることに。

候補はいろいろあるけれど、Linux Mintが好評だけどベースのひとつはubuntuだしちょっと気がひけるので、Debianに乗り換えてみる。
まあubuntuの親みたいなものだし、これまでの知識も無駄ではないだろうということで....

ちょうどDebian 7 wheezyもリリースされたところだし!!

んで、インストールしてみた印象は、やっぱりubuntuはデフォルトでも十分使いやすいディストリなんだなと改めて感心した。
Debianだといろいろと追加しないと行けないので、デフォルトでは使いづらい。

でも安定感はこちらも負けていない。
少し前まではsidのリポジトリを追加しないと心もとなかったけど、今のところは不自由していない。
Debianは初心者なので、これからぼちぼち勉強していこう。

さて、CUDA
ubuntuのリポジトリにはCUDA5が上がっているみたいだけど、wheezyはおそらくcuda4のままみたい。sidにはCUDA5が登録されているようだけど.....

ということで、ソースからインストールすることにしたのだけど......

しかし.....

ただドラーバーをインストールするだけで二日間も消耗してしまった....

ちょっとばてたので、インストール方法は元気になったら載せようと思う。

2013年6月11日火曜日

Terrain Ruggedness Indexなるもの

Gdalツールのgdaldemには簡易ながら地形解析をするオプションがある。

hillshade
slope
aspect
color-relief
Terrain Ruggedness Index (TRI)
Topographic Position Index (TPI)
roughness

使い方はマニュアルの通り

gdaldem [オプション] [入力するDEMファイル名] [出力ファイル名]

オプションにはslopeとかaspectとか、TPIとかを指定する。

これは、QGISの[ラスタ]→[解析手法]→[DEM（テリアンモデル）]にも組み込まれているので、コマンドラインだけでなくQGISのGUIからでも操作できる。

ところで、TPIとか他のオプションは馴染みがあるのだが、Terrain Ruggedness Index (TRI)なるものの正体がわからなかったので調べてみた。

Wilson et al. (2007)によると、3セル×3セルのウィンドウの場合、対象となる中心セル$Z_{(0,0)}$のTRI値は

$TRI = \frac{(|Z_{(-1,1)}-Z_{(0,0)}| + |Z_{(0,1)}-Z_{(0,0)}| + |Z_{(1,1)}-Z_{(0,0)}| + |Z_{(-1,0)}-Z_{(0,0)}| + |Z_{(1,0)}-Z_{(0,0)}| + |Z_{(1,-1)}-Z_{(0,0)}| + |Z_{(1,0)}-Z_{(0,0)}| + |Z_{(1,1)}-Z_{(0,0)}|)}{8}$

となる。

任意のウィンドウサイズ（n×n）の場合には、

$TRI (n) = \frac{\sum_{i=-N}^{N}\sum_{j=-N}^{N}|Z_{ij}-Z_{00}|}{(n^{2}-1)}$

となる。
このとき $N=\frac{(n-1)}{2}$である。

式をみればわかる通り、これは標高の二次微分、つまりデジタルラプラシアンに相当する。

2013年6月9日日曜日

Garmin GPSmap62sで国土地理院の地形図を表示する

ガーミンGPSで国土地理院の地形図を表示させる手順についての防備録

1. まずQGISのOpenlayersプラグインを使って電子国土をQIGS上で表示できるようにする
方法は以下を参照

地図はたいへん　QGISで電子国土基本図

ちょっと手間はかかるけど比較的簡単に終了。

2.QGISプラグインからGarminCustomMapをインストール

以上で準備は完了

あとはOpenlayersプラグインから電子国土をロードして、欲しい範囲の地図をGarminCustomMapで切り取るだけ。

1.電子国土基本図をロード

２．GarminCustomMapのアイコンをクリック

適当に名前をつけて保存をクリック

３．そうすると次のような情報画面が表示される

ここで大事なのがrowsとcolumsの部分。

この場合は　595 rows, 1278 colums

これを次の画面で入力する

４．rowsとcolumsに先ほどの値を入力

５．そしてOKをクリック

できた画像をGPSのカスタムマップフォルダにコピー
これで地形図がGPSに表示される。

かんたんかんたん

2013年6月6日木曜日

白髪山に行ってきた

改めて白髪山に行ってきた。
今日はふるさと林道から避難小屋までのルート

白髪山のブナ林

林床は刈りこまれたミヤマクマザサ

ウラジロモミの純林

ウラジロモミは食害がひどいので、みんな腹巻してます。
林床では表土の流亡が目立ちます。

ガリーも徐々に発達しているようで、斜面浸食がいたるところで見られます。

うぅ....

2013年5月26日日曜日

大山沢に行ってきた

久しぶりに大山沢に行ってきた

2013年5月26日秩父市

いつもの通りモニタリングサイト1000の作業のためです。
梅雨前の気持ちの良い時期でした。

相変わらずシカの食害がひどいですが....

2013年5月3日金曜日

さおりが原に行ってきた

昨日に引き続き、山嶺のさおりが原に行ってきた。
こちらも同じように、シカに食い荒らされています。

うぅ....

2013年5月2日木曜日

白髪山（三嶺）にいってきた

今日は研究室のメンバー＆地学の先生と白髪山登山に行ってきた。

ルートは白髪山登山口から白髪山、白髪分岐、ふるさと林道。

白髪山山頂からカヤハゲ、三嶺、アオザレを望む

今回の登山の目的は、シカによる食害の状況の観察。

状況はかなり悪く、山域ほぼ食い荒らされている状況。
上の写真は白髪山山頂の枯れ上がったダケカンバなどの低木林
林床植生は、ミヤマクマザサはまだ残っているけど、その他はほぼ壊滅的な状態。
ウラジロモミは樹皮剥ぎの被害がひどく、枯死木も目立つ。

どこもかしこも....

さらに悪いことに三嶺では林床植生の衰退による表土の流亡、斜面崩壊が進んでいる。

山域全体で土砂災害の発生リスクが高まっている
とにもかくにもシカの個体数をなんとかしないと....

2013年4月19日金曜日

トレンド

Biodiversity and Conservationに面白い論文が出ていた。

Mccallum and Bury,
Google search patterns suggest declining interest in the environment,
Biodiversity and Conservation,
DOI:10.1007/s10531-013-0476-6

内容は、タイトルにある通りGoogleで環境に関連する検索のトレンドを調べたというもの。
データの信ぴょう性はさておき、面白そうだったので、自分でも試してみた。

眉唾だけど、COP10 あたりで生物多様性絡みの検索がピークを示しているは納得できる。

地域性とかあるとおもうんだけどどうだろうか...

2013年4月9日火曜日

Rcpp+inlineを使ってみた

まえから気になっていたRcpp+inlineを使ってみた

単純な繰り返しループの集計計算がかなり時間がかかるので、Rcppでどれくらい改善されるか挑戦

まずは通常のコード

非常に単純で以下のxとBeta,Alphaの線形結合のquantileを求めるというもの。

x <- rnorm(100)
Beta <- Alpha <- rnorm(1000000)
res1 <- sapply(1:length(x), function(i) 
               quantile(exp(Beta*x[i]+Alpha)/(1+exp(Beta*x[i]+Alpha)),probs=0.025))

次にRcpp+inline

library(Rcpp)
library(inline)
src <- "NumericVector xx(x);
  NumericVector beta(b);
  NumericVector alpha(a);
  Function quantile(\"quantile\");
  int xsize=xx.size(), bsize=beta.size();
  NumericVector mu(bsize);
  NumericVector res(xsize);
  for(unsigned long int i=0; i < xsize; i++){
   for(unsigned long int j=0; j < bsize; j++){
    double dummy=exp(beta[j]*xx[i]+alpha[j]);
    mu[j]=0.0;
    mu[j]=dummy/(1.0+dummy);
   }
   res[i]=0.0;
   res[i]=as(quantile(mu, _[\"probs\"]=0.025));
  }
  return res;"

fun.R <- cxxfunction(signature(x="numeric",a="numeric",b="numeric"),
       src,plugin="Rcpp")

結果

x <- rnorm(100)
Beta <- Alpha <- rnorm(1000000)
system.time(res1 <- sapply(1:length(x), function(i) 
    quantile(exp(Beta*x[i]+Alpha)/(1+exp(Beta*x[i]+Alpha)),probs=0.025)))
   ユーザ   システム       経過  
    19.162      1.008     20.169 
system.time(res2 <- fun.R(x,Beta,Alpha))
   ユーザ   システム       経過  
     8.936      0.932      9.870

計算結果も同じだった

summary(res1-res2)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       0       0       0       0

RcppでRの関数を使うのではなく、GNU Scientific Library (GSL)のquantile(gsl_stats_quantile_from_sorted_data)も試してみたけどこちらは通常のRの1.5倍ぐらいの時間がかかった....
何がボトルネックなんだろう...

2013年4月8日月曜日

R　メジャーバージョンアップ

知りませんでしたが、Rのバージョンがアップされたようです。

統計解析ツール「R」、8年半ぶりのメジャーバージョンアップ版「R 3.0.0」リリース

ニュースをみるといろいろと変わっているようです。

2013年1月8日火曜日

アクセプト

昨日出してた論文がEcological Indicatorsにアクセプトされました♪

本当に嬉しい！

突っ込みどころ満載の内容（専門家の人にはいろいろと批判をされていた）ので、嬉しいよりも安心したという気持ちのほうが強いですね~

今回は、1回目の論文改訂にかなり時間をかけてしまいました（この悪い癖を直さなければ....）

これがなければ年内には受理されていたかも

6月25日　投稿

6月28日　受付

8月27日　１回目のコメント

11月23日　改訂原稿返送

12月18日　2回目のコメント

1月7日　改訂原稿返送

1月8日　受理

2013年1月7日月曜日

これで最後？

年末にレビューアーから返ってきていた論文を、修正して返しました。

以前から投稿していたやつです。

今回は、こまごまとした修正だけだったのですが、エディターに例のごとく英文校閲をうけろと言われました。

しかたがないですね。

今回もテキストチェックに出しました。同一原稿で２回目のチェックだったので、かなり割り引いてもらえました。

さて、良い便りがくるとよいですが.....

生態学会の要旨登録

生態学会の要旨を登録しました.....

が、今年はなんとも中途半端な内容で、方法までしか書けませんでした....

これから解析を頑張らねば!

登録: 投稿 (Atom)

A Big Fish in a Little Pond (もきログ)

ページ