基于GDAL的栅格数据/遥感影像IO (非分块)

本文链接：https://blog.csdn.net/whuyao/article/details/51227144

本文介绍了一个用于简化GDAL使用的封装库，支持并行读取影像数据，并提供了投影变换功能。该库支持常见的栅格影像格式，利用OpenMP实现多线程处理。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

这个封装的目的主要是为了简化GDAL的读写文件的使用，支持投影变换，支持openMP。
需要的库：GDAL 2.x, Proj
库的网址：http://www.gdal.org
支持常见栅格影像

GDALRead.h

/************************************************************************/
/* This class is wrote for supercomputer to read image use parallel.
/* Not support Block read & write
/* But support multi-thread
/* Author: Y. Yao
/* E-mail: whuyao@foxmail.com
/* Version: v4.0
/************************************************************************/

#ifndef CLASS_GDAL_READ
#define CLASS_GDAL_READ

#include "gdal_priv.h"
#include "ogr_core.h"
#include "ogr_spatialref.h"

class CGDALRead
{
public:
    CGDALRead(void);
    ~CGDALRead(void);

public:
    bool loadFrom(const char* _filename);
    unsigned char* read(size_t _row, size_t _col, size_t _band);
    unsigned char* readL(size_t _row, size_t _col, size_t _band);   //extension read-data

    //get data via bilinear interpolation
    template<class TT> double linRead(double _row, double _col, size_t _band);

public:
    void close();
    bool isValid();

public:
    GDALDataset* poDataset();
    size_t rows();
    size_t cols();
    size_t bandnum();
    size_t datalength();
    double invalidValue();
    unsigned char* imgData();
    GDALDataType datatype();
    double* geotransform();
    char* projectionRef();
    size_t perPixelSize();

public:
    bool world2Pixel(double lat, double lon, double *pcol, double *prow);
    bool pixel2World(double *lat, double *lon, double col, double row);
    bool pixel2Ground(double col,double row,double* pX,double* pY);
    bool ground2Pixel(double X,double Y,double* pcol,double* prow);

protected:
    template<class TT> bool readData();

protected:
    GDALDataset* mpoDataset;    //=>
    size_t mnRows;                  //
    size_t mnCols;                  //
    size_t mnBands;             //
    unsigned char* mpData;      //=>
    GDALDataType mgDataType;    //
    size_t mnDatalength;            //=>
    double mpGeoTransform[6];   //
    char msProjectionRef[2048]; //
    char msFilename[2048];      //
    double mdInvalidValue;
    size_t mnPerPixSize;            //=>

public:
    OGRSpatialReferenceH srcSR;
    OGRSpatialReferenceH latLongSR;
    OGRCoordinateTransformationH poTransform;       //pixel->world
    OGRCoordinateTransformationH poTransformT;      //world->pixel
};

#endif

GDALRead.cpp

#include "GDALRead.h"
#include "ogrsf_frmts.h"
#include <iostream>
using namespace std;

CGDALRead::CGDALRead(void)
{
    mpoDataset = NULL;
    mpData = NULL;
    mgDataType = GDT_Byte;
    mnRows = mnCols = mnBands = -1;
    mnDatalength = -1;
    mpData = NULL;
    memset(mpGeoTransform, 0, 6*sizeof(double));
    strcpy(msProjectionRef, "");
    strcpy(msFilename, "");
    mdInvalidValue = 0.0f;
    mnPerPixSize = 1;

    //
    srcSR = NULL;
    latLongSR = NULL;
    poTransform = NULL;
    poTransformT = NULL;
}


CGDALRead::~CGDALRead(void)
{
    close();
}

