2017年2月27日月曜日

PythonでGeoTIFFの地形開度を計算する(並列化済)

概要

DEMから地形開度(横山ほか, 1999)を求めるPythonスクリプトを作成した。

背景

「地形開度」は(横山ほか, 1999)で提案されている地形の表現方法である。地上開度と地下開度の2つからなり、それぞれ尾根線と谷線を抽出・強調する効果がある。

なお(地上開度-地下開度)÷2という演算を行うことで有名な赤い地図にも用いられる「尾根谷度」(千葉ほか, 2007)を生成することが可能であるが、こちらは某社の特許に触れる懸念があるので扱わない。

この地形開度の計算はSaga GISでMPI (Morphometrical Protection Index)という名前で可能であり、QGISからもツールボックスを介して実行することができる。ところがこの機能、非常に挙動不審でおかしな出力をすることもあり、使っていてとてもストレスが溜まるものだった。加えて進捗状況の表示がないので、クラッシュしているのかどうかの判別もつかない。そこで先述の文献を参考に、DEMから地上開度と地下開度の双方を一度に計算するPythonスクリプトを作成した。

手法

プログラムの概要

  • スクリプトはPython 2.7.11で動作を確認
  • 基本モジュールに加えて、numpy(数値演算)・gdal(GIS)が必要
  • (自分が使ったのは、ここで公開されている64bit向けPythonのライブラリ)
  • 並列化処理により、スクリプト言語のPythonながらSaga GISに近い処理速度を実現

ワークフロー

  1. gdalモジュールでGeoTiffをnumpyの配列として読み出す (class GeoTIFF)
  2. 各ノード(メッシュ)に対して地上開度・地下開度を計算 (search)
  3. 受け取った計算結果を配列に書き込む (receive)
  4. 計算結果をGeoTiffとして書き出す (GeoTIFF.gdal_save)

入出力データ

  • Pythonのgdalモジュールが扱えるラスタ形式のファイル(デフォルトはGeoTiff)
  • ... 基本的に平面直角座標系(UTMや日本平面直角座標系)である必要がある。
  • パラメータ(radius)。地形開度を計算する上での考慮範囲で、単位は座標系に依存(多くの場合メートル)
  • ... 小さすぎると効果が薄く、大きすぎても無駄な計算時間が増える。現実的にはグリッド幅×(10~20)程度。

免責と注意点

作者はスクリプトを実行した際に生じたいかなる損害も感知しない。また結果の正確性についても保証しない。例えば幾つか明らかにSaga GISと挙動が異なるケースがある事を確認しているが、深く検証していない(というか本スクリプトの方が品質が良いのはなぜ…)。

スクリプト

import os, numpy, multiprocessing, gdal
from math import atan, sqrt, pi
from time import time

### presets ###
os.chdir(r'E:\ ')
infile ='test1.tif' # Input file (Any raster data GDAL able to read)
nodata = 0
radius = 10.0
mpi_file1 = infile.split('.')[0]+'_mpi1.tif'
mpi_file2 = infile.split('.')[0]+'_mpi2.tif'
max_proc = 4 #multiprocessing.cpu_count() ## To all logical processsor cores

class GeoTIFF:
    __c__ = 0
    def nanize(i0,i1): self.Z[i0,i1] = numpy.nan
    def __init__(self, srcfile='tmp.tif', band=1, nodata=0):
        print "[Open] < %s" % srcfile,
        gdal.UseExceptions(); gdal.AllRegister()
        ds = gdal.Open(srcfile)
        band = ds.GetRasterBand(band)
        self.GT = ds.GetGeoTransform() #x0, dx, dxdy, y0, dydx, dy
        self.shape, self.nodata = (band.XSize,band.YSize), nodata
        self.x0, self.y0, self.dx, self.dy = self.GT[0], self.GT[3], self.GT[1], self.GT[5]        
        self.Z = band.ReadAsArray() ## Load as array
        [self.nanize(i[0],i[1]) for i in numpy.rot90(numpy.where(self.Z==nodata))] #Apply nodata
        print "(%d x %d)" % self.shape
    def gdal_save(self, outfile='tmp2.tif', driver='GTiff', band=1):
        if os.path.exists(outfile): os.remove(outfile)
        ds = gdal.GetDriverByName(driver).Create(outfile, self.shape[0],self.shape[1], band, gdal.GDT_Float32)
        ds.SetGeoTransform(self.GT)
        band = ds.GetRasterBand(band)
        band.SetNoDataValue(self.nodata)
        band.WriteArray(self.Z)
        print "[Saved] > %s" % outfile
    def counter(self, reset=False):
        if reset: self.__c__=reset
        else:
            self.__c__ += 1
            if(self.__c__%100000==0): print "%d%% ->" % (100*self.__c__/self.shape[0]/self.shape[1]),
            if(self.__c__%1000000==0): print "(%d / %d)" % (self.__c__, self.shape[0]*self.shape[1])

