[PR] この広告は3ヶ月以上更新がないため表示されています。
ホームページを更新後24時間以内に表示されなくなります。

memo: How to make watershed map

流域界や河道のデジタルマップを作る方法のメモ

1. DEMをダウンロードする。

国土地理院の基盤地図情報ダウンロードサービス： http://fgd.gsi.go.jp/download/ から、必要なあたりのDEMをダウンロードする（JPGIS(GML)形式の奴がいい）。とりあえず、すべてのデータをダウンロードしてみた。（結構大変）

2. GeoTIFF形式への変換

まず、FG-GML-*-DEM10B.zipを解凍展開（FG-GML-*-*-dem10b-*.xmlができる）

---
#!/bin/sh
for zipfile in FG-GML-*-DEM10B.zip
do
  unzip ${zipfile}
done
---

つぎに、XMLをTIFに変換

---
for xmlfile in FG-GML-*.xml
do
  gdal_translate -of GTiff ${xmlfile}.bin ${xmlfile}.tif
done
---

3. 適当な領域でマージする

オリジナルの地図の区分だとちょっと細かすぎてわけわかんないので、ひとつのzipファイル分でマージしてやる。

---
for x in 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
do
for y in 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
do
gdal_merge.py -o ${x}${y}.tif -ps 0.00011111111111 0.00011111111111 -init -9999 FG-GML-${x}${y}-??-dem10b-????????.xml.tif
done
done
----

とやると、こんな感じのタイルになる→ 日本全体（一部地域を除く）

4. 目的の流域周辺のDEMを取り出す

とりあえず、GRASSでtifを読み込んで表示させてあたりをつける。

----
r.in.gdal --o input=${x}${y}.tif output=${x}${y}
r.null map=${x}${y} setnull=-9999
r.colors map=${x}${y} rules=srtm
g.region rast=${x}${y}
d.rast ${x}${y}
----

この結果が上の日本全体（一部地域を除く）という訳です。
これで当たりをつけて、必要な部分を適用にマージする（悲しいことに由良川流域ｊは４枚にまたがってる）

----
gdal_merge.py -o yura.tif -ps 0.00011111111111 0.00011111111111 -init -9999 5334.tif 5335.tif 5234.tif 5235.tif
----

出来たyura.tifをr.in.gdalして、d.zoomで必要部分だけズームアップして、さらに、その部分だけのファイルを作る。

----
r.in.gdal --o input=yura.tif output=tmp
r.null map=tmp setnull=-9999
r.colors map=tmp rules=srtm
g.region rast=tmp
d.rast tmp
d.zoom
（適当にズームアップ）
r.out.ascii input=tmp output=tmp.dat
r.in.ascii input=tmp.dat output=yura
g.remove tmp
----

という感じで表示部分だけ切り取る。（r.outしてr.inするあたりは、チョット・ヲイヲイだけど、他の方法がよく分からん。きっと、スマートな方法があるに違いないが）

5. いよいよ流域界と河道の抽出

GRASSのコマンドは、r.watershed を使うのが良さそう。さてどう使うんだ？・・・つづく（20110408)
以下、ちょっと練習（20110408)

d.rast ashu
r.watershed elevation=ashu stream=ashu_stream threshold=10000
r.thin input=ashu_stream output=ashu_stream_thin
r.to.vect input=ashu_stream_thin output=ashu_stream_vect feature=line
d.vect ashu_stream_vect

というわけで、つづき。(20110411)

-----
r.watershed elevation=ashu accumulation=ashu_accu basin=ashu_basin threshold=10000
r.to.vect input=ashu_basin output=ashu_basin_vect feature=area
d.vect ashu_basin_vect type=boundary
-----

由良川上流部

由良川流域全体

めでたしめでたし。