参考李民录的gdal源代码剖析那本书写的,由于代码主要是用C的库,下面我使用C++相应的库进行重写,有的函数用法需要稍作修改,直接上代码:
//生成等高线
///C++
int CreateContourDlg::Createontour(const char* pszSrcDEM, const char* pszDstShp, int iBandIndex, double dInterval,
const char * pszFormat, CProgressBase* pProcess)
{
if(pProcess != NULL)
{
pProcess->ReSetProcess();
pProcess->SetProgressTip("CreateContour...");
}
GDALAllRegister();
OGRRegisterAll();
GDALDataset* poDataset;
poDataset = (GDALDataset*) GDALOpen(pszSrcDEM, GA_ReadOnly);
GDALRasterBand* poBand;
poBand = poDataset->GetRasterBand(iBandIndex);
// 获取原始DEM数据中的NODATA值
int bNoDataSet = FALSE, bIgnoreNoData = FALSE;
double dfNoData = 0.0;
if( !bNoDataSet && !bIgnoreNoData )
dfNoData = poBand->GetNoDataValue(&bNoDataSet);
OGRSpatialReference oSRS;
const char* pszWKT;
pszWKT = poDataset->GetProjectionRef();
char* psz;
psz = const_cast<char*>(pszWKT);
if( pszWKT != NULL && strlen(pszWKT) != 0 )
oSRS.importFromWkt(&psz);
//OGRSpatialReference oSRS(pszWKT);
GDALDriver* poDriverDem;
poDriverDem = GetGDALDriverManager()->GetDriverByName(pszFormat);
GDALDataset* poDS;
poDS = poDriverDem->Create(pszDstShp, 0, 0, 0, GDT_Unknown, 0);
OGRLayer* poLayer;
poLayer = poDS->CreateLayer("Contour", &oSRS, wkbLineString ,NULL);
OGRFieldDefn oFieldID( "ID", OFTInteger);
oFieldID.SetWidth(8);
poLayer->CreateField(&oFieldID, false);
OGRFieldDefn oFieldEle( "Elevation", OFTReal);
oFieldEle.SetWidth(12);
oFieldEle.SetPrecision(3);
poLayer->CreateField(&oFieldEle, false);
int nElevField = 1;
// 调用GDAL库中的函数生成等高线
MyGDALProgressFunc pfnProgress = ALGTermProgress;
CPLErr eErr = GDALContourGenerate( (GDALRasterBandH)poBand, dInterval, 0.0,
0, NULL,
bNoDataSet, dfNoData,
poLayer, 0, nElevField,
pfnProgress, pProcess );
oSRS.Release();
GDALClose(poDataset);
return 1;
}