ページ

2011年3月3日木曜日

GRASSのr.univarの出力をラスター演算(r.mapcalc)で使うには

今回もマニアックな内容です。

やりたいこと

GRASSのr.univarが出力する,あるラスターデータの最大値や最小値を使ってラスター演算(r.mapcalc)をする

方法
ターミナルユーザーであれば,awkを使うことで,問題なく処理をすることができる

あるラスターデータを "dummy" とする。
r.univarを実行すると,以下のように出力される。

$r.univar map=dummy
total null and non-null cells: 5157441
total null cells: 0
Of the non-null cells:
----------------------
n: 5157441
minimum: 0
maximum: 1
range: 1
mean: 0.0155387
mean of absolute values: 0.0155387
standard deviation: 0.123682
variance: 0.0152973
variation coefficient: 795.961 %
sum: 80140

出力結果を見てみると,最小値は7行目,最大値は8行目のそれぞれ2列目に出力される。

ということで,最大値のみを取り出すには次のようにする。
NRで取り出す行数を,\$で列数を指定している。
r.univar map=dummy | awk '{if(NR==8){printf(\$2); printf(RS)}}'

最小値の場合にはNRの部分を7にすればよい。
r.univar map=dummy | awk '{if(NR==7){printf(\$2); printf(RS)}}'

さて,値自体は取り出すことができたので,次はラスター演算(r.mapcalc)で使えるようにしたい。
単純に,取り出した値をシェル変数に代入する。
試しに,割合(%)を示すデータdummy2を作成してみる。
i=`r.univar map=dummy | awk '{if(NR==8){printf(\$2); printf(RS)}}'`
r.mapcalc "dummy2=(dummy/$i)*100"


これ以外にも方法はあると思う。単純にr.mapcalcの構文内に`r.univar map=dummy | awk '{if(NR==8){printf(\$2); printf(RS)}}'`を入れるだけでも良いかもしれない。

応用
次のようにすることで,r.univarの出力結果をテキストファイルに1行で書き出すことができる。
r.univar map=dummy | awk '{
if(NR==6){printf(\$2); printf(FS);}
else if(NR==7){printf(\$2); printf(FS);}
else if(NR==8){printf(\$2); printf(FS);}
else if(NR==9){printf(\$2); printf(FS);}
else if(NR==10){printf(\$2); printf(FS);}
else if(NR==11){printf(\$5); printf(FS);}
else if(NR==12){printf(\$3); printf(FS);}
else if(NR==13){printf(\$2); printf(FS);}
else if(NR==14){printf(\$3); printf(FS);}
else if(NR==15){printf(\$2); printf(FS);}}
END{printf(RS)}' > dummy2

これは非常に使える方法で,例えば沢山のラスターデータの統計量をひとつのテキストファイルにまとめることができる。
結果を書き込むファイルを "result.txt" とする
#結果を書き込むファイルを作成し,ヘッダを書き込む
echo "n minimum maximum range mean MeanOfAbsoluteValues
StandardDeviation  Variance VariationCcoefficient Sum" > result.txt
#マップセット内のラスターデータすべてのr.univarの出力をresult.txtに書きこむ
for i in `g.mlist rast`
do
r.univar map=\$i | awk \'{if(NR==6){printf(\$2); printf(FS);}
else if(NR==7){printf(\$2); printf(FS);}
else if(NR==8){printf(\$2); printf(FS);}
else if(NR==9){printf(\$2); printf(FS);}
else if(NR==10){printf(\$2); printf(FS);}
else if(NR==11){printf(\$5); printf(FS);}
else if(NR==12){printf(\$3); printf(FS);}
else if(NR==13){printf(\$2); printf(FS);}
else if(NR==14){printf(\$3); printf(FS);}
else if(NR==15){printf(\$2); printf(FS);}}
END{printf(RS)}' >> result.txt
done

GRASSでもこの辺りを充実してくれると嬉しいのですが...