void CGDALRead::close()
{
    if (mpoDataset != NULL)
    {
        GDALClose(mpoDataset);
        mpoDataset = NULL;      
    }

    if (mpData != NULL)
    {
        delete []mpData;
        mpData = NULL;
    }

    mgDataType = GDT_Byte;
    mnDatalength = -1;
    mnRows = mnCols = mnBands = -1;
    mpData = NULL;
    memset(mpGeoTransform, 0, 6*sizeof(double));
    strcpy(msProjectionRef, "");
    strcpy(msFilename, "");
    mdInvalidValue = 0.0f;
    mnPerPixSize = 1;

    //destory
    if (poTransform!=NULL)
        OCTDestroyCoordinateTransformation(poTransform);
    poTransform = NULL;

    if (poTransformT!=NULL)
        OCTDestroyCoordinateTransformation(poTransformT);
    poTransformT = NULL;

    if (latLongSR != NULL)
        OSRDestroySpatialReference(latLongSR);
    latLongSR = NULL;

    if (srcSR!=NULL)
        OSRDestroySpatialReference(srcSR);
    srcSR = NULL;


}

bool CGDALRead::isValid()
{
    if (mpoDataset == NULL || mpData == NULL)
    {
        return false;
    }

    return true;

}

GDALDataset* CGDALRead::poDataset()
{
    return mpoDataset;
}

size_t CGDALRead::rows()
{
    return mnRows;
}

size_t CGDALRead::cols()
{
    return mnCols;
}

size_t CGDALRead::bandnum()
{
    return mnBands;
}

unsigned char* CGDALRead::imgData()
{
    return mpData;
}

GDALDataType CGDALRead::datatype()
{
    return mgDataType;
}

double* CGDALRead::geotransform()
{
    return mpGeoTransform;
}

char* CGDALRead::projectionRef()
{
    return msProjectionRef;
}

size_t CGDALRead::datalength()
{
    return mnDatalength;
}

double CGDALRead::invalidValue()
{
    return mdInvalidValue;
}

size_t CGDALRead::perPixelSize()
{
    return mnPerPixSize;
}

bool CGDALRead::loadFrom( const char* _filename )
{
    //close fore image
    close();

    //register
    if(GDALGetDriverCount() == 0)
    {
        GDALAllRegister();
        OGRRegisterAll();
        CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
    }

    //open image
    mpoDataset = (GDALDataset*)GDALOpenShared(_filename, GA_ReadOnly);

    if (mpoDataset == NULL)
    {
        cout<<"CGDALRead::loadFrom : read file error!"<<endl;
        return false;
    }

    strcpy(msFilename, _filename);

    //get attribute
    mnRows = mpoDataset->GetRasterYSize();
    mnCols = mpoDataset->GetRasterXSize();
    mnBands = mpoDataset->GetRasterCount();
    mgDataType = mpoDataset->GetRasterBand(1)->GetRasterDataType();
    mdInvalidValue = mpoDataset->GetRasterBand(1)->GetNoDataValue();

    //mapinfo
    mpoDataset->GetGeoTransform(mpGeoTransform);
    strcpy(msProjectionRef, mpoDataset->GetProjectionRef());

    srcSR = OSRNewSpatialReference(msProjectionRef); // ground
    latLongSR = OSRCloneGeogCS(srcSR);  //geo
    poTransform =OCTNewCoordinateTransformation(srcSR, latLongSR);
    poTransformT =OCTNewCoordinateTransformation(latLongSR, srcSR);



    //get data
    bool bRlt = false;
    switch(mgDataType)
    {
    case GDT_Byte:
        mnPerPixSize = sizeof(unsigned char);
        bRlt = readData<unsigned char>();
        break;
    case GDT_UInt16:
        mnPerPixSize = sizeof(unsigned short);
        bRlt = readData<unsigned short>();
        break;
    case GDT_Int16:
        mnPerPixSize = sizeof(short);
        bRlt = readData<short>();
        break;
    case GDT_UInt32:
        mnPerPixSize = sizeof(unsigned int);
        bRlt = readData<unsigned int>();
        break;
    case GDT_Int32:
        mnPerPixSize = sizeof(int);
        bRlt = readData<int>();
        break;
    case GDT_Float32:
        mnPerPixSize = sizeof(float);
        bRlt = readData<float>();
        break;
    case GDT_Float64:
        mnPerPixSize = sizeof(double);
        bRlt = readData<double>();
        break;
    default:
        cout<<"CGDALRead::loadFrom : unknown data type!"<<endl;
        close();
        return false;
    }

    if (bRlt == false)
    {
        cout<<"CGDALRead::loadFrom : read data error!"<<endl;
        close();
        return false;
    }


    return true;
}

