AvalancheStartingZone②(雪崩発生地点の算出)

斜度・植生・積雪深を合計。AvalancheStartingZoneという名前は最近使ってない。論文格上ではわかりにくいと思って。

斜度と植生だけのときと一緒。



QGISでGoogleMapと重ねてみた。QGISは表示の変更をするとすぐ落ちるのが嫌。



なんかずれてる。

拡大・縮小してたら直った。

積雪深データ作成

時間が無いのでスクリーンショットだけ。

積雪深フィールドワークデータをエクセルに格納してArcに追加。shpにエクスポートして空間内挿。



なぜかクリギングは一度しかうまくいかない。IDWでやる。



オプションいじればポイントデータの外でも空間内挿できるよ。



ほら。

AvalancheStartingZone①―植生と斜度の合計

植生データ作成⑤―ラスタ変換

ラスタ変換の前に座標系をWGS84UTM54に変換。

ラスタ変換。

入力フィーチャ→植生shpデータ、値フィールド→point、出力ラスタデータセット→適当な名前、セルの割り当て方法→CELL_CENTER、優先フィールド→NONE、セルサイズ→傾斜ラスタデータかDEMにしてOK。

結果はssとってない！

植生データ作成④

Pythonでも辞書型内に改行があっても大丈夫じゃないの？と思って以下のようにしました。

[sourcecode language="python"]
import csv

def csv2dic(readed):
for row in readed:
print ','.join(row)+','

test = csv.reader(open('csvファイル', 'rb'))
csv2dic(test)
[/sourcecode]

結果↓


出力結果をターミナル上でコピーしてテキストエディタかなんかに貼付け。前と後を編集。

[sourcecode language="python"]
dic={1001:9,
1002:9,
1003:9,
（中略）
9998:10,
9999:10}
[/sourcecode]

最後にprint dic[9925]とかなんとかつけて.pyで保存すれば実行した時に10とか返ってきます。

ArcMapを起動してフィールド演算をします。
フィールド演算ウィンドウはレイヤウィンドウでshpファイルを右クリックして「属性テーブル」を開き、「新規フィールド」(だったかな？)を選択、型をshort、名前はpointとかの任意の名前してOK。作成されたフィールドの上の名前を右クリックして「フィールド演算」

形式で"Python"にチェックを入れて、"コードブロックを表示"にもチェック。上のテキストボックスに先ほど作った辞書型をコピペ、下のテキストボックスにはdic [!MAJOR1!]を入力して「OK」
計算が始まります。

できた。

色をつけてみました。

---

12/7 追記
入門 Unix for Mac OS X 第4版読んでて知ったんですが、別にターミナルの出力をコピペなんてことしなくても
python csv2dic.py > XXXX.txt
とかすればいいんですね。

植生データ作成①

新しい評価基準に基づき、植生ファイルを修正。ArcPyで県別植生shpファイルを1次メッシュで切り取る。


植生shpファイルは例によって生物多様性センターから新潟と群馬をダウンロード。

生物多様性センター(環境省 自然環境局)
http://www.biodic.go.jp/

1次メッシュはESRIのサイトから。

標準地域メッシュ作成ユーティリティ
http://www.esrij.com/products/arcview3/faq/file/make/makemesh.html

ダウンロードした1次メッシュshpファイルはなぜか座標系が指定されてないのでArcCatalogでファイルを右クリック>「プロパティ」からWGS1984に指定してあげる。


1次メッシュshpファイル、都道府県界shpファイル、新潟県植生shpファイルを表示させたとこ。グンマーいらなかったですね。


1次メッシュshpファイルの必要なところだけ選択してエクスポート。



ArcGISでの作業はここまで。本当はこういうのもArcPyでできるんだろうなあ。
次はいよいよArcPyの出番。コードはESRIのサポートページにログインしてQ&Aで"サンプルコード"とかで検索すれば出てくると思います。"Python スクリプトを ArcGIS Desktop を起動せずに実行する方法"という記事です。

Esri製品サポート|サポート|ESRIジャパン株式会社
http://www.esrij.com/support/esri/

