ラスタ演算をGDALとOpenMPで高速に実行する方法

10mメッシュ標高DEMなどの高解像度のラスタデータを広域に計算しようとすると、とても時間がかかって待ってられません。
そこで、マルチコアCPU(Core2Duoとか)をフル活用して、短時間でラスタ演算できるGDALとOpenMPを使った並列プログラミングの方法を紹介します。

手順

1. Visual c++ 2008 Express Edition をインストール(プログラミング環境)
2. Windows SDK for Windows Server 2008 and .NET Framework 3.5をインストール(OpenMP環境)
3. GDALライブラリをインストール(ラスタ演算環境)
4. GDALとOpenMPを使ってプログラミング
5. 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;
}

5.OpenMPオプションでコンパイル、そして実行

OpenMPオプションとgdalのライブラリをリンクしてコンパイル

cl /EHsc /openmp -IC:\warmerda\bld\include C:\warmerda\bld\lib\gdal_i.lib dem.cpp