基盤地図情報25000の等高線で計曲線・主曲線・補助曲線を描き分ける

GIS

自然環境調査などで使う印刷用3次メッシュマップのQGISを使った作り方（http://d.hatena.ne.jp/kooi/20110619）の補足です。
QGISの地図プロジェクトに読み込んだベクタレイヤの記号は今のところ「共通シンボル」で設定されています。これを「固有値」を使って属性によって記号を変化させるという話です。
基盤地図情報25000の等高線データ(Cntr.shp)の属性ファイル(Cntr.dbf)には「標高」という列があります。この標高フィールドを使って等高線を色分けしてみます。プロパティのシンボルで凡例タイプに「固有値」を選び，分類フィールドを「標高」にします。「分類」のボタンを押すと自動で全ての属性値に個別の凡例をつけることができるようになります。（「標高」のフィールドは数値ではなくテキストでデータが入っています。）

しかしこれを1つ1つ設定するのはかなり面倒です。10，20，30，40，60…と主曲線にあたる標高は「クラスを削除」してデフォルトで一括指定するとしても，50mおきの計曲線だけでも20本以上あります。緩傾斜のところで5mごとに入ってくる補助曲線も何本もあり，それぞれ定義していたらやはり大変です。そこで，ちょうど属性ファイルの「表示区分」フィールドは空欄になっているようなので，ここに等高線の種類を入れて表示区分で凡例を分類することにします。

  • QGISの属性テーブルを編集モードにしたときに使える「フィールド計算機」では，まだ複雑な計算は難しいようです。http://www.osgeo.jp/user_guide/user_guidech3.html#x10-1750003.7
  • Excel2003まではdbf形式を扱うことができたので，Excelで開いて表示区分のG2セルに以下の式を入れて列へコピーすることが考えられます。しかし，Excel2003までは65535行しか扱えないので，福岡県全部（約29万行）だとオーバーフローします。Excel2007ではdbf形式で保存できないのでこの方法はできません。基盤地図情報閲覧・コンバートソフトでシェープファイルにするときに少しずつ（ダウンロードしたzipファイル1つずつ）別のシェープファイルに変換すればExcel2003でオーバーフローせずにできるでしょう。

=if(right(I2,1)="5"),"2",if(or(right(I2,2)="00",right(I2,2)="50"),"1","0"))

今回は，約29万行のdbfファイルを一発で扱うために，またしても統計解析ソフトウェアRを使います。コードは以下のようになります。dbfはいちおう別名で保存していますので，できたdbfファイルを確認してから，名前をCntr.dbf→Cntr_ori.dbf，Cntr_R.dbf→Cntr.dbfというように変更してください。

library(foreign)
CntrType <- function(x) {
y <- as.numeric(as.character(x))
ifelse(floor(y/50)==y/50,"1",ifelse(floor(y/10)==y/10,"0","2"))
}
dd <- read.dbf("C:/GIS/Cntr.dbf")
dd$表示区分 <- CntrType(dd$標高)
write.dbf(dd, "C:/GIS/Cntr_R.dbf")

Cntr.dbfを置き換えたら，表示区分のフィールドに主曲線は"0"，計曲線（50mごと）は"1"，補助曲線(5m）は"2"が入っているので，それぞれ適切な記号を定義します。

  • R（参考RjpWiki http://www.okada.jp.org/RWiki/）のベクトル演算は強力で，上記のコードではデータの読み込みや書き出しに比べたらあっという間に計算は終わります。関数中でif文ではなくifelseを使ってベクトル演算できるようにしておくのがミソです。
  • 個人的には一番手間取ったのはdbfを読み込んだデータフレームで標高の数値がただのテキストではなくてfactorだったので，as.characterをかまさないと期待した結果を返してくれなかったことでした。

このようにして作った地図（pdfに出力したもの）の一部です。