这个封装的目的主要是为了简化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;
}