これまでに傾斜量の算出方法を紹介し，Pyhonプログラムも掲載しました．これらの記事で用いた標高データ（傾斜量計算の入力データ）は比較的狭い範囲でしたので，計算時間はPython実装でも数分程度で完了しました．しかし，土木分野の実務では，より広域な範囲に対して地形解析を行う場合もあるため，このような場合では，業務の進捗を考えると計算時間は非常に重要なファクタとなります．
そこで，傾斜量の計算時間について，今回はC実装（最適化あり）と簡単に比較してみようと思います．
1．システム環境
傾斜量計算のPython実装とC実装を比較する環境は以下の通りです．
- CPU：Intel(R) Core(TM) i7-9700K CPU @ 3.60GHz
- メモリ：64GB
- OS：Ubuntu 18.04.4 LTS
- Pythonのバージョン：Python 3.6.9
- GCCのバージョン：gcc version 7.5.0
2．条件
傾斜量の計算方法は，「傾斜量の算出方法（その２）」で紹介した戸田の計算方法としました．Pyhonのプログラムは，そこに掲載しています．そして，ほぼ同等のC実装を以下に掲載します．計算のボトルネックとなる部分は，傾斜量計算の2重ループになります．
入力となる標高データは，前回と同様，兵庫県が公開している1mピッチの標高データをGeoTIFFファイル（131MB）に変換したものです．C実装のコンパイル時にはO3オプションを付加しました．そして，計算時間の測定にはtimeコマンドを利用しました．したがって，計算時間はファイルI/Oと傾斜量計算の時間を足した時間となります．
- 戸田堅一郎，”安全な路網計画のための崩壊危険地ピンポイント抽出技術 -CS立体図を用いた崩壊危険地形判読技術の開発-”，長野県林業総合センター，32号，pp.1-16，2018年
#include <math.h>
#include <gdal_priv.h>
#include <cpl_conv.h>
int main(int argc, char *argv[]){
//===========================================
// GeoTIFFを開く
//===========================================
GDALAllRegister();
GDALDataset *pInTif = (GDALDataset *)GDALOpen(argv[1], GA_ReadOnly);
if(NULL == pInTif){ // openに失敗した場合
printf("Error : GDALOpen\n");
exit(1);
}
double dTransform[6];
pInTif->GetGeoTransform(dTransform);
int iSizeX = pInTif->GetRasterBand(1)->GetXSize();
int iSizeY = pInTif->GetRasterBand(1)->GetYSize();
double dNoData = pInTif->GetRasterBand(1)->GetNoDataValue();
// GeoTIFFのデータを取得
float **fElevMat = (float**)CPLMalloc(sizeof(float*) * iSizeY);
for(int i = 0 ; i < iSizeY ; ++i){
fElevMat[i] = (float*)CPLMalloc(sizeof(float) * iSizeX);
CPLErr Ret = pInTif->GetRasterBand(1)->RasterIO(GF_Read, 0, i, iSizeX, 1, fElevMat[i], iSizeX, 1, GDT_Float32, 0, 0);
}
//===========================================
// 傾斜量の計算
//===========================================
// 結果を格納する2次元配列
float **fRetMat = (float**)CPLMalloc(sizeof(float*) * iSizeY);
for(int i = 0 ; i < iSizeY ; ++i)
fRetMat[i] = (float*)CPLMalloc(sizeof(float) * iSizeX);
// 傾斜量の計算
float fDisttance = 1.;
for(int i = 1 ; i < iSizeY-1 ; ++i){
for(int j = 1 ; j < iSizeX-1 ; ++j){
if(fabs(fElevMat[i][j] - dNoData) <= __DBL_EPSILON__) // NoDataの場合
fRetMat[i][j] = fElevMat[i][j];
else {
float fDzDx = ((fElevMat[i-1][j+1] + 2 * fElevMat[i][j+1] + fElevMat[i+1][j+1]) * 4
- (fElevMat[i-1][j-1] + 2 * fElevMat[i][j-1] + fElevMat[i+1][j-1]) * 4)
/ (8 * fDisttance);
float fDzDy = ((fElevMat[i+1][j-1] + 2 * fElevMat[i+1][j] + fElevMat[i+1][j+1]) * 4
- (fElevMat[i-1][j-1] + 2 * fElevMat[i-1][j] + fElevMat[i-1][j+1]) * 4)
/ (8 * fDisttance);
float fTmp01 = sqrtf(fDzDx * fDzDx + fDzDy * fDzDy);
fRetMat[i][j] = atanf(fTmp01) * (180 / M_PI);
}
}
}
// 計算結果をGeoTIFFとして出力
GDALDriver *pDriverTiff = GetGDALDriverManager()->GetDriverByName("GTiff");
char **papszOptions = NULL;
papszOptions = CSLSetNameValue(papszOptions, "COMPRESS", "LZW");
GDALDataset *pOutTif = pDriverTiff->Create("Slop.tif", iSizeX, iSizeY, 1, GDT_Float32, papszOptions);
pOutTif->GetRasterBand(1)->SetNoDataValue(dNoData);
pOutTif->SetGeoTransform(dTransform);
pOutTif->SetProjection(pInTif->GetProjectionRef());
for(int i = 0 ; i < iSizeY ; ++i)
CPLErr Ret = pOutTif->GetRasterBand(1)->RasterIO(GF_Write, 0, i, iSizeX, 1, fRetMat[i], iSizeX, 1, GDT_Float32, 0, 0);
GDALClose(pInTif);
GDALClose(pOutTif);
return 0;
}
3．計測結果
以下が，Python実装とC実装の計算時間となります．
Python実装は少し手を加えることで高速化でき，C実装では最適化オプションを付けてコンパイルしたため，Python実装にはかなり不利な条件ではありますが，user比で約190倍もC実装の方が早くなりました．
GDALの特殊なメモリ管理は必要ではありますが，実務ではC実装を採用した方がいいと考えられます．
Python実装の計算時間
real 10m43.567s
user 10m42.611s
sys 0m0.656s
C実装の計算時間
real 0m3.583s
user 0m3.339s
sys 0m0.244s