概要
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に近い処理速度を実現
ワークフロー
- gdalモジュールでGeoTiffをnumpyの配列として読み出す (class GeoTIFF)
- 各ノード(メッシュ)に対して地上開度・地下開度を計算 (search)
- 受け取った計算結果を配列に書き込む (receive)
- 計算結果を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(self,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ピクセル毎秒に達した。
0 件のコメント:
コメントを投稿