エクスポートした1次メッシュshpファイルと新潟の植生shpファイルを同じディレクトリに入れてpyファイルにパスを指定、それぞれのshpファイルを指定してメインメニュー「Run」>「Run Module」

失敗したときはエラーメッセージが表示されますが、成功時は何も表示されません。結構時間かかります。

できたファイルをArcGISで確認。

クリップできてました。

次は植生の分類をpythonで書いて楽にやりたい。
DSMが手に入らなかったから○○群生と○○・△△群生とかで点数変えるかー。

植生の内訳まとめ
基礎調査目次
http://www.biodic.go.jp/kiso/fnd_f.html
3次メッシュ植生データのリンクから。

傾斜データ修正

傾斜データに抜けがあるので作成しなおし。

概ね前回と一緒。

今回は　“再分類されなかった値をNoDataに変更”のチェックを入れて　の手順を行わない。









四隅がじゃまだけどとりあえずデータの抜けはなし。

Surferを試してみる

既存研究がSurferを使っていたというので試してみる。

Surfer 10 | Golden Software
http://www.goldensoftware.com/products/surfer/surfer.shtml

GISの文字が見られないんだけどなんでだろー。
Free Demoをダウンロードしてインストール。立ち上げてみる。

見た目はArcGISやQuantumGISに似てますね。

shpファイルを追加してみたとこ。
メインメニュー「File」>「Import」でshpマップを追加できます。ここでtiffファイルを選択したらエラーが出ました。ラスターデータはダメなのかも。


.tiffのDEMは「Map」>「New」>「New 3D Surface」で追加できました。

「Scale」タブの"Z Scale"の"Length"をいじったところ。かっこいい。

「Map」>「New」>「Image file」なら2Dで表示されます。


ざっと見た感じArcGISのArcToolBoxやQuantumGISのプラグインに該当するものはなさそう。ビューワーに近い感じだからホームページにはGISと書かれてはなかったのかな。

傾斜データ作成

分類基準を変えたのでArcGISで再計算。tifファイル形式のDEMを作るまでは一緒。

「データの追加」からDEM(marge.tif)を追加。灰色一色なのでレイヤプロパティの「シンボル」タブで"ストレッチ"内の"タイプ"を「割合クリップ」に変更して「OK」。一応座標系が想定通りか確認。(WGS84UTM54N)





次に傾斜角を計算。

ArcToolboxから「Special Analyst ツール」>「サーフェス」>「傾斜角(Slope)」で"入力ラスタ"を「(marge.tif)」、"出力ラスタ"に適当な名前(Slope)を入力して他はそのままで「OK」





新しい分類に基づき再分類。

ArcToolBoxから「Special Analyst ツール」>「再分類」>「再分類(Reclassify)」で"入力ラスタ"を「(Slope)」、"再分類フィールド"はそのままで「Value」、「分類」ボタンをクリック。



"分類""クラス"を「3」、"閾値"に(30,40,90)と入力、"データの除外"の「除外」ボタンをクリックして"除外する値"に(0)を入力して「OK」、「OK」



"新しい値"に(4,7,10)と入力、"出力ラスタ"に適当な名前(slope_rec2)、("再分類されなかった値をNoDataに変更"のチェックを入れて)「OK」




右クリックで「データ」>「データのエクスポート」でimgファイル(slope_rec2.img)でエクスポート。Arcで作ると余計なファイルがたくさん生成されるなあ。

QGISで見てみたところ。レイヤプロパティで「スタイル」タブで"カラーマップ"を「カラーマップ」に、「カラーマップ」タブで"エントリー数"を(3)、「等間隔」で「分類」して「OK」にしたところ。ArcでNoDataにしたところがnullとなって描画されていない。

QGISからKMLをエクスポートする(多分続きます)

研究会でポリゴンshpをGoogleMapに表示させたいのだけど、ArcGISではうまくKMLでエクスポートできないと言っている人がいたので、昔やってみようとしたQGISからKMLをエクスポートする方法を再度調べてみる。

