斜度と植生だけのときと一緒。
QGISでGoogleMapと重ねてみた。QGISは表示の変更をするとすぐ落ちるのが嫌。
なんかずれてる。
拡大・縮小してたら直った。
2012年1月6日金曜日
AvalancheStartingZone②(雪崩発生地点の算出)
斜度・植生・積雪深を合計。AvalancheStartingZoneという名前は最近使ってない。論文格上ではわかりにくいと思って。
斜度と植生だけのときと一緒。
QGISでGoogleMapと重ねてみた。QGISは表示の変更をするとすぐ落ちるのが嫌。
なんかずれてる。
拡大・縮小してたら直った。
斜度と植生だけのときと一緒。
QGISでGoogleMapと重ねてみた。QGISは表示の変更をするとすぐ落ちるのが嫌。
なんかずれてる。
拡大・縮小してたら直った。
2011年12月27日火曜日
[Python]lxmlのインストール
通常通りにインストールできなかったのでMacPortsを使ってのインストール。
lxmlをインストール - YAMAGUCHI::weblog:
http://d.hatena.ne.jp/ymotongpoo/20090203/1233674446
上記どおり。
lxmlをインストール - YAMAGUCHI::weblog:
http://d.hatena.ne.jp/ymotongpoo/20090203/1233674446
上記どおり。
MacPortsのインストール
SnowLeopardのときにXcodeをインストールしてLion似した人はXcodeの再インストール。今のところ4.2.1。MacAppStoreでンストールできない！と思ってたけど多分インストーラーのダウンロードで止まってた。
MacPortsのインストールはLion用のMacPorts2.0.3をインストールしようとすると先に1.7.1をインストールしてくれと言われるので1.7.1を探してきてダウンロード。拡張子が.tar.gzのものを解凍してターミナルでそのフォルダに移動。
$ ./configure
$ make
$ sudo make install
でMacPorts1.7.1インストール完了。
再びLion用のMacPorts2.0.3を起動すればインストールできるはず。
MacPortsのインストールはLion用のMacPorts2.0.3をインストールしようとすると先に1.7.1をインストールしてくれと言われるので1.7.1を探してきてダウンロード。拡張子が.tar.gzのものを解凍してターミナルでそのフォルダに移動。
$ ./configure
$ make
$ sudo make install
でMacPorts1.7.1インストール完了。
再びLion用のMacPorts2.0.3を起動すればインストールできるはず。
2011年12月4日日曜日
植生データ作成④
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
とかすればいいんですね。
[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
とかすればいいんですね。
2011年12月3日土曜日
植生データ作成③
csvファイルをPythonの辞書型にするモジュールなんてなかった。ので自分で作る。こんなに手間なら手動でフィールド演算した方が早いよねなんて言わない。
まずcsvファイルの区切りを:に変える。xlsxファイルをLibreOfficeでcsvファイルで保存する際に区切り文字をコロンにして保存。
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="bash"]
>>> testdic={9999:10,9953:4}
>>> testdic
{9953: 4, 9999: 10}
>>> testdic[9999]
10
[/sourcecode]
あとは一行にしてそれぞれカンマで区切りたいわけなんだけど。また明日。
まずcsvファイルの区切りを:に変える。xlsxファイルをLibreOfficeでcsvファイルで保存する際に区切り文字をコロンにして保存。
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="bash"]
>>> testdic={9999:10,9953:4}
>>> testdic
{9953: 4, 9999: 10}
>>> testdic[9999]
10
[/sourcecode]
あとは一行にしてそれぞれカンマで区切りたいわけなんだけど。また明日。
2011年11月30日水曜日
QGISでジオリファレンス
スキャンしたjpg画像などに座標をつけてGISで重ねられるようにすることをジオリファレンスと呼んでます。多分正確な定義は違います。ArcGISではやったことがあったのでQGISでもできるだろうってことでやってみる。
QGISでジオリファレンスをするにはGDALジオリファレンサープラグインを使います。
[PDF] QGIS によるジオリファレンス（幾何補正について）
http://kanagawa-hozen-igaku.cocolog-nifty.com/blog/files/QGIS.pdf
[caption id="attachment_2068" align="alignnone" width="519" caption="いただいた画像なのでモザイクをかけています"]
[/caption]
ジオリファレンサープラグインを入れたらアイコンをクリックして起動。「ラスタを開く」アイコンで重ねたい画像を選択、「ポイントの追加」アイコンをクリックして画像をクリックするとその場所の座標入力を求められます。
「マップキャンパスより追加」ボタンではQGIS本体に戻って座標をクリックして指定できるらしいのですがなぜかできませんでした。
なので手動で座標を入力していきます。
国内ならウォッ地図で座標を知ることができるのですが、国外だとそうもいかないのでGoogleMapを利用しました。まじGoogle先生。
お決まりのGoogleMapで右上の歯車アイコンをクリックして「マップ Labs」>"緯度経度マーカー"を"オンにする"にチェック、「保存」
マップに戻って座標を調べたいところで右クリック>「緯度経度マーカーを設置」
座標が表示されます。
この座標をジオリファレンサーに入力していきます。Y,Xの順になってるので注意。
4ヶ所以上入力します。終わったとこ。
「ジオリファレンスの開始」アイコンをクリックします。ダイアログが現れますので設定画面にいきます。
"変換タイプ"とかもうわかりません。とりあえず"出力ラスタ"に適当な名前をいれて"ターゲットSRS"を(EPSG:4326(WGS84))にして"実行された時にQGISにロードします"のチェックを入れて「OK」
暫く待つとジオリファレンスが完了します。
非常にわかりにくいですが、できてます。
次は道をshpファイル化(デジタイズ)。
あと結局等高線もDEMも2枚のデータにまたがってるからマージとモザイク？しないと。
QGISでジオリファレンスをするにはGDALジオリファレンサープラグインを使います。
[PDF] QGIS によるジオリファレンス（幾何補正について）
http://kanagawa-hozen-igaku.cocolog-nifty.com/blog/files/QGIS.pdf
[caption id="attachment_2068" align="alignnone" width="519" caption="いただいた画像なのでモザイクをかけています"]
[/caption]ジオリファレンサープラグインを入れたらアイコンをクリックして起動。「ラスタを開く」アイコンで重ねたい画像を選択、「ポイントの追加」アイコンをクリックして画像をクリックするとその場所の座標入力を求められます。
「マップキャンパスより追加」ボタンではQGIS本体に戻って座標をクリックして指定できるらしいのですがなぜかできませんでした。
なので手動で座標を入力していきます。
国内ならウォッ地図で座標を知ることができるのですが、国外だとそうもいかないのでGoogleMapを利用しました。まじGoogle先生。
お決まりのGoogleMapで右上の歯車アイコンをクリックして「マップ Labs」>"緯度経度マーカー"を"オンにする"にチェック、「保存」
マップに戻って座標を調べたいところで右クリック>「緯度経度マーカーを設置」
座標が表示されます。
この座標をジオリファレンサーに入力していきます。Y,Xの順になってるので注意。
4ヶ所以上入力します。終わったとこ。
「ジオリファレンスの開始」アイコンをクリックします。ダイアログが現れますので設定画面にいきます。
"変換タイプ"とかもうわかりません。とりあえず"出力ラスタ"に適当な名前をいれて"ターゲットSRS"を(EPSG:4326(WGS84))にして"実行された時にQGISにロードします"のチェックを入れて「OK」
暫く待つとジオリファレンスが完了します。
非常にわかりにくいですが、できてます。
次は道をshpファイル化(デジタイズ)。
あと結局等高線もDEMも2枚のデータにまたがってるからマージとモザイク？しないと。
2011年11月23日水曜日
ASTER GDEMからGRASSフォーマットのDEMを作る①
どうも数値地図10mメッシュ標高データだとTitan2Dが動かないので他のDEMデータでやってみることにする。
30m解像度だけど扱いやすいASTER GDEMをダウンロードしてきてQGISで座標系変換(←今回ここまで)、GRASSにインポートしてDEMを作る。
まず、ASTER GDEMのダウンロード。こちらから。日本語もあるから問題ないはず。ログインしたら左メニュー「検索」からです。一時期ダウンロードできなくなってたけど直ってました。
解凍して所定のフォルダ(DETA5)に入れる。
次にQGISで座標系をWGS84からWGS84UTMzone54Nに変換する。もしかしてWGS84のままでもTitan2D動く？
_demの方に雪崩発生地点を重ねたとこ。どっちもWGS84です。
座標系を変換します。
メインメニュー「ラスタ」>「ワープ(再投影)」
"出力ファイル"は元の_demとおなじtif形式がいいんじゃないかな。"ソースSRS"がEPSG:4326(WGS84)になっていることを確認、"ターゲットSRS"をEPSG:32654(WGS84UTMzone54N)に「選択」して「OK」
じょーずにできましたー。
ワープが実行できない人はこちらのエントリーを参照するといいかも。
30m解像度だけど扱いやすいASTER GDEMをダウンロードしてきてQGISで座標系変換(←今回ここまで)、GRASSにインポートしてDEMを作る。
まず、ASTER GDEMのダウンロード。こちらから。日本語もあるから問題ないはず。ログインしたら左メニュー「検索」からです。一時期ダウンロードできなくなってたけど直ってました。
解凍して所定のフォルダ(DETA5)に入れる。
次にQGISで座標系をWGS84からWGS84UTMzone54Nに変換する。もしかしてWGS84のままでもTitan2D動く？
_demの方に雪崩発生地点を重ねたとこ。どっちもWGS84です。
座標系を変換します。
メインメニュー「ラスタ」>「ワープ(再投影)」
"出力ファイル"は元の_demとおなじtif形式がいいんじゃないかな。"ソースSRS"がEPSG:4326(WGS84)になっていることを確認、"ターゲットSRS"をEPSG:32654(WGS84UTMzone54N)に「選択」して「OK」
じょーずにできましたー。
ワープが実行できない人はこちらのエントリーを参照するといいかも。
2011年11月16日水曜日
[GRASS]ロケーション、マップセットの作成、ラスタデータのインポート
じゃあQGISからGRASSに変換するとこにデータが抜けたんだ。(そういえばGRASSツールボックスを使ったときに警告マークがでていた)と思ってGRASSオンリーでtifファイル(DEM)をインポートする。
GRASSを起動。スタート画面から"GIS Data Directory"を設定して(GRASS)、「Location wizard」をクリック。
"Project Location"に新しいロケーションの名前(AvalancheProject2)を入力して「次へ」
上から2番目の"Select EPSG code of spatial reference system"にチェックを入れて「次へ」
マップの"EPSG code"を選択(今回は32654 WGS84 UTM zone 54N)。「次へ」
確認して「完了」
次にマップセットの作成。スタート画面から「Create mapset」をクリック。
マップセットの名前を入力(tanigawa_dem)。「OK」
確認して「Start GRASS」
続いてラスタデータのインポート。(今回はtifファイル)
メインメニュー>「File」>「Import raster data」>「Common import formats」
"Source name"の"File:"で「Browse」をクリックしてインポートしたいファイルを選択。"Format:"が正しいか確認して「Import」
いいですね。
GRASSを起動。スタート画面から"GIS Data Directory"を設定して(GRASS)、「Location wizard」をクリック。
"Project Location"に新しいロケーションの名前(AvalancheProject2)を入力して「次へ」
上から2番目の"Select EPSG code of spatial reference system"にチェックを入れて「次へ」
マップの"EPSG code"を選択(今回は32654 WGS84 UTM zone 54N)。「次へ」
確認して「完了」
次にマップセットの作成。スタート画面から「Create mapset」をクリック。
マップセットの名前を入力(tanigawa_dem)。「OK」
確認して「Start GRASS」
続いてラスタデータのインポート。(今回はtifファイル)
メインメニュー>「File」>「Import raster data」>「Common import formats」
"Source name"の"File:"で「Browse」をクリックしてインポートしたいファイルを選択。"Format:"が正しいか確認して「Import」
いいですね。
(未完)VHubからTITAN2Dを動かす②
※この投稿は最後までうまくいってないのであまり参考にしないでください。
Titan2Dの使い方を教わったので谷川のDEMで試してみたのですが…
(24467.2,85070.4)で何か問題があるそうです。(あとから考えたらこの座標もおかしかった)
GRASSで見てみても特に問題はありません。
--
QGISで斜度を細分類したデータを眺めていたらデータの抜けがありました。
考えられること：
DEMを作る際、周りの余分なところを省いたときにNoDataにしてしまった？
Titan2Dの使い方を教わったので谷川のDEMで試してみたのですが…
(24467.2,85070.4)で何か問題があるそうです。(あとから考えたらこの座標もおかしかった)
GRASSで見てみても特に問題はありません。
--
QGISで斜度を細分類したデータを眺めていたらデータの抜けがありました。
考えられること：
DEMを作る際、周りの余分なところを省いたときにNoDataにしてしまった？
2011年11月12日土曜日
Google App Engineを使ってみる。
Google App Engine - Google Code
http://code.google.com/intl/ja/appengine/
右の文中の"登録"から「Create Application」をクリック、"Country and Carrier"は"Japan"、キャリアを選択して、@以前のアドレスを入力。携帯に認証コードが送られてくるので入力。携帯からGmailに転送設定をしてても転送されてこないみたい。
↑必要ないかも
http://code.google.com/intl/ja/appengine/downloads.htmlからMacOSX用のGoogleAppEngineLauncher-1.6.0.dmgをダウンロード。
GoogleAppEnginLauncherをアプリケーションフォルダに入れる。マウントしたままだとだめ。
ダブルクリック起動するといろいろダイアログが出るけどOKで大丈夫。
あとはマニュアル通りです。ローカルでファイル作成すればｲｲﾝﾀﾞﾖｰ
あとうまくいかなかったらコードはコピペでおｋ
Hello, World! - Google App Engine - Google Code
http://code.google.com/intl/ja/appengine/docs/python/gettingstarted/helloworld.html
最後まで行くとハロワと表示されたテキストがダウンロードされるはず。おぺん。
http://code.google.com/intl/ja/appengine/
右の文中の"登録"から「Create Application」をクリック、"Country and Carrier"は"Japan"、キャリアを選択して、@以前のアドレスを入力。携帯に認証コードが送られてくるので入力。携帯からGmailに転送設定をしてても転送されてこないみたい。
↑必要ないかも
http://code.google.com/intl/ja/appengine/downloads.htmlからMacOSX用のGoogleAppEngineLauncher-1.6.0.dmgをダウンロード。
GoogleAppEnginLauncherをアプリケーションフォルダに入れる。マウントしたままだとだめ。
ダブルクリック起動するといろいろダイアログが出るけどOKで大丈夫。
あとはマニュアル通りです。ローカルでファイル作成すればｲｲﾝﾀﾞﾖｰ
あとうまくいかなかったらコードはコピペでおｋ
Hello, World! - Google App Engine - Google Code
http://code.google.com/intl/ja/appengine/docs/python/gettingstarted/helloworld.html
最後まで行くとハロワと表示されたテキストがダウンロードされるはず。おぺん。
2011年11月2日水曜日
傾斜データ作成
分類基準を変えたので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となって描画されていない。
「データの追加」から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となって描画されていない。
2011年10月27日木曜日
MacにGRASSをインストールする
前回作ったDEMをGRASSで表示させようとしても表示されず、なんでだろうなーと思っていたら結局必要なものがインストールされてないだけでした。QGISをインストールすればGRASSも使えるようになると思っている人は要注意。自分です。
http://www.kyngchaos.com/software/frameworksの上からGDAL CompleteまでをOSに応じてダウンロードしてインストール。具体的にはFreeType framework v2.4.6-1、cairo framework v1.10.2-3a、GSL framework v1.15-1、Tcl/Tk Aqua 8.5.8-2、GDAL 1.8 Completeです。本当は必要ないのもあるかもしれません。インストールはダウンロードしたdmgファイルを解凍して中のpkgファイルをダブルクリックして指示に従うだけです。
そして最後にhttp://www.kyngchaos.com/software/grassからGRASS本体をダウンロード。同じようにしてインストール。
「アプリケーション」から「GRASS-6.4」を起動します。
暫くするとGRASSが起動しますので「Browse」をクリックしてロケーションの一つ上のフォルダを選択、ロケーションが表示されるので、"Accessible mapsets"で「PERMANENT」ではないマップセットを選択して「Start GRASS」。
ウィンドウが２つ表示されると思います。
GRASS GIS Layer Managerの左から2番目のアイコンをクリック。
"Load map layers into workspace"をクリック。
追加したいレイヤにチェックが入ってることを確認して「Load」
追加されたレイヤの左のチェックボックスにチェックを入れます。
もう一つのウィンドウにレイヤが表示されました。
ちゃんとQGISからGRASSに変換できてましたね。
http://www.kyngchaos.com/software/frameworksの上からGDAL CompleteまでをOSに応じてダウンロードしてインストール。具体的にはFreeType framework v2.4.6-1、cairo framework v1.10.2-3a、GSL framework v1.15-1、Tcl/Tk Aqua 8.5.8-2、GDAL 1.8 Completeです。本当は必要ないのもあるかもしれません。インストールはダウンロードしたdmgファイルを解凍して中のpkgファイルをダブルクリックして指示に従うだけです。
そして最後にhttp://www.kyngchaos.com/software/grassからGRASS本体をダウンロード。同じようにしてインストール。
「アプリケーション」から「GRASS-6.4」を起動します。
暫くするとGRASSが起動しますので「Browse」をクリックしてロケーションの一つ上のフォルダを選択、ロケーションが表示されるので、"Accessible mapsets"で「PERMANENT」ではないマップセットを選択して「Start GRASS」。
GRASS GIS Layer Managerの左から2番目のアイコンをクリック。
"Load map layers into workspace"をクリック。
追加したいレイヤにチェックが入ってることを確認して「Load」
追加されたレイヤの左のチェックボックスにチェックを入れます。
もう一つのウィンドウにレイヤが表示されました。
ちゃんとQGISからGRASSに変換できてましたね。
登録:
投稿 (Atom)