template<class TT> bool CGDALRead::readData()
{
    if (mpoDataset == NULL)
        return false;

    //new space
    mnDatalength = mnRows*mnCols*mnBands*sizeof(TT);
    mpData = new unsigned char[(size_t)mnDatalength];

    //raster IO
    CPLErr _err= mpoDataset->RasterIO(GF_Read, 0, 0, mnCols, mnRows, mpData, \
                mnCols, mnRows, mgDataType, mnBands, 0, 0, 0, 0);

    if (_err != CE_None)
    {
        cout<<"CGDALRead::readData : raster io error!"<<endl;
        return false;
    }

    return true;
}

unsigned char* CGDALRead::read( size_t _row, size_t _col, size_t _band )
{
    return &(mpData[(_band*mnRows*mnCols + _row*mnCols + _col)*mnPerPixSize]);
}

unsigned char* CGDALRead::readL( size_t _row, size_t _col, size_t _band )
{
    //if out of rect, take mirror
    if (_row < 0)
        _row = -_row;
    else if (_row >= mnRows)
        _row = mnRows - (_row - (mnRows - 1));

    if (_col < 0)
        _col = -_col;
    else if (_col >= mnCols)
        _col = mnCols - (_col - (mnCols - 1));

    return &(mpData[(_band*mnRows*mnCols + _row*mnCols + _col)*mnPerPixSize]);
}

template<class TT> double CGDALRead::linRead( double _row, double _col, size_t _band )
{
    TT val[4];
    double t1, t2, t, tx1, ty1;

    //calculate the excursion
    float xpos = _row - size_t(_row);
    float ypos = _col - size_t(_col);

    //get the pixel value of 4-neighbour
    tx1 = _row+1.0; ty1 = _col+1.0;

    val[0] =*(TT*)ReadL(size_t(_row), size_t(_col), bands); //band
    val[1] =*(TT*)ReadL(size_t(tx1), size_t(_col), bands);
    val[2] =*(TT*)ReadL(size_t(_row), size_t(ty1), bands);
    val[3] =*(TT*)ReadL(size_t(tx1), size_t(ty1), bands);

    //y-direction size_terpolation
    t1 = (1-ypos)*(double)val[0]+ypos*(double)val[2];
    t2 = (1-ypos)*(double)val[1]+ypos*(double)val[3];

    //x-direction size_terpolation
    t = (1-xpos)*t1 + xpos*t2;

    return (double)t;
}

bool CGDALRead::world2Pixel( double lat, double lon, double *x, double *y )
{
//  if (poTransformT==NULL)
//  {
//      poTransformT =OCTNewCoordinateTransformation(latLongSR, srcSR);
//  }

    if(poTransformT != NULL)
    {
        double height;
        OCTTransform(poTransformT,1, &lon, &lat, &height);

        double  adfInverseGeoTransform[6];
        GDALInvGeoTransform(mpGeoTransform, adfInverseGeoTransform);
        GDALApplyGeoTransform(adfInverseGeoTransform, lon,lat, x, y);

        return true;
    }
    else
    {
        return false;
    }

}

bool CGDALRead::pixel2World( double *lat, double *lon, double x, double y )
{
    if (poTransform!=NULL)
    {
        OCTDestroyCoordinateTransformation(poTransform);
        poTransform =OCTNewCoordinateTransformation(latLongSR, srcSR);
    }

    GDALApplyGeoTransform(mpGeoTransform, x, y, lon, lat);

    if(poTransform != NULL)
    {
        double height;
        OCTTransform(poTransform,1, lon, lat, &height);
        return true;
    }
    else
    {
        return false;
    }
}

bool CGDALRead::pixel2Ground( double x,double y,double* pX,double* pY )
{
    GDALApplyGeoTransform(mpGeoTransform, x, y, pX, pY);

    return true;
}

bool CGDALRead::ground2Pixel( double X,double Y,double* px,double* py )
{
    double  adfInverseGeoTransform[6];

    GDALInvGeoTransform(mpGeoTransform, adfInverseGeoTransform);
    GDALApplyGeoTransform(adfInverseGeoTransform, X, Y, px, py);

    return true;
}