すると、9/15にどうやらqgis2googleがQGIS1.8で使えるようになったみたいだ。QGIS1.8? 最新verは1.7だけど…。とりあえず書いてあるとおりにやってみる。
例によってMacは非対応でWindows/Linux限定なので注意。今回はWindowsVistaでやってみた。

GIS-Lab: Открытие объектов из QGIS в Google Maps/Earth:
http://gis-lab.info/qa/qgis2google.html

まず、最新のOSGeo4Wをダウンロードするのだけれども、QGISに問題があるとよろしくないらしいのでアンインストールする。QGIS単体でインストールしている場合はコントロールパネルから削除。OSGeo4Wのアンインストールは面倒で、以下の3ステップ必要らしい。

GIS-Lab: Установка QGIS с помощью OSGeo4W
http://gis-lab.info/qa/qgis-osgeo4w.html

まず、C:OSGeo4Wフォルダを削除。これは簡単。
次がややこしくてスタート→検索バーに"Regedit"と入力。一番上に出てくるRegeditを実行。

HKEY_USERS\S-1-5-21-1294419703-2643849633-1075583040-4140\Software\QuantumGIS を削除。

最後に C:\Users\USERNAME\.qgis を削除。
これでQGISが消えました。

ようやくOSGeo4Wをインストールします。

インストールが完了したらQGISを起動。
OSGeo4Wを起動してQGISと入力。

できた？

ところが…。

QGISのバージョンは1.7でした。
本来C:\OSGeo4W\apps\qgis-dev\pluginsにqgis2googleを入れるわけですが、qgis-devフォルダがないのでqgis-dev\pluginsに入れてみましたが、地物選択するときにエラー落ちするだけでした。
QGIS1.8が一般公開されるまで様子見かなー。

多分続きます。

StrokePointsの作成

雪崩発生地点をベクタポイントデータで追加する。
以前と同じようにtxtファイルを作っていたら面倒くさくなったのでLibreOfficeで作る。

ここまでは良かったんだけど・・・。

QGISでxlsデータをインポートする方法を知らないのでArcGISでインポートする。
ArcGISと手順が違っていたので注意。
データ追加→「すべてのファイル」にしてxlsファイルを開く。
で、レイヤを右クリックしてファイル→データのエクスポートでshpファイルにする。
できあがり。シンボルやラベリングを調整。

と思ったら座標系がWGS84のままらしく、後々面倒なことに。
iGISに転送しても表示が変に。座標系の変換を試みるもうまくいかず、今に至る。
あとx,yが空白なのはまずいみたい。

LANDSAT画像のダウンロード→データセット作成

GoogleEarthの背景画像だと積雪期のがないよなーということでLANDSAT画像を見てみる。

Landsat Search and Download:
http://landsat.usgs.gov/Landsat_Search_and_Download.php


Glovisの方が検索は簡単。EarthExplorerの方が詳細に検索がかけられるかな？
でもほとんどの画像が3バンドまでのJPEGファイル。座標はついてるけど。
8バンドまであるGeotiffファイルはなかなか見つけられない。

GlovisならDL画面に行く前にわかる。変なマーク。

Level1Productでないと各Bandのダウンロードができない。

"Natural Color"JPEGファイルだとこんなの。

---

つぎにGeotiffファイルを複数バンドのデータセットにまとめる。
ArcGISでもできた記憶があるけど、EARDAS IMAGINEでもできると教えてもらったのでやってみる。
ちなみに学内ネットワークに接続しないとEARDASが起動しないのは本当だった。なにこのライセンス。

tiffファイルを開くときに多分ピラミッドを作るか聞いてくる。OKするしかない。
ピラミッドは大きいデータを拡大縮小するときにあらかじめ縮小したときの画像を作っておくことなんだって！

Raster→Spectral→Layer Stackで、重ねる.tiffファイルを選択→Add→重ねるtiffファイルを選択→Add→...を繰り返して適当な名前をつけておk

RGBのバンドを変えたくなったけど眠いから寝る。