2011年12月27日火曜日

[Python]lxmlのインストール

通常通りにインストールできなかったのでMacPortsを使ってのインストール。
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を起動すればインストールできるはず。

2011年12月21日水曜日

VHubからTitan2D

VHubでエラーの質問をしたら、GRASSデータが圧縮されているのかも。という回答をいただきました。

ということでGRASSで

$ r.compress -u < map name >

を実行。



VHubにアップロードしなおして再計算。



成功?

TitanViewerでみてみる。/vhub/ユーザー名/titan2dで作成したフォルダ名/結果のタイムスタンプ を開いて、Standard Titan Viewだったかな?でRun Tool。そしたらResultのドロップダウンからMap of maximum flow depthを選んで地図を表示。



うーん。いじってたら地図も出なくなってしまった。

AvalancheStartingZone①―植生と斜度の合計





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

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

ラスタ変換。

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

結果はssとってない!

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
とかすればいいんですね。

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]

あとは一行にしてそれぞれカンマで区切りたいわけなんだけど。また明日。

2011年12月2日金曜日

植生データ作成②

ArcGIS10でフィールド演算をPythonコードでできるようになったというので、そのための下ごしらえ。

ダウンロードした植生フォルダにどの番号が何の植生かが書かれているcsvファイルがあるので、これをLibreOfficeで開いて.xlsxで保存(←あとあと考えたら.csvのままで良かった)。辞書で調べながら点数付けをしていく。



余計な列を削除して.csvで保存。



vimで開いてみる。あら。



余計な列をさらに削除して再度vimで開く。



なぜか;になってましたがこれはすぐに直せるでしょう。

次は.csvをpythonの辞書型にしようと思う。

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枚のデータにまたがってるからマージとモザイク?しないと。

2011年11月28日月曜日

植生データ作成①

新しい評価基準に基づき、植生ファイルを修正。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次メッシュ植生データのリンクから。

2011年11月27日日曜日

twitterBOTの作成③

タイムゾーンを東京(GMT+09:00)に直す。

datetime.datetime.now()はタイムゾーンを引数にとれるみたい。でも引数がちょっと面倒。

以下を参考にさせていただきました。

Pythonでタイムゾーンを扱うメモ « taichino.com
http://taichino.com/programming/1876

タイムゾーンが設定できたところでdatetime.datetime(2012,1,21) - datetime.datetime.now(JST())をすると今度はnaiveがどーたらこーたら。

今度は以下を参考にdatetime.datetime(2012,1,21)の引数にtzinfo=JST()を加える。

あるエンジニアの不定期勉強記録
http://abalone.ununu.org/archives/705



できました。この画像からはよくわかりませんが…。

現在のsotsuprobot.py

[sourcecode language="python"]
# -*- coding: utf-8 -*-
import twitter
import datetime

api = twitter.Api(
cache = None,
consumer_key = '****',
consumer_secret = '****',
access_token_key = '****',
access_token_secret = '****')

# タイムゾーンの設定
class JST(datetime.tzinfo):
def utcoffset(self, dt):
return datetime.timedelta(hours=9)
def dst(self, dt):
return datetime.timedelta(0)
def tzname(self, dt):
return 'JST'

# 残り時間の算出
teisyutsubi = datetime.datetime(2012,1,21,23,59,tzinfo=JST())
ima = datetime.datetime.now(JST())
nokorijikan = teisyutsubi - ima

# tweet
text = u"現在時刻"+ str(ima) + u" 卒プロ成果登録締め切りまで あと" + str(nokorijikan.days) + u"日 #sfc_sotsupro"

api.PostUpdate(status=text)
[/sourcecode]

cron.yaml

[sourcecode language="text"]
cron:
- description: tweet
url: /sotsuprobot
schedule: every 3 hours
timezone: Asia/Tokyo
[/sourcecode]

[twitter-follow screen_name='sotsupro2012' show_count='no']

twitterBOTの作成②

今日はBOT用アカウントを作成してアプリを刺してGoogleAppEngineで自動的に実行させる。

まずBOT用アカウントを作成。これは一般登録と同じ。同じメールアドレスは使えないんだね。昨年のBOT名を参考にしたら卒論じゃなくて卒プロでした。

次にBOT用アカウントにログインした状態で前回oauth2認証に使ったPythonを実行する。すると別のAccess tokenとAccess token secretが取得できる。

GoogleAppEngineに載せる。

以下を参考にさせて頂きました。

ほげおメモ: Python Twitter API python-twitter の使い方 (Part4): Google App Engine
http://blog.hogeo.jp/2011/06/python-twitter-api-python-twitter-part4.html

twitter.pyやhttplib2/やoauth2/は/Library/Python/2.7/site-packagesからコピーしたりPython-twitterをインストールする際にダウンロードしたフォルダからコピーしたり。