GDALWrite.h

/************************************************************************/
/* This class is wrote for supercomputer to write image use parallel.
/* Not support Block read & write
/* But support multi-thread
/* Author: Y. Yao
/* E-mail: whuyao@foxmail.com
/* Version: v4.0
/************************************************************************/

#ifndef CLASS_GDAL_WRITE
#define CLASS_GDAL_WRITE

#include "gdal_priv.h"
#include "ogr_core.h"
#include "ogr_spatialref.h"

class CGDALRead;

class CGDALWrite
{
public:
    CGDALWrite(void);
    ~CGDALWrite(void);

public:
    bool init(const char* _filename, size_t _rows, size_t _cols, size_t _bandnum,\
                double _pGeoTransform[6], const char* _sProjectionRef, \
                GDALDataType _datatype = GDT_Byte, \
                double _dInvalidVal = 0.0f);

    bool init(const char* _filename, CGDALRead* pRead);

    bool init(const char* _filename, CGDALRead* pRead, size_t bandnum, \
                GDALDataType _datatype = GDT_Byte, \
                double _dInvalidVal = 0.0f);

    void write(size_t _row, size_t _col, size_t _band, void* pVal);

public:
    void close();

public:
    GDALDriver* poDriver();
    GDALDataset* poDataset();
    size_t rows();
    size_t cols();
    size_t bandnum();
    size_t datalength();
    double invalidValue();
    unsigned char* imgData();
    GDALDataType datatype();
    double* geotransform();
    char* projectionRef();
    size_t perPixelSize();

public:
    bool world2Pixel(double lat, double lon, double *pcol, double *prow);
    bool pixel2World(double *lat, double *lon, double col, double row);
    bool pixel2Ground(double col,double row,double* pX,double* pY);
    bool ground2Pixel(double X,double Y,double* pcol,double* prow);

protected:
    template<class TT> bool createData();

protected:
    GDALDriver* mpoDriver;      //can not release this, maybe cause some memory error!
    GDALDataset* mpoDataset;    //=>
    size_t mnRows;                  //
    size_t mnCols;                  //
    size_t mnBands;             //
    unsigned char* mpData;      //=>
    GDALDataType mgDataType;    //
    size_t mnDatalength;            //=>
    double mpGeoTransform[6];   //
    char msProjectionRef[2048]; //
    char msFilename[2048];      //
    double mdInvalidValue;      //
    size_t mnPerPixSize;            //=>

public:
    OGRSpatialReferenceH srcSR;
    OGRSpatialReferenceH latLongSR;
    OGRCoordinateTransformationH poTransform;   //pixel->world
    OGRCoordinateTransformationH poTransformT;  //world->pixel
};

#endif

GDALWrite.cpp

#include "GDALWrite.h"
#include "ogrsf_frmts.h"
#include "GDALRead.h"
#include <iostream>
using namespace std;


CGDALWrite::CGDALWrite(void)
{
    mpoDriver = NULL;
    mpoDataset = NULL;
    mnRows = mnCols = mnBands = -1;
    mpData = NULL;
    mgDataType = GDT_Byte;
    mnDatalength = 0;
    memset(mpGeoTransform, 0, 6*sizeof(double));
    strcpy(msProjectionRef, "");
    strcpy(msFilename, "");
    mdInvalidValue = 0.0f;
    mnPerPixSize = 1;

    //
    srcSR = NULL;
    latLongSR = NULL;
    poTransform = NULL;
    poTransformT = NULL;
}


CGDALWrite::~CGDALWrite(void)
{
    close();
}