2011年2月22日火曜日

GPSゲット!!


新しいGPS(GPSmap 62S)を購入しました。
右側ののやつが先代,左側のフェイスが最新型!!

英語版なので,4万2000円!!
お手頃価格ですね~!!

早速今週末の調査から導入しよう!!

2011年2月14日月曜日

誤ってフォーマットしたデジカメのメモリーカードを復旧!!

もう半年も昔に誤ってフォーマットしてしまったデジカメのメモリーカードをようやく復旧しました!!

これぐらい復元できることはわかっていたので,いつか復元しようと思っていたのですが,気がついたらもう半年もたってしまいました...

またいつか必要になるかもしれないので,念のため手順を書いておきます。
本当は必要にならないと良いのだけど...
このページを参考にされる方は,自己責任でお願いします。



まずは必要なソフトのインストール!!
今回はPhotoRecを使用した。
このソフトはubuntuのリポジトリにも登録されているので端末から以下のコマンドでOK!


sudo aptitude install testdisk

インストールしたら起動

photorec

すると以下の画面が表示される。
なになに,管理者権限で実行しなければならないとのこと
実行するには[sudo]に進む。











するとパスワードの入力が求められる。
パスワードを入力しEnter














復元したいメディアを選択する。
内蔵ハードディスクと,USBフラッシュメモリが表示されている。
ここは必ず慎重に確認してから進む












次にパーティションテーブルのタイプを選択する。今回は「Intel」を選択した。














次に復元したい領域を選択する。














ファイルシステムを選択する。FATなので「Other」を選択。














次に,解析対象を空き領域のみにするか,それともディスク全体にするか選択する。ここでは全体を選択。












次に復元したファイルを格納するディレクトリを選択する。













そしたら後は自動で解析が進む















ファイルサイズにもよるが,今回は結構早く終わりました














実際に指定した復元用のディレクトリを見ると,何かしら復元されているのがわかります。












ディレクトリを開いてみるとこのとおり。
ちゃんと復元されています。
ただ,ファイル名は変わってしまっていますが....












でもファイルのプロパティを確認してみると画像情報はちゃんと残っていました。



















すごいですね~
作業も楽で,とても便利です。

さて,今回は説明していませんが,作業は実際にはddこまんどでメモリーカードのデータをUSBフラッシュに移してから,USBフラッシュにたいして復元作業を行っています。
いろいろ検索してみるとわかりますが,そのほうが安全です。

2011年2月2日水曜日

RとFortranの演算速度の比較

そろそろFortranもお勉強しようかとおもい、プログラムをしてみることに!!

作業は、2地点間の距離と角度を求めるというシンプルなもの
ただ、データ数が多く、組み合わせが全部で22,232,390もある。

まずRでプログラミング
特に変なことはしていないと思いますが...

Tree<-read.csv("TreePlot",header=T)
IncW3<-read.table("IncW3Th10S3",header=T,sep=" ")
IncW5<-read.table("IncW5Th10S3",header=T,sep=" ")
IncW11<-read.table("IncW11Th10S3",header=T,sep=" ")

Tree<-as.matrix(Tree[,1:5])
IncW3<-as.matrix(IncW3)
IncW5<-as.matrix(IncW5)
IncW11<-as.matrix(IncW11)

DistW3<-matrix(rep(NA,nrow(Tree)*nrow(IncW3)*7),ncol=7)
DistW5<-matrix(rep(NA,nrow(Tree)*nrow(IncW5)*7),ncol=7)
DistW11<-matrix(rep(NA,nrow(Tree)*nrow(IncW11)*7),ncol=7)

