点群データを格子データにするもう一つの方法：GRASS r.in.xyz #Python3

はじめに

昨日Pythonで点群データを格子データにする手法について紹介した。
本記事はそこで言及した「もう一つの方法」、GRASS GISのr.in.xyzについて。

実はGRASS GISは殆ど使ったことがないのだが、プロジェクトをセットアップしないと立ち上がらないなど色々とっつきにくい雰囲気を漂わせている(※個人の印象です)。
しかしマニュアルを見てもらえれば分かるように、r.in.xyzに限らずともかなり魅力的なツールが揃っていることがわかる。しかしこれらは外部からはアクセスできない。

そこでお行儀が悪いが、以下のようなアプローチを取ってGRASSのGUIを起動せずに中のツールを使う。
1. GRASSに捨てプロジェクトを作らせる
2. GRASSに再帰的にスクリプトを読み込ませてGRASS内のPythonでr.in.xyzを走らせる

コード

まずr.in.xyzのための設定と、GRASSの関連パスを指定する。今回はQGISに内包されているGRASSを使うので、後者は全てQGISのディレクトリ内にある。
捨てプロジェクトはユーザフォルダの中にGRASSのデフォルトらしき名前で作る。場所を変えたら動かなくなったので触らない方が良さそう。

GIS_Grass-r-in-xyz.py

import os
### Configuration for r.in.xyz ###
selfile = os.path.abspath(__file__) ## Path to this file
xyzfile = r'G:\GIS\SampleData_Mariana\Mariana_GoodQ.xyz' ## File to load
outfile = r'G:\GIS\SampleData_Mariana\test.tif' ## File to export
bound = [13, 10.8, 140, 145] ## Grid extent: N-S-W-E
dx = 0.001 ## Grid resolution
epsg = 4326 ## 4326:WGS84

### Configuration for GRASS GIS ###
grass7x = r'C:\Program Files\QGIS 2.18\bin\grass72.bat' ## GRASS executable
grassdir= r'C:\Program Files\QGIS 2.18\apps\grass\grass-7.2.0'
gisdb = os.path.join(os.environ['HOME'], 'grassdata') ##Do not change
locname, mapset = 'newLocation', 'PERMANENT' ##Do not change

続いてGRASS内部のPythonではgrass.scriptモジュールが読み込めることを利用して、try~節で分岐させるというトリックを使う。

読み込めない場合(GRASS外の場合)、grass7x.batを呼び出してプロジェクトを生成する。GUIを使わずにプロジェクトを作成する方法についてはWorking with GRASS without starting it explicitlyで説明されている。
プロジェクトをセットアップ後、コマンドプロンプトを立ち上げgrass7x.bat --execute pythonでこのPythonスクリプト自身を読ませる。

読み込める場合(grass7x.bat --execute pythonで実行された場合)、
1. r.regionで範囲と解像度を指定
2. r.in.xyzを実行
3. r.out.gdalでGeoTIFFとしてエクスポート
という手順を踏む。

GIS_Grass-r-in-xyz.py

try:
    import grass.script as gscript
except ModuleNotFoundError:
    #os.environ['GRASS_MESSAGE_FORMAT'] = 'plain'
    os.environ['GISDBASE'], os.environ['LOCATION_NAME'], os.environ['MAPSET'] = gisdb, locname, mapset
    if not os.path.exists(gisdb): ## Create a new one if grassdata is not present
        cmd = '"%s" -c epsg:%s -e %s' % (grass7x, epsg, os.path.join(gisdb,locname))
        print(cmd); os.system(cmd)
    ## Using exec | See https://grass.osgeo.org/grass73/manuals/grass7.html#exec-interface-example
    cmd = 'cmd /c "%s" --exec python %s' % (grass7x, selfile)
    print(cmd); os.system(cmd)
else:
    def main():
        gscript.use_temp_region()
        ## See https://grass.osgeo.org/grass72/manuals/raster.html
        gscript.run_command('g.region', n=bound[0], s=bound[1], w=bound[2], e=bound[3], res=dx, flags='p')
        gscript.run_command('r.in.xyz', input=xyzfile, output='tmp', separator='tab', method='mean', overwrite=True)
        gscript.run_command('r.out.gdal', input='tmp', output=outfile, format='GTiff', overwrite=True)
    if __name__ == '__main__': main()

実行するとこんな感じでちゃんと動いた。 (Python 3.6.1 x64 / QGIS 2.18.7 / GRASS 72 on Windows 10)

なお相変わらず進捗表示はあてにならないが、できるだけコードを減らしたい場合はこちらの方が良いかもしれない。
ただし大変お行儀が悪いアプローチなのでずっと使えるかどうかは不明。