j3259の日記： 地形 (terrain)

脳味噌ぶら～ん link集/GIS系というページを見つける。モデルのファイルフォーマットでもお世話になったけど、地形データもカバーしてる。この管理人の You＆I という人はなんかすごい人だ。

USGS（米地質調査所）の DEM という高度地図のファイルフォーマットを使って地形をレンダリングするのにGDALというのが使えそうだなと思ってます。Windows でのバイナリが vterrain.org で配布されてる。libpng13.dll も一緒にダウンロードする必要あり。

データには10-meter DEM files in Victoriaを使ってみる。

USGS DEM の仕様は概要と仕様を一応目を通した。10-meter DEM は幸い "north up" だったけど、目の粗い普通の DEM は地球の北と UTM（ユニバーサル横メルカトル）地図の北がずれるので注意。あと、-32767 というのは void もしくはエラー容疑の値。地面の真ん中だとフィルタかけて補正する必要があると思うけど、海面っぽい場合は 0。
チュートリアルに従って DEM を読み込んでみたコード

    HeightMap * DEMLoader::load(const string& a_fileName)
    {
        GDALDataset * dataset;
        GDALAllRegister();
        dataset = (GDALDataset *) GDALOpen(a_fileName.c_str(), GA_ReadOnly);
 
        cout << "Driver: "
            << dataset->GetDriver()->GetDescription() << "/"
            << dataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME) << endl;
        cout << "Size is "
            << dataset->GetRasterXSize() << "x"
            << dataset->GetRasterYSize() << "x"
            << dataset->GetRasterCount() << endl;
        if (dataset->GetProjectionRef() != NULL) {
            cout << "Projection is " << dataset->GetProjectionRef() << endl;
        } // if
        double geoTransform[6];
        if (dataset->GetGeoTransform(geoTransform) == CE_None) {
            cout << "Origin = ("
                << geoTransform[0] << ", " << geoTransform[3] << ")" << endl;
            cout << "Pixel Size = ("
                << geoTransform[1] << ", " << geoTransform[5] << ")" << endl;
            cout << "Incline = ("
                << geoTransform[2] << ", " << geoTransform[4] << ")" << endl;
        } // if
 
        GDALRasterBand  * band;
        band = dataset->GetRasterBand(1);
        int xSize = band->GetXSize();
        int ySize = band->GetYSize();
        cout << "Band Type: " << GDALGetDataTypeName(GDALGetRasterDataType(band)) << endl;
 
        HeightMap * retval = new HeightMap();
        retval->setSize(xSize, ySize);
 
        cout << "Reading... ";
        for (int i = 0; i < xSize; i++) {
            band->RasterIO(GF_Read, i, 0, 1, ySize,
                retval->getColumn(i), 1, ySize, GDT_Float32,
                0, 0);
        } // for i
        cout << "Done" << endl;
 
        for (int i = 0; i < 12; i++) {
            cout << "(";
            for (int j = 0; j < 100; j++) {
                cout << retval->getHeight(i, j) << ", ";
            } // for j
            cout << ")," << endl;
        } // for i
 
        GDALClose(dataset);
 
        return retval;
    }

戻り値の HeightMap というクラスは float の二次元配列をクラスにしただけの簡単なもの。
DEMLoader はファイル読み込みに専念して LOD は別クラスで実装する。

#ifndef HEIGHT_MAP_H_
#define HEIGHT_MAP_H_
 
#include <iostream>
#include <string>
#include <cmath>
 
using namespace std;
 
namespace eed3si9n
{
class HeightMap
{
public:
    HeightMap();
    ~HeightMap();
    void setSize(int a_xSize, int a_ySize);
    float * getColumn(int a_xOffset);
    float getHeight(int a_xOffset, int a_yOffset);
private:
    float ** m_heights;
    int m_xSize;
    int m_ySize;
};
 
inline HeightMap::HeightMap():
    m_heights(NULL),
    m_xSize(0),
    m_ySize(0)
{
}
 
inline HeightMap::~HeightMap()
{
    if (m_heights == NULL) {
        return;
    } // if
 
    for (int i = 0; i < m_xSize; i++) {
        float * buffer = m_heights[i];
        delete buffer;
    } // for i
 
    delete [] m_heights;
}
 
inline void HeightMap::setSize(int a_xSize, int a_ySize)
{
    m_xSize = a_xSize;
    m_ySize = a_ySize;
 
    m_heights = new float *[m_xSize];
    for (int i = 0; i < m_xSize; i++) {
        m_heights[i] = new float[m_ySize];
    } // for i
}
 
inline float * HeightMap::getColumn(int a_xOffset)
{
    return m_heights[a_xOffset];
}
 
inline float HeightMap::getHeight(int a_xOffset, int a_yOffset)
{
    return m_heights[a_xOffset][a_yOffset];
}
 
//-------------------------------------------------
}
 
#endif /*HEIGHT_MAP_H_*/