Rのテキスト処理


ご無沙汰です。
Mwクwリ a.k.a プレゼン作りはデザインから派 です。

 支部会で使うデータで地形の凹凸度を出す必要があったのですが、出しようがなかったのでRで書いてみることにしました。

 前に九州全部のDEMとかで組んだときはPerlだったんですが、コードが長くなっちゃうのと、ループのループのループみたいな構造になっちゃうのがちと嫌。

 んでどうにか完成。半日くらいかかっちゃいました。やっぱ、sum()とかで、ベクトルの合計とかができるのは便利。
 コードもだいぶ短く書けますし、無駄なループもないので処理も速い(多分)。
 デメリットは、前にも書きましたが、データ構造の行き来がイマイチわからん。

 例えば今回もデータフレームから取り出した数値をつなげていってベクトルにしても、ベクトルにならないという、とんちのようなトリックを理解できるのに2時間ほどかかってしまいました。

 これはどゆことかいうと、データフレームから数値を取ってきてもそれはベクトルでなく、データフレームになっていて、データフレームをつなげていくと、それはベクトルではなく、リストになってしまうのですね~~
 たいした違いないかと思うと、コレが恐ろしいくらい違うわけで、例えば・・・・

   hogehoge <- hogehoge + 1

 このオブジェクトがベクトルだと、ベクトルの要素全部に1が足されるのですが、リストだとエラーを出しやがります。
というわけで、今回のスクリプトの肝は、リストをベクトルに変換するunlist()。この関数で変換すれば、普通のベクトルのように使えます。
オブジェクト指向型の功罪というか、なんというか・・・
こういった類のトラップは往々にしてあるわけでそゆのを回避しながらコーディングするのがなかなか一苦労なわけでして・・・ま、こういうコードを財産にしていけばそういう悩みも減っていくと信じて、しばらくはRでいろいろ書いてみようかと思います。

##凹凸度計算スクリプト
##周囲10セルの高度を参照して計算
##最初6行はヘッダなのでスキップ
########################################
setwd(“d:/【研究】葉枯れ/地図データ”)
dat<-read.table(“dem.txt”,header=F,skip=6)

##最初の10行はバッファなので、すべてを0にする。
roughness <- dat[FALSE,]
for( i in 1:10){
  col <- dat[i,]
  col <- col*0
  roughness <-rbind(roughness,col)
}

##周囲10セルの値から凹凸度を計算
##行頭10と行末10はバッファなので0にする。
for( j in 11:(nrow(dat)-10)){
v.point.rough <- c(0,0,0,0,0,0,0,0,0,0)
  for( p in 11:(ncol(dat)-10)){
    v.rough <- NULL
    rel.rough <- NULL
    point.rough <- NULL
    v.rough <- c(dat[c(j-10),(p-10):(p+10)],dat[c(j-9),(p-10):(p+10)],dat[c(j-8),(p-10):(p+10)],dat[c(j-7),(p-10):(p+10)],dat[c(j-6),(p-10):(p+10)],dat[c(j-5),(p-10):(p+10)],dat[c(j-4),(p-10):(p+10)],dat[c(j-3),(p-10):(p+10)],dat[c(j-2),(p-10):(p+10)],dat[c(j-1),(p-10):(p+10)],dat[c(j),(p-10):(p-1)],dat[c(j),(p+1):(p+10)],dat[c(j+1),(p-10):(p+10)],dat[c(j+2),(p-10):(p+10)],dat[c(j+3),(p-10):(p+10)],dat[c(j+4),(p-10):(p+10)],dat[c(j+5),(p-10):(p+10)],dat[c(j+6),(p-10):(p+10)],dat[c(j+7),(p-10):(p+10)],dat[c(j+8),(p-10):(p+10)],dat[c(j+9),(p-10):(p+10)],dat[c(j+10),(p-10):(p+10)])
    v.rough <- unlist(v.rough)
    target <-unlist(dat[j,p])
    rel.rough <- v.rough – target
    point.rough <- sum(rel.rough)
    v.point.rough <- c(v.point.rough,point.rough)
  }
buff <- c(0,0,0,0,0,0,0,0,0,0)
v.point.rough <- c(v.point.rough,buff)
roughness <- rbind(roughness,v.point.rough)
}

##最初の10行はバッファなので、すべてを0にする
for( i in 1:10) {
  col <- dat[i,]
  col <- col*0
  roughness <-rbind(roughness,col)
}
 write.table(roughness, “d:/【研究】葉枯れ/地図データ/roughness.txt”, quote=F,col.names=F, append=T,row.names=F)

  • イニシャルB

    SECRET: 0
    PASS: 74be16979710d4c4e7c6647856088456
    エクセルとかいうソフト、意外と便利ですよ。

    今研究室も空前の露出度ブームです。このコマンド頂くかもしれません。

  • myuhe

    SECRET: 0
    PASS: 74be16979710d4c4e7c6647856088456
    エクセルなんてヒューマンエラーの温床な屎ソフト。ざっと見るくらいで使うが吉かと。
    今回のスクリプトは肝心の凹凸度出す辺りがかっこ良くなくて嫌だったんすが、良い知恵も浮かばず、そんまま強引にやっちまいました。
    露出度やるときは、判定範囲をデータフレームで切り取って、判定かけると格好いいのが書けるかも。
     また、なんかの機会に作ってみます。