for (i in 1:nrow(Tree)){
 for (j in 1:nrow(IncW3)){
 DistW3[i*j,1:4]<-c(Tree[i,4:5],IncW3[j,1:2])
 DistW3[i*j,5]<-sqrt((Tree[i,4]-IncW3[j,1])^2+(Tree[i,5]-IncW3[j,2])^2)
 DistW3[i*j,6]<-atan((Tree[i,5]-IncW3[j,2])/(Tree[i,4]-IncW3[j,1]))
 DistW3[i*j,7]<-atan((Tree[i,5]-IncW3[j,2])/(Tree[i,4]-IncW3[j,1]))*180/pi
 }
 for (j in 1:nrow(IncW5)){
 DistW5[i*j,1:4]<-c(Tree[i,4:5],IncW5[j,1:2])
 DistW5[i*j,5]<-sqrt((Tree[i,4]-IncW5[j,1])^2+(Tree[i,5]-IncW5[j,2])^2)
 DistW5[i*j,6]<-atan((Tree[i,5]-IncW5[j,2])/(Tree[i,4]-IncW5[j,1]))
 DistW5[i*j,7]<-atan((Tree[i,5]-IncW5[j,2])/(Tree[i,4]-IncW5[j,1]))*180/pi
 }
 for (j in 1:nrow(IncW11)){
 DistW11[i*j,1:4]<-c(Tree[i,4:5],IncW11[j,1:2])
 DistW11[i*j,5]<-sqrt((Tree[i,4]-IncW11[j,1])^2+(Tree[i,5]-IncW11[j,2])^2)
 DistW11[i*j,6]<-atan((Tree[i,5]-IncW11[j,2])/(Tree[i,4]-IncW11[j,1]))
 DistW11[i*j,7]<-atan((Tree[i,5]-IncW11[j,2])/(Tree[i,4]-IncW11[j,1]))*180/pi
 }
}

結果
ユーザ   システム       経過  
  1522.730      1.510   1528.272 

約25分ぐらいですか

次にFortranの場合
コードはこちら
program DistCalc
 implicit none
 real,dimension(1565,5)::Tree
 real,dimension(7146,3)::IncW3
 real,dimension(4952,3)::IncW5
 real,dimension(2108,3)::IncW11
 real pi
 real DistW3(1565*7146),AngW3(1565*7146)
 real DistW5(1565*4952),AngW5(1565*4952)
 real DistW11(1565*2108),AngW11(1565*2108)
 integer i,j 

 open(10, file='Tree')
 open(20, file='IncW3')
 open(30, file='IncW5')
 open(40, file='IncW11')
 read(10,*) ((Tree(i,j),j=1,5),i=1,1565) !iは行、jは列
 read(20,*) ((IncW3(i,j),j=1,3),i=1,7146)
 read(30,*) ((IncW5(i,j),j=1,3),i=1,4952)
 read(40,*) ((IncW11(i,j),j=1,3),i=1,2108)

 do i=1,1565
  do j=1,7146
  DistW3(i*j) = sqrt((Tree(i,4)-IncW3(j,1))**2 + (Tree(i,5)-IncW3(j,2))**2)
  AngW3(i*j) = atan2((Tree(i,5)-IncW3(j,2)),(Tree(i,4)-IncW3(j,1)))
  end do
 end do

 do i=1,1565
  do j=1,4952
  DistW5(i*j) = sqrt((Tree(i,4)-IncW5(j,1))**2 + (Tree(i,5)-IncW5(j,2))**2)
  AngW5(i*j) = atan2((Tree(i,5)-IncW5(j,2)),(Tree(i,4)-IncW5(j,1)))
  end do
 end do

 do i=1,1565
  do j=1,2108
  DistW11(i*j) = sqrt((Tree(i,4)-IncW11(j,1))**2 + (Tree(i,5)-IncW11(j,2))**2)
  AngW11(i*j) = atan2((Tree(i,5)-IncW11(j,2)),(Tree(i,4)-IncW11(j,1)))
  end do
 end do

 pi=4*atan(1.0)

 open(50,file="DistW3")
 do i=1,1565
  do j=1,7146
  write(50,*) Tree(i,4), Tree(i,5), IncW3(j,1), IncW3(j,2),&
   DistW3(i*j), AngW3(i*j), (AngW3(i*j)*180)/pi
  end do
 end do

 open(70,file="DistW5")
 do i=1,1565
  do j=1,4952
  write(70,*) Tree(i,4), Tree(i,5), IncW5(j,1), IncW5(j,2),&
   DistW5(i*j), AngW5(i*j), (AngW5(i*j)*180)/pi
  end do
 end do

 open(80,file="DistW11")
 do i=1,1565
  do j=1,2108
  write(80,*) Tree(i,4), Tree(i,5), IncW11(j,1), IncW11(j,2),&
   DistW11(i*j), AngW11(i*j), (AngW11(i*j)*180)/pi
  end do
 end do