cron.yamlはとりあえず2分毎にツイートするように設定。これらを同じフォルダに格納。

GoogleAppEngineLauncherを起動してメインメニュー>「File」>「Add Existing Application」でフォルダを追加。あとhttps://appengine.google.com/で保存?するアプリケーション?を作成。「Deploy」



いろいろあったけどとりあえずこれで2分おきに発言するBOTができました。

Deploy中にいろいろいじるとDeployされなくなってしまうので、そのときはDeployしてる上のフォルダにいってrollbackする。

アプリをappcfgでrollbackせよ - CaMiL
http://d.hatena.ne.jp/beso/20110521/rollback

appcfg.py rollback /フォルダ名 でよかったはず。

結果がこちら。



タイムゾーンがおかしくなってる…。GoogleAppEngineは米国のサービスだからか。

つづく

2011年11月26日土曜日

twitterBOTの作成

そうだ、卒業論文提出日までの日数を自動で発言するtwitterBOTをPythonで作成すれば勉強になるんじゃないかと思って。

Pythonでtwitterを扱うPython-twitterのインストールは以下。

Python Twitter のインストール
http://www.yukun.info/blog/2011/03/python-twitter-install.html

コードの作成はほとんど以下を参考にさせていただきました。

ほげおメモ: How to make Twitter Bot by Python (Part 1) python-twitter
http://blog.hogeo.jp/2011/05/how-to-make-twitter-bot-python-twitter.html

part1~5まであります。part1はなぜか英語なのでTwitter API ポケットリファレンス (POCKET REFERENCE)も参考にしました。

[sourcecode language="python"]
# -*- coding: utf-8 -*-
import twitter
import datetime

api = twitter.Api(
consumer_key = '****',
consumer_secret = '****',
access_token_key = '****',
access_token_secret = '****')

teisyutsubi = datetime.date(2012,1,21)
kyou = datetime.date.today()
nokorijikan = teisyutsubi - kyou
text = u"卒業論文提出日まで あと" + str(nokorijikan.days) + u"日"

api.PostUpdate(status=text)
[/sourcecode]

2011年11月24日木曜日

QGISでDEMのカラーマップを変更する。陰影図を作る。

パタゴニアの地図をもっとわかりやすくする。

プラグイン — GeoPacific.org
http://www.geopacific.org/opensourcegis/gcngisbook/GCN_book/4ed89332/4ed89332b/qgis_plugin

 

カラーマップの変更→1-Band Raster Color Table V1.x

一度作ったカラーマップを変更するのはカラーマップのデータを削除しなければエラーがでる。

陰影図→relief shader

使い方は上記URLを参照。

 



ほとんど氷河なんだなーと。氷河も入れたい。どうしたらいいんだろう。眠い。

2011年11月23日水曜日

QGISでDEMから等高線データを作る

卒業旅行でパタゴニア(特にフィッツロイ&セロトーレ)に行きたいんだけど地図が手に入らなーいということで自作する。まず等高線を作ろう。

まずGoogleMapからフィッツロイの座標を調べる。なんでフィッツロイなのかというとセロトーレより有名なので。

座標を調べるツールは右上歯車アイコンから「マップLabs」>"経緯度ツールチップ"で"オンにする"を選択して「保存」。これでShiftキーを押すことでだいたいの座標を知ることができます。





こんなかんじ。

DEMはASTER GDEMから。解凍してQGISで表示。心配なので2枚とってきてしまった。ArcGISでいうモザイクをしてもいいかもそれません。



等高線を作るにはメインメニュー「ラスタ」>「等高線」。やってることは過去記事と同じです。

"入力ファイル"に元のDEM、"等高線の出力ディレクトリ"には適当な名前、"等高線の間隔"は今回20.00(m)にしました。残り2つのチェックボックスにチェックして「OK」。



しばらくかかりますが等高線が作られます。それぞれやっても等高線は違わないみたいですね。



元のDEMはもう洋ナシなのでリストから削除。



GoogleMapと重ねたいので右下から座標系をGoogle Mercator(EPSG:900913)にして'オンザフライ'CRS変換を有効にする。

メインメニュー「プラグイン」>「OpenLayers plugin」>「OpenLayers Overview」で表示しているところのGoogleMapが表示されるので、がんばって対象地域を探して+印のボタンをクリック。



今日はここまで。

次はDEMから標高による色分けをしようか陰影図もつくろうか。でもトレッキングルートを入れないとなあ。山頂のポイントデータもほしい。

ASTER GDEMからGRASSフォーマットのDEMを作る②

エントリー2つに分けたけどやること前のと同じだったw

ワープしたASTER GDEMのtifファイルをGRASSにインポートする。



今回ロケーション、マップセットはこのようにしました。もちろん座標系はEPSG:32654

インポートの手順は前回と同じ。



Titan2Dは精神的に疲れるから次回にしよう。