10mメッシュ標高DEMなどの高解像度のラスタデータを広域に計算しようとすると、とても時間がかかって待ってられません。
そこで、マルチコアCPU(Core2Duoとか)をフル活用して、短時間でラスタ演算できるGDALとOpenMPを使った並列プログラミングの方法を紹介します。
手順
- Visual c++ 2008 Express Edition をインストール(プログラミング環境)
- Windows SDK for Windows Server 2008 and .NET Framework 3.5をインストール(OpenMP環境)
- GDALライブラリをインストール(ラスタ演算環境)
- GDALとOpenMPを使ってプログラミング
- OpenMPオプションでコンパイル、そして実行
1.Visual c++ 2008 Express Edition をインストール(プログラミング環境)
↓こちらからダウンロード＆インストール
http://www.microsoft.com/japan/msdn/vstudio/Express/
2.Windows SDK for Windows Server 2008 and .NET Framework 3.5をインストール(OpenMP環境)
↓こちらからダウンロード＆インストール
http://www.microsoft.com/downloads/details.aspx?FamilyId=E6E1C3DF-A74F-4207-8586-711EBE331CDC&displaylang=en
ドキュメントとかはインストールしなくてもいいです。
3.GDALライブラリをインストール(ラスタ演算環境)
↓こちらから最新のソースをダウンロード
http://trac.osgeo.org/gdal/wiki/DownloadSource
Set PATH=C:\Program Files\Microsoft Visual Studio 9.0\VC\bin;C:\Program Files\Microsoft Visual Studio 9.0\Common7\IDE;C:\Program Files\Microsoft SDKs\Windows\v6.0A\Bin;%PATH%
Set INCLUDE=C:\Program Files\Microsoft Visual Studio 9.0\VC\INCLUDE;C:\Program Files\Microsoft SDKs\Windows\v6.0A\Include;%INCLUDE%
Set LIB=C:\Program Files\Microsoft Visual Studio 9.0\VC\LIB;C:\Program Files\Microsoft SDKs\Windows\v6.0A\Lib;%LIB%
nmake.optをメモ帳で開いて、MSVC_VER= のところを MSVC_VER=1500に変更
コマンドプロンプトで、コンパイル＆インストール
nmake /f makefile.vc
nmake /f makefile.vc install
nmake /f makefile.vc devinstall
4.GDALとOpenMPを使ってプログラミング
gdalとopenmpを使ってプログラミング。詳細は、また今度。
とりあえず、傾斜を出すプログラムを張っておきます。
#include <stdio.h> #include <stdlib.h> #include <math.h> #include "gdal_priv.h" #include <omp.h> int main(int argc, char **argv) { if (argc < 4) { printf("###############################\n"); printf("USAGE:\n"); printf("slope_mp.exe input_dem.tif output.tif div\n"); printf("###############################\n"); exit(1); } GDALDataset *poDataset; const float degrees_to_radians = 3.14159 / 180.0; const float radians_to_degrees = 180.0 / 3.14159; double adfGeoTransform[6]; const char *pszFilename = argv[1]; const char *pszSlopeFilename = argv[2]; int divmax = atoi(argv[3]); int R = 10; GDALAllRegister(); poDataset = (GDALDataset *) GDALOpen(pszFilename, GA_ReadOnly); GDALRasterBand *poBand; poBand = poDataset->GetRasterBand(1); poDataset->GetGeoTransform(adfGeoTransform); const double cellsizeY = adfGeoTransform[5]; const double cellsizeX = adfGeoTransform[1]; const float nullValue = -9999; const int nXSize = poBand->GetXSize(); const int nYSize = poBand->GetYSize(); const int celln = ceil(R * 2 / cellsizeX); const int pcell = (celln + 1) / 2; GDALDriver *poDriver; poDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); GDALDataset *poSlopeDS; GDALRasterBand *poSlopeBand; char **papszOptions = NULL; papszOptions = CSLSetNameValue(papszOptions, "PROFILE", "GeoTIFF"); poSlopeDS = poDriver->Create(pszSlopeFilename, nXSize, nYSize, 1, GDT_Float32, papszOptions); poSlopeDS->SetGeoTransform(adfGeoTransform); poSlopeDS->SetProjection(poDataset->GetProjectionRef()); poSlopeBand = poSlopeDS->GetRasterBand(1); poSlopeBand->SetNoDataValue(-9999); int percent[11]; for (int i = 0; i <= 10; i++) { percent[i] = nXSize * nYSize * 0.1 * i; } int start[100], end[100]; for (int i = 0; i < divmax; i++) { start[i] = i * nYSize / divmax - celln; end[i] = (i + 1) * nYSize / divmax - celln; if (i == 0) start[i] = 0; } float win[10000]; float SlopeBuf[50000]; int p = 0; int strp = 0; for (int div = 0; div < divmax; div++) { float *bigwin = (float *) VSIMalloc3((end[div] - start[div] + celln), nXSize, sizeof(float)); if (bigwin == NULL) { printf("### Please gain div count. ###\n"); exit(1); } poBand->RasterIO(GF_Read, 0, start[div], nXSize, (end[div] - start[div] + celln), bigwin, nXSize, (end[div] - start[div] + celln), GDT_Float32, 0, 0); #pragma omp parallel for firstprivate(win,SlopeBuf) for (int i = start[div]; i < end[div]; i++) { for (int j = 0; j < nXSize; j++) { #pragma omp critical { if (p >= percent[strp]) { printf("%d%s", strp * 10, "..."); strp++; } p++; } if (i > nYSize - celln || j > nXSize - celln) { SlopeBuf[j] = nullValue; continue; } for (int ii = 0; ii < celln; ii++) { for (int jj = 0; jj < celln; jj++) { win[ii * celln + jj] = bigwin[((i - start[div] + ii) * nXSize + j) + jj]; } } float dx = ((win[celln * (pcell - 2) + pcell - 2] + win[celln * (pcell - 1) + pcell - 2] + win[celln * (pcell) + pcell - 2]) - (win[celln * (pcell - 2) + pcell] + win[celln * (pcell - 1) + pcell] + win[celln * (pcell) + pcell - 2])) / (6 * cellsizeX); float dy = ((win[celln * (pcell - 2) + pcell - 2] + win[celln * (pcell - 2) + pcell - 1] + win[celln * (pcell - 2) + pcell]) - (win[celln * (pcell) + pcell - 2] + win[celln * (pcell) + pcell - 1] + win[celln * (pcell) + pcell])) / (6 * cellsizeY); SlopeBuf[j] = atan(sqrt(dx * dx + dy * dy)) * radians_to_degrees; } poSlopeBand->RasterIO(GF_Write, pcell - 1, pcell - 1 + i, nXSize - (pcell - 1), 1, SlopeBuf, nXSize - (pcell - 1), 1, GDT_Float32, 0, 0); } VSIFree(bigwin); } delete poBand; delete poSlopeDS; return 0; }
