GDAL生成等高线——等值线
本文主要介绍：利用gdal的函数，根据DEM图像，生成等高线或等值线，介绍两种方法：一种是利用GDAL自带的exe文件，一种是利用GDAL函数。
说明：GDAL使用版本为Gdal2.0.0。
1 利用GDAL自带exe，生成等高线
1.1 参数说明
具体说明可见：GDAL实用工具简介
官网说明：gdal_contour.exe
1.2 调用exe程序
//数字 转 string
string Float2String(float num, string format)
{
char temp[100];
sprintf_s(temp, format.c_str(), num);//保留2位小数
string str = temp;
return str;
}
/********************************************************************
函数功能：调用gdal自带exe生成等高线
输入参数：
imgPath ：dem文件路径
shpPath ：等高线矢量图保存路径
dfContourInterval ：等高线间隔
备注：
需要头文件： #include <windows.h>
*********************************************************************/
bool CalContourByExe(string imgPath, string shpPath, double dfContourInterval)
{
string temp = "gdal_contour.exe -a Height "; //海拔高度对应的字段名为Height 也可以改为其它名字
temp = imgPath + " " + shpPath + " " + "-i " + Float2String(dfContourInterval, "%.1f");;
STARTUPINFO si;
memset(&si,0,sizeof(STARTUPINFO));//初始化si在内存块中的值（详见memset函数）
si.cb=sizeof(STARTUPINFO);
si.dwFlags=STARTF_USESHOWWINDOW;
si.wShowWindow=SW_SHOW;
PROCESS_INFORMATION pi;//必备参数设置结束
WCHAR wcharTemp[256]; //string 转 LPWSTRT
MultiByteToWideChar(CP_ACP,0,temp.c_str(),-1,wcharTemp,sizeof(wcharTemp)/sizeof(wcharTemp[0]));
if(!CreateProcess(NULL, wcharTemp,NULL,NULL,FALSE,CREATE_NO_WINDOW,NULL,NULL,&si,&pi))
{
cout<<"运行等高线exe失败"<<endl;
exit(1);
}
else
{
}
//不使用的句柄最好关掉
CloseHandle(pi.hThread);
CloseHandle(pi.hProcess);
return true;
}
1.3 exe程序存在问题
（1）DEM文件路径和生成的shp文件路径中，不能含有中文；
（2）不能处理含有异常值的DEM文件；
2 利用GDAL函数生成等高线
2.1 主要用到GDAL函数
CPLErr GDALContourGenerate (GDALRasterBandH hBand,
double dfContourInterval,
double dfContourBase,
int nFixedLevelCount,
double * padfFixedLevels,
int bUseNoData,
double dfNoDataValue,
void * hLayer,
int iIDField,
int iElevField,
GDALProgressFunc pfnProgress,
void * pProgressArg
) 该函数官网说明：GDALContourGenerate
理解：
CPLErr GDALContourGenerate (GDALRasterBandH hBand, //图像波段
double dfContourInterval, //等高线间隔
double dfContourBase, //等高线起始高度
int nFixedLevelCount, //固定高度值，如果为非0值，ContourInterval and ContourBase将不起作用
double * padfFixedLevels, //一组固定高度值，
int bUseNoData, //是否使用无效值 true使用
double dfNoDataValue, //如果使用，则生成等值线时，会忽略无效值，该数值为无效值对应的数值
//因此如果dem中有无效值，需要先确定无效值大小
void * hLayer, //矢量图层
int iIDField, //如果非-1，ID号会保存在指定属性字段
int iElevField, //如果非-1，高程值会保存在指定属性字段
GDALProgressFunc pfnProgress,
void * pProgressArg
) 2.2 生成等高线完整程序
/********************************************************************
函数功能：调用gdal自带exe生成等高线
输入参数：
imgPath ：dem文件路径
shpPath ：等高线矢量图保存路径
dfContourInterval ：等高线间隔
备注：
需要头文件： #include "gdal_priv.h"
#include "gdal_alg.h"
#include "gdal_priv.h"
#include "ogrsf_frmts.h"
*********************************************************************/
bool CalContourByGdal(string imgPath, string shpPath, double dfContourInterval)
{
OGRRegisterAll();
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");
GDALDataset * pInDataset = (GDALDataset * )GDALOpen(imgPath.c_str(),GA_ReadOnly);
if(pInDataset==NULL)
{
printf("读取图像失败！");
return FALSE;
}
int nDemWidth = pInDataset->GetRasterXSize(); //获取图像宽
int nDemHeight = pInDataset->GetRasterYSize(); //获取图像高
int Count = pInDataset->GetRasterCount(); //波段数
//读取图像数据波段
GDALRasterBand *pInRasterBand = pInDataset->GetRasterBand(1);
float *pData = new float[nDemWidth*nDemHeight]();
CPLErr err = pInRasterBand->RasterIO(GF_Read, 0, 0, nDemWidth, nDemHeight, pData, nDemWidth, nDemHeight, GDT_Float32, 0, 0);
if(err == CE_Failure)
{
cout << "读取DEM图像数据时出错！" << endl;
GDALDestroyDriverManager();
return 0;
}
//判断图像中是否有异常值，并获取异常值实际值
float fNoData = 0;
int nIdx;
for(int i=0; i<nDemHeight; i++)
{
for(int j=0; j<nDemWidth; j++)
{
nIdx = i*nDemWidth + j;
if(pData[nIdx] <= -9999)
{
fNoData = pData[nIdx];
}
}
}
//创建矢量图
const char *pszDriverName = "ESRI Shapefile";
GDALDriver *poDriver;
GDALAllRegister();
poDriver = GetGDALDriverManager()->GetDriverByName(pszDriverName );
if( poDriver == NULL )
{
printf( "%s driver not available.\n", pszDriverName );
exit( 1 );
}
GDALDataset *poDS;
poDS = poDriver->Create( shpPath.c_str(), 0, 0, 0, GDT_Unknown, NULL ); //创建数据源
if( poDS == NULL )
{
printf( "Creation of output file failed.\n" );
exit( 1 );
}
OGRLayer *poLayer;
OGRSpatialReference *poSpatialRef = new OGRSpatialReference(pInDataset->GetProjectionRef());
poLayer = poDS->CreateLayer( "Elevation", poSpatialRef, wkbLineString, NULL ); //创建图层
if( poLayer == NULL )
{
printf( "Layer creation failed.\n" );
exit( 1 );
}
OGRFieldDefn ofieldDef1("Elevation", OFTInteger); //在矢量图中创建高程值字段
if(poLayer->CreateField(&ofieldDef1) != OGRERR_NONE)
{
cout<<"创建矢量图层属性表失败！"<<endl;
GDALClose((GDALDatasetH)poDS);
GDALClose(pInDataset);
return 0;
}
//根据图像波段生成矢量图等高线
if(fNoData == 0)
GDALContourGenerate(pInRasterBand, dfContourInterval, 0, 0, NULL, false, 0, (OGRLayerH)poLayer, -1, 0, NULL, NULL);
else //有异常值时，不对异常值进行处理
GDALContourGenerate(pInRasterBand, dfContourInterval, 0, 0, NULL, true, fNoData, (OGRLayerH)poLayer, -1, 0, NULL, NULL);
GDALClose(poDS);
GDALClose(pInDataset);
if(pData != NULL)
{
delete[] pData;
pData = NULL;
}
}3 调用主函数
int _tmain(int argc, _TCHAR* argv[])
{
string imgPath ="F:\\2.tif";
string shpPath = "C:\\contour.shp";
CalContourByGdal(imgPath, shpPath, 500);
cout<<"处理完成"<<endl;
getchar();
return 0;
}注：完整程序下载 http://download.csdn.net/detail/hong__fang/9543727
转载自：https://blog.csdn.net/hong__fang/article/details/51605030