void CGDALWrite::close()
{
    //write into data
    if (mpoDataset!=NULL && mpData!=NULL)
    {
        mpoDataset->RasterIO(GF_Write, 0, 0, mnCols, mnRows, \
            mpData, mnCols, mnRows, mgDataType, mnBands, 0, 0, 0, 0);
        mpoDataset->FlushCache();
    }


    ////release memory
    if (mpoDataset!=NULL)
    {
        GDALClose(mpoDataset);
        mpoDataset = NULL;
    }

    mnRows = mnCols = mnBands = -1;

    if (mpData!=NULL)
    {
        delete []mpData;
        mpData = NULL;
    }



    mgDataType = GDT_Byte;
    mnDatalength = 0;
    memset(mpGeoTransform, 0, 6*sizeof(double));
    strcpy(msProjectionRef, "");
    strcpy(msFilename, "");
    mdInvalidValue = 0.0f;
    mnPerPixSize = 1;

    //destory
    if (poTransform!=NULL)
        OCTDestroyCoordinateTransformation(poTransform);
    poTransform = NULL;

    if (poTransformT!=NULL)
        OCTDestroyCoordinateTransformation(poTransformT);
    poTransformT = NULL;

    if (latLongSR != NULL)
        OSRDestroySpatialReference(latLongSR);
    latLongSR = NULL;

    if (srcSR!=NULL)
        OSRDestroySpatialReference(srcSR);
    srcSR = NULL;

//  if (mpoDriver!=NULL)
//  {
//      //GDALDestroyDriver(mpoDriver);
//      delete mpoDriver;
//      mpoDriver = NULL;
//  }

}

GDALDriver* CGDALWrite::poDriver()
{
    return mpoDriver;
}

GDALDataset* CGDALWrite::poDataset()
{
    return mpoDataset;
}

size_t CGDALWrite::rows()
{
    return mnRows;
}

size_t CGDALWrite::cols()
{
    return mnCols;
}

size_t CGDALWrite::bandnum()
{
    return mnBands;
}

size_t CGDALWrite::datalength()
{
    return mnDatalength;
}

double CGDALWrite::invalidValue()
{
    return mdInvalidValue;
}

unsigned char* CGDALWrite::imgData()
{
    return mpData;
}

GDALDataType CGDALWrite::datatype()
{
    return mgDataType;
}

double* CGDALWrite::geotransform()
{
    return mpGeoTransform;
}

char* CGDALWrite::projectionRef()
{
    return msProjectionRef;
}

size_t CGDALWrite::perPixelSize()
{
    return mnPerPixSize;
}

bool CGDALWrite::init( const char* _filename, size_t _rows, size_t _cols, size_t _bandnum, double _pGeoTransform[6], const char* _sProjectionRef, GDALDataType _datatype /*= GDT_Byte*/, double _dInvalidVal /*= 0.0f*/ )
{
    close();

    //register
    if(GDALGetDriverCount() == 0)
    {
        GDALAllRegister();
        OGRRegisterAll();
        CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
    }

    //load
    mpoDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
    if (mpoDriver == NULL)
    {
        cout<<"CGDALWrite::init : Create poDriver Failed."<<endl;
        close();
        return false;
    }

    //
    strcpy(msFilename, _filename);
    mnRows = _rows;
    mnCols = _cols;
    mnBands = _bandnum;

    for (size_t i=0; i<6; i++)
        mpGeoTransform[i] = _pGeoTransform[i];

    strcpy(msProjectionRef, _sProjectionRef);
    mgDataType = _datatype;
    mdInvalidValue = _dInvalidVal;

    //create podataset
    char** papseMetadata = mpoDriver->GetMetadata();
    mpoDataset = mpoDriver->Create(msFilename, mnCols, mnRows, mnBands, mgDataType, papseMetadata);
    if (mpoDataset == NULL)
    {
        cout<<"CGDALWrite::init : Create poDataset Failed."<<endl;
        close();
        return false;       
    }

    //create others
    srcSR = OSRNewSpatialReference(msProjectionRef); // ground
    latLongSR = OSRCloneGeogCS(srcSR);  //geo
    poTransform =OCTNewCoordinateTransformation(srcSR, latLongSR);
    poTransformT =OCTNewCoordinateTransformation(latLongSR, srcSR);

    //add projection and coordinate
    poDataset()->SetGeoTransform(mpGeoTransform);
    poDataset()->SetProjection(msProjectionRef);
    for (size_t i =0; i<mnBands; i++)
    {
        poDataset()->GetRasterBand(i+1)->SetNoDataValue(mdInvalidValue);
    }

    //create data
    bool bRlt = false;
    switch(mgDataType)
    {
    case GDT_Byte:
        bRlt = createData<unsigned char>();
        break;
    case GDT_UInt16:
        bRlt = createData<unsigned short>();
        break;
    case GDT_Int16:
        bRlt = createData<short>();
        break;
    case GDT_UInt32:
        bRlt = createData<unsigned int>();
        break;
    case GDT_Int32:
        bRlt = createData<int>();
        break;
    case GDT_Float32:
        bRlt = createData<float>();
        break;
    case GDT_Float64:
        bRlt = createData<double>();
        break;
    default:
        cout<<"CGDALWrite::init : unknown data type!"<<endl;
        close();
        return false;
    }

    if (bRlt == false)
    {
        cout<<"CGDALWrite::init : Create data error!"<<endl;
        close();
        return false;
    }

    return true;
}

