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

/************************************************************************/
/* 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

#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 CGDAL
12a38
Read::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;
}

/************************************************************************/
/* 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

#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;
}