timestart = time()
## main () ##############################################################
G = GeoTIFF(srcfile=infile, band=1, nodata=nodata)
dx,dy,dxy = abs(G.dx), abs(G.dy), sqrt(G.dx**2+G.dy**2)
idx,idy,idxy = int(float(radius)/dx), int(float(radius)/dy), int(float(radius)/dxy)
G.Z0, G.Z1, G.Z2 = G.Z[:], numpy.zeros(G.Z.shape),numpy.zeros(G.Z.shape)

def search(ixiy):
    ix,iy = ixiy[0], ixiy[1]
    slopes = {}
    slopes['E'] = [atan((G.Z[iy,ix+i]-G.Z[iy,ix])/(dx*i)) for i in numpy.arange(idx)+1 if (ix+i < G.shape[0])]
    slopes['W'] = [atan((G.Z[iy,ix-i]-G.Z[iy,ix])/(dx*i)) for i in numpy.arange(idx)+1 if (ix-i >=0)]
    slopes['N'] = [atan((G.Z[iy+i,ix]-G.Z[iy,ix])/(dy*i)) for i in numpy.arange(idy)+1 if (iy+i < G.shape[1])]
    slopes['S'] = [atan((G.Z[iy-i,ix]-G.Z[iy,ix])/(dy*i)) for i in numpy.arange(idy)+1 if (iy-i >=0)]
    slopes['NE']= [atan((G.Z[iy+i,ix+i]-G.Z[iy,ix])/(dxy*i)) for i in numpy.arange(idxy)+1 if (ix+i < G.shape[0]) and (iy+i < G.shape[1])]
    slopes['NW']= [atan((G.Z[iy+i,ix-i]-G.Z[iy,ix])/(dxy*i)) for i in numpy.arange(idxy)+1 if (ix-i >=0) and (iy+i < G.shape[1])]
    slopes['SE']= [atan((G.Z[iy-i,ix+i]-G.Z[iy,ix])/(dxy*i)) for i in numpy.arange(idxy)+1 if (ix+i < G.shape[0]) and (iy-i >=0)]
    slopes['SW']= [atan((G.Z[iy-i,ix-i]-G.Z[iy,ix])/(dxy*i)) for i in numpy.arange(idxy)+1 if (ix-i >=0) and (iy-i >=0)]
    up = [2/pi - max(slopes[key]) for key in slopes.keys() if not slopes[key]==[]] #Morphological Protection Index (Subaerial)
    down=[2/pi + min(slopes[key]) for key in slopes.keys() if not slopes[key]==[]] #Morphological Protection Index (Subsurface)
    if up==[] or down==[]: return False
    return ix,iy,numpy.average(up), numpy.average(down)

def receive(results):
    G.counter()
    if not results: return None
    ix,iy, subaerial, subsurface = results
    G.Z1[iy,ix], G.Z2[iy,ix] = subaerial, subsurface
    

if __name__ == '__main__':
    ### Enumerate target nodes ###
    args = [(ix,iy) for ix in range(G.shape[0]) for iy in range(G.shape[1]) if numpy.isnan(G.Z[iy,ix])==False and (G.Z[iy,ix]<>G.nodata)]

    ###Single processing ###
    ### 4000 pixel/s with radius=10
    #Results = map(search, args); map(receive, Results)
    
    ### Multi-processing Method 1 ###
    ### 14000 pixel/s with radius=10 | Progress not displayed ###
    #Pool = multiprocessing.Pool(max_proc)
    #Results = Pool.map(search,args); map(receive, Results)
    #Pool.close()

    ### Multi-processing Method 2 ###
    ### 11000 pixel/s with radius=10 | Progress displayed ###
    qsize = 1000*max_proc  ## Size of queue sent processes (Usually no need to change)
    Pool = multiprocessing.Pool(max_proc)
    while len(args)>qsize:
        Results = Pool.map(search,args[:qsize]); map(receive, Results)
        del args[:qsize]
    Results = Pool.map(search,args); map(receive, Results)
    Pool.close()
    