bool CGDALWrite::init( const char* _filename, CGDALRead* pRead )
{
    if (pRead == NULL)
    {
        cout<<"CGDALWrite::init : CGDALRead Point is Null."<<endl;
        return false;
    }

    return init(_filename, pRead->rows(), pRead->cols(), pRead->bandnum(), \
        pRead->geotransform(), pRead->projectionRef(), pRead->datatype(), 
        pRead->invalidValue());
}

bool CGDALWrite::init( const char* _filename, CGDALRead* pRead, size_t bandnum, GDALDataType _datatype /*= GDT_Byte*/, double _dInvalidVal /*= 0.0f*/ )
{
    if (pRead == NULL)
    {
        cout<<"CGDALWrite::init : CGDALRead Point is Null."<<endl;
        return false;
    }

    return init(_filename, pRead->rows(), pRead->cols(), bandnum, \
                pRead->geotransform(), pRead->projectionRef(), _datatype, 
                _dInvalidVal);

}

template<class TT> bool CGDALWrite::createData()
{
    if (mpoDataset == NULL)
        return false;

    if (mpData!=NULL)
        delete mpData;
    mpData = NULL;

    mnPerPixSize = sizeof(TT);
    mnDatalength = mnRows*mnCols*mnBands*mnPerPixSize;
    mpData = new unsigned char[mnDatalength];
    memset(mpData, 0, mnDatalength);
    return true;
}

void CGDALWrite::write( size_t _row, size_t _col, size_t _band, void* pVal )
{
    size_t nloc = (_band*mnRows*mnCols + _row*mnCols + _col)*mnPerPixSize;
    memcpy(mpData+nloc, pVal, mnPerPixSize);    
}

bool CGDALWrite::world2Pixel( double lat, double lon, double *x, double *y )
{
//  if (poTransformT==NULL)
//  {
//      poTransformT =OCTNewCoordinateTransformation(latLongSR, srcSR);
//  }


    if(poTransformT != NULL)
    {
        double height;
        OCTTransform(poTransformT,1, &lon, &lat, &height);

        double  adfInverseGeoTransform[6];
        GDALInvGeoTransform(mpGeoTransform, adfInverseGeoTransform);
        GDALApplyGeoTransform(adfInverseGeoTransform, lon,lat, x, y);

        return true;
    }
    else
    {
        return false;
    }
}

bool CGDALWrite::pixel2World( double *lat, double *lon, double x, double y )
{
    if (poTransform!=NULL)
    {
        OCTDestroyCoordinateTransformation(poTransform);
        poTransform =OCTNewCoordinateTransformation(latLongSR, srcSR);
    }

    GDALApplyGeoTransform(mpGeoTransform, x, y, lon, lat);

    if(poTransform != NULL)
    {
        double height;
        OCTTransform(poTransform,1, lon, lat, &height);
        return true;
    }
    else
    {
        return false;
    }
}

bool CGDALWrite::pixel2Ground( double x,double y,double* pX,double* pY )
{
    GDALApplyGeoTransform(mpGeoTransform, x, y, pX, pY);

    return true;
}

bool CGDALWrite::ground2Pixel( double X,double Y,double* px,double* py )
{
    double  adfInverseGeoTransform[6];

    GDALInvGeoTransform(mpGeoTransform, adfInverseGeoTransform);
    GDALApplyGeoTransform(adfInverseGeoTransform, X, Y, px, py);

    return true;
}