end program DistCalc

で結果はというと、一目瞭然ですね!!
$ f95 DistCalc.f90
$ ./a.out

real 3m20.994s
user 3m14.140s
sys 0m5.330s

今回の場合には、Fortranのコードを書くだけで結構時間がかかった(まだ初心者なので)ので、Rだけで作業しても良かったかもしれません。
ですが、今後のことを考えてもFortranをお勉強したほうが良いですね〜
できればC++もお勉強したいのですが...

2011年1月21日金曜日

GRASSのr.mapcalcのifの取り扱い

r.mapcalcにおけるifの取り扱いについて
結構忘れやすいのでメモしておく

マニュアルには以下の4つが記載されている

1)     $if(x)$
2)     $if(x, a)$
3)     $if(x, a, b)$
4)     $if(x, a, b, c)$
  1. もし$x$が0(ゼロ)でない場合は1を、それ以外は0(ゼロ)とする
  2. もし$x$が0(ゼロ)でない場合は$a$を、それ以外は0(ゼロ)とする
  3. もし$x$が0(ゼロ)でない場合は$a$を、それ以外は$b$とする
  4. もし$x$が0(ゼロ)より大きい場合場合は$a$を、0(ゼロ)の場合は$b$を、0(ゼロ)未満の場合は$c$とする

$x$の部分には比較演算子が利用可能。
また、ifのなかにifを使うこともできる。

1)     $if(x \textgreater 5, 1)$
2)     $if(x \textgreater 5, 1, if(x \textgreater 10, 2))$

  1. もし$x$が5より大きい場合は1を、それ以外は0(ゼロ)とする
  2. もし$x$が5より大きい場合は1を、10より大きい場合は2を、それ以外は0(ゼロ)とする

2011年1月20日木曜日

シカ襲来

ここ数日、(妄想の中で)シカの襲来にあっていました
というのも、締切り間際になって森林学会の要旨を登録するための解析をしていたため!!
今度の森林学会(ちなみに初参加)ではシカの食害の発表をします。

要旨を書き始めたのは、昨日のアフターファイブ(←もしかして死語?)を過ぎてから...
でも、結果はすでに出ていたので、なんとか昨日のうちに書き終えて共同研究者に送ることができました。

にわかシカ研究ですが、よかったら聞きに来てください。

2011年1月7日金曜日

subversionで一括でファイルを追加または削除

Subversionで一括でファイルを追加または削除についての防備録

ぐぐるといろいろと情報がでてくるけどファイル名に空白が含まれる場合には対応できない

ということで自作する。

追加する場合
svn status | grep ^? | sed 's/?//g' | awk '{sub(/^[ ]+/,"");print "#",$0,"#"}' | sed -e "s/ #/'/g" | sed -e "s/# /'/g" | xargs svn add
svn commit -m "Data added!!"

削除する場合
svn status | grep ^! | sed 's/!//g' | awk '{sub(/^[ ]+/,"");print "#",$0,"#"}' | sed -e "s/ #/'/g" | sed -e "s/# /'/g" | xargs svn delete
svn commit -m "Data deleted!!"

やたらパイプが多いですね....
まあ,初めからファイル名に空白を使わなければ良いのですが...