GDAL——生成等值线

时间:2022-09-14 03:29:46

参考李民录的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;
}