OSGeo4WのPythonでR、GRASS、GDALを連携させる方法 - 自然環境保全のための周辺技術

windows7 64bit PCで、OSGeo4Wを使ったpython、r、grass、gdalを連携させるための環境設定メモと、その他もろもろです。
ようするにQGISのプロセッシング(SEXTANTE)をQGISを起動させずにpython スクリプトでやりたいということです。

※この記事はまだ途中書きですが、忘れそうなので、とりあえずアップします。

RとPythonの連携

1. OSGeo4W 32bit版をインストール。64bit版だとGRASSが上手く動作しない場合があったので。http://trac.osgeo.org/osgeo4w/
2. OSGeo4Wからeasy_installとpipを利用できるように設定。https://trac.osgeo.org/osgeo4w/wiki/ExternalPythonPackages
3. R3.0.2をインストール。http://www.r-project.org/
4. OSGeo4WでRの環境変数を設定

set PATH=C:\Program Files\R\R-3.0.2\bin;%PATH%
set R_HOME=C:\Program Files\R\R-3.0.2

5. rpy2のバイナリをダウンロードしてインストール。http://www.lfd.uci.edu/~gohlke/pythonlibs/

easy_install rpy2-2.3.8.win32-py2.7.exe

6. pythonで以下のコードを打って動作確認

import rpy2.robjects as robjects
robjects.r("""library("MASS")""")
stepAIC = robjects.r("stepAIC")
robjects.r("""data("iris")""")
rslt = robjects.r("""lm(Sepal.Length~.,data=iris)""")
print(stepAIC(rslt))

http://www.slideshare.net/gepuro/rpython

ついでに

Rとは関係ないけど、そのうち使いそうな以下のパッケージもインストール。
ただし、OSGeo4Wで入るものとバージョンが違うものがあるので、そのあとQGISとかGRASSが正しく動くかどうかは？
（いまのところ大丈夫だけど。。。）

easy_install scipy-0.13.2.win32-py2.7.exe
easy_install pandas-0.13.0.win32-py2.7.exe
easy_install numpy-MKL-1.8.0.win32-py2.7.exe
easy_install matplotlib-1.3.1.win32-py2.7.exe
pip install ipython #これもeasy_installでいいとおもうけど、pipで入ったのでpipでいれた

参考

https://sites.google.com/site/aslugsguidetopython/data-analysis/pandas/calling-r-from-python
http://rpy.sourceforge.net/rpy2/doc-2.1/html/introduction.html#getting-started

補足

pyr2ではなくpyperというのもあるらしい
http://mia-0032.hatenablog.jp/entry/2013/08/30/000000

OSGeo4Wでのpipを使ったコンパイル＆インストール

単にpip install rpy2とかだと、コンパイラがないため"Unable to find vcvarsall.bat"のエラーがでる。
VC2008Express Editionを入れてコンパイルしても相性が悪い？せいかうまくいかない。
記事によるとMinGWのgccを使えばうまくいくらしい。けど試していない。

http://ultrainfinitum.blogspot.jp/2012/12/python-error-unable-to-find-vcvarsallbat.html
http://lists.faunalia.it/pipermail/animov/2013-August/001125.html

pythonのディストリビューションを利用する

OSGeo4Wのpythonにこだわらなければ、windows用のpython ディストリビューションを利用すればよいのかも。
https://www.enthought.com/downloads/
https://store.continuum.io/
ただし、rpy2は入ってないみたい。

GRASSとPythonの連携

OSGeo4Wで以下の環境変数を設定する。バッチファイルにしておくと便利かも

REM インストーラーバージョンの場合はOSGEO4W_ROOTは↓に置き換えてください。
REM set OSGEO4W_ROOT=C:\Program Files (x86)\QGIS Dufour\
set OSGEO4W_ROOT=C:\OSGeo4W
set GISBASE=%OSGEO4W_ROOT%\apps\grass\grass-6.4.3
set GISRC=%APPDATA%\GRASS6\grassrc6
set LD_LIBRARY_PATH=%GISBASE%\lib;%LD_LIBRARY_PATH%
set PATH=%GISBASE%\bin;%GISBASE%\etc;%GISBASE%\etc\python;%GISBASE%\lib;%GISBASE%\extralib;%OSGEO4W_ROOT%\apps\msys\bin;%PATH%
set GRASS_SH=%OSGEO4W_ROOT%\apps\msys\bin\sh.exe
set PYTHONPATH=%GISBASE%\etc\python

動作確認

import grass.script as grass
import grass.script.setup as gsetup
import os,sys

gisbase = os.environ['GISBASE']
gisdb = "C:/Users/hogehoge/grassdata"
location = "newLocation"
mapset = "testmapset"
gsetup.init(gisbase, gisdb, location, mapset)
  
r = grass.read_command("g.list", type='rast')
print r

参考

http://www.giscience.sakura.ne.jp/blog/?p=102
http://grasswiki.osgeo.org/wiki/GRASS_and_Python
http://grass.osgeo.org/programming6/pythonlib.html
http://grasswiki.osgeo.org/wiki/GRASS_Python_Scripting_Library
locationの作成もできるといいな↓を参考に
http://geoinformatics.fsv.cvut.cz/gwiki/On_scripting_GRASS_GIS:_Building_location-independent_command_line_tools

ポイントでラスタ値を抽出してグラフにする

その1　RとRGDALとPythonで

http://www.r-bloggers.com/extracting-raster-values-from-points-in-r-and-grass/
Rのコード

require(rgdal)
r<-readGDAL("dem5.tif")
r<-readGDAL("tci.tif")
p<-readOGR("test","sansyouuo")
proj4string(p) <- proj4string(r)
print(over(p,r))

R & Pythonのコード（途中、グラフ部分追加）

import rpy2.robjects as robjects
r = robjects.r
r("""require("rgdal")""")
rast = r.readGDAL("dem5.tif")
point = r.readOGR("test","sansyouuo")
r.proj4string(point)[0] = r.proj4string(rast)[0]
data = r.over(point,rast)
print data[0]

その2　PythonとGDALで

http://wiki.gis-lab.info/w/Extract_raster_values_at_vector_points_location_with_Python/GDAL

- そのままだとエラーなので書き直し

気が向いたらコードを修正

その3 GRASSとPythonで

やらない

ipython notebookを使ってみる

http://blog.godo-tys.jp/2013/07/10/3073/

GRASSとPythonで累積日射量を計算する

コードを書く

http://www.grassbook.org/examples_menu2nd.php
http://help.arcgis.com/ja/arcgisdesktop/10.0/help/index.html#/na/009z000000t9000000/
http://grass.osgeo.org/grass65/manuals/r.sun.html

感想など

QGISのSEXTANTEってすごいなーと思っていたのですが、つまりPythonの守備範囲がすごいということが分かってきました。
http://chezou.wordpress.com/2014/01/18/%E7%A7%91%E5%AD%A6%E8%A8%88%E7%AE%97%E3%81%AB%E3%81%8A%E3%81%91%E3%82%8B%E5%9D%87%E8%B3%AA%E5%8C%96%E3%80%81%E3%81%82%E3%82%8B%E3%81%84%E3%81%AF%E3%81%AA%E3%81%9Cpython%E3%81%8C%E7%9D%80%E5%AE%9F/