#Save to GeoTIFF files
G.Z = G.Z1; G.gdal_save(mpi_file1) # Export subaerial MPI
G.Z = G.Z2; G.gdal_save(mpi_file2) # Export subsurface MPI
print "...OK"

##########################################################################
print "Finished in %d seconds!!" % (time()-timestart)

効果と検証

Pythonはスクリプト言語であるため、シングルスレッドで処理した場合はCで書かれたSaga GISの処理モジュールに比べて一桁ほど遅い。このスクリプトでは同期型の並列化処理(multiprocessing.Pool.map)を導入したことにより、それを埋め合わせるだけの大幅な高速化を果たした。

具体的には、Intel MKLを使ってコンパイルされた64bit版のNumpyにおいてIntel Core-i7 6700HQ (実コア数4・仮想コア数8・全コア負荷時3.1GHz)で実行した場合、考慮距離radius=10で11000ピクセル毎秒の速度(Multi-processing Method 2)が得られた。これは1000x1000のDEMが1分半で処理できることを意味する。なお1スレッドのSingle processingでは4000ピクセル毎秒、進捗表示のないMulti-processing Method 1ではオーバーヘッドが減るために更に高速化し、14000ピクセル毎秒に達した。


2016年11月29日火曜日

NZ・カイコウラ地震(2016) 画像まとめ

ニュージーランド南島北部の太平洋岸で2016/11/14に起きたカイコウラ*地震(Kaikoura Earthquake Mw=7.8)について。来月行く予定なのでその為のまとめ。
*カイコウラ=カイクラ=Kaikoura


地震諸元

USGSの第一報。その後M7.8に更新された。


2016年11月22日火曜日

パホイホイ溶岩とアア溶岩の違いが生まれる理由について(動画あり)

玄武岩質溶岩に特徴的な「パホイホイ溶岩」と「アア溶岩」の違いが生まれるプロセスについて、日本語で手頃な資料がなかったので作成してみた。



2016年11月11日金曜日

学生の自費フィールドワーク文化は何が問題か?


発端


自費フィールドワーク文化とは何か?

  • 野外調査や学会参加など教員なら研究費として処理するようなコストを、カリキュラムで課された活動の一部であるにもかかわらず自費で補填する文化。
  • 何十年も前から続いてきた伝統的文化。現在の教員たちもそうやって育ってきた。
  • 基盤経費が減らされているので、教員も使える金が無い(当事者でないのでよく知らない)

(11/13追記) 反響ツイートがあったので、Twitterのモーメント(まとめ機能)でまとめました。当事者の経験、安く抑えるコツ、教員側の努力など。

2016年8月13日土曜日

BlenderGISでDEMを3Dレンダリング

DEM(標高ラスタ)を高品質な3Dでインタラクティブに描画できるソフトウェアはあまり多くない。商用ソフトでは色々あり、今まで使った中ではSurferが手軽で高品質、また動画作成に関してはFledermausが良かった。しかしやはりOSS(オープンソース・ソフトウェア)で実現したいところである。そこで、3DモデリングソフトBlenderのアドオンであるBlenderGISを試してみることにした。

2016年4月17日日曜日

Time Managerプラグインを使ったQGISによる時系列アニメーションの作成

4/14の熊本地震(M6)に端を発する地震活動は、16日に別府島原地溝の東へ飛び火し以降広い範囲で地震が群発し始めた。これを「震源が東に伸びていっている」と言った人がいたので、Hinet自動処理震源をQGISでプロットすることでどんなものか検証することにした。
必要な物は、震源データ、ベースマップ、そしてそれらをアニメーションとしてプロットする環境である。30秒の長きにわたるGoogle検索の結果、Time ManagerというQGISプラグインを見つけたのでそれを用いることにした。

2016年3月20日日曜日

MSI GS40 レビュー

2016年2月、2kg以下のノートPCとしては最もハイスペックな製品を買った。
MSI GS40 6QEという製品で、ジャンルとしてはいわゆるゲーミングノートの扱いになる。お値段は気を失っていたので覚えていない。購入したのはオーストラリアモデルで、日本国内でARKが販売している上位モデル(6QE-039JP)と同等スペックである。(値段もそんな感じなので察してほしい)

手にした実機の構成は以下のとおり。なかなか凄まじいスペックである。
サンプル袋に囲まれるGS40(手前・14インチ)とSVZ1311(奥・12インチ)