使用GDAL下载并转换SRTM的DEM数据（二）

之前写了一篇使用GDAL下载SRTM的数据，只是大概说明了下怎么下载，文章后面说要写代码实现一下，最近太慢，一直拖着，这事肯定不会忘记的，呵呵。今天就把上次剩下的尾巴处理一下。之前的博客地址：http://blog.csdn.net/liminlu0314/article/details/8068715。

在上篇博客中大概分析了一下，首先要通过代码构造一个VRT文件。函数代码如下：

int CreateSRTMVRTFile(const char* pszSrcFile, const char* pszDstFile, const char* pszWkt,
					  int iWidth, int iHeight, int iDataType, double *dGeoTransform,
					  int iImageOffset, int iLineOffset, int iPixelOffset, 
					  const char* pszFormat, CProcessBase* pProcess)
{
	if (pProcess != NULL)
		pProcess->SetMessage("正在执行SRTM数据转换");

	GDALAllRegister();

	//首先将裸数据转换为一个VRT文件格式，然后用VRT导出为指定的类型
	GDALDriver *poDriver = (GDALDriver *) GDALGetDriverByName( "VRT" );
	if (poDriver == NULL)
	{
		if (pProcess != NULL)
			pProcess->SetMessage("不支持VRT格式创建！");
		return RE_FILENOTSUPPORT;
	}

	string strVrt = GetTemporaryPath();
	string strImageOffset = lexical_cast<string>(iImageOffset);
	string strLineOffset = lexical_cast<string>(iLineOffset);
	string strPixelOffset = lexical_cast<string>(iPixelOffset);

	GDALDataset *poVRTDS;
	poVRTDS = poDriver->Create( strVrt.c_str(), iWidth, iHeight, 0, (GDALDataType)iDataType, NULL );
	if (poVRTDS == NULL)
	{
		if (pProcess != NULL)
			pProcess->SetMessage("不支持VRT格式创建！");
		return RE_FILENOTSUPPORT;
	}

	poVRTDS->SetGeoTransform(dGeoTransform);
	poVRTDS->SetProjection(pszWkt);

	char** papszOptions = NULL;
	// 指定为RawRasterBand，if not specified, default to VRTRasterBand
	papszOptions = CSLAddNameValue(papszOptions, "subclass", "VRTRawRasterBand");
	// 指定原始数据
	papszOptions = CSLAddNameValue(papszOptions, "SourceFilename", pszSrcFile); 
	// optionnal. default = false
	papszOptions = CSLAddNameValue(papszOptions, "RelativeToVRT", "true");
	// optionnal. default = machine order
	papszOptions = CSLAddNameValue(papszOptions, "ByteOrder", "MSB"); 
	// optionnal. default = 0 
	papszOptions = CSLAddNameValue(papszOptions, "ImageOffset", strImageOffset.c_str()); 
	// optionnal. default = size of band type 
	papszOptions = CSLAddNameValue(papszOptions, "PixelOffset", strPixelOffset.c_str());
	// optionnal. default = size of band type * width 
	papszOptions = CSLAddNameValue(papszOptions, "LineOffset", strLineOffset.c_str()); 
	poVRTDS->AddBand((GDALDataType)iDataType, papszOptions);
	CSLDestroy(papszOptions);

	GDALDriver *poDstDriver = (GDALDriver *) GDALGetDriverByName(pszFormat);
	if (poDstDriver == NULL)
	{
		if (pProcess != NULL)
			pProcess->SetMessage("不支持输出格式创建！");
		return RE_FILENOTSUPPORT;
	}

	GDALDataset *poDstDS = poDstDriver->CreateCopy(pszDstFile, poVRTDS, FALSE, NULL, GDALTermProgress, pProcess);
	if (poDstDS == NULL)
	{
		if (pProcess != NULL)
			pProcess->SetMessage("输出文件创建失败！");
		return RE_CREATEFAILED;
	}

	poDstDS->FlushCache();
	delete poVRTDS;
	RasterDelete(strVrt.c_str());	//移除临时vrt文件

	GDALClose(poDstDS);

	if (pProcess != NULL)
		pProcess->SetMessage("SRTM数据转换完成！");

	return RE_SUCCESS;
}

调用的代码，如下：

	const char* pszInFile = "/vsicurl/ftp://topex.ucsd.edu/pub/srtm30_plus/srtm30/data/e020n40.Bathymetry.srtm";
	const char* pszDstFile = "D:\\e020n40_SRTM.tif";
	const char* pszWkt = SRS_WKT_WGS84;
	const char* pszFormat = "GTiff";
	double dGeoTransform[6] = {20, 0.008333333333333,  0.0,  40,  0.0, -0.008333333333333};

	iRev = CreateSRTMVRTFile(pszInFile, pszDstFile, pszWkt, 4800, 6000, 2, dGeoTransform, 0, 9600, 2, pszFormat, pPro);

上面的函数中用到了Boost库，当然必不可少的GDAL库，这两个库包进去应该可以编译过去。

刚测试的时候发现这个SRTM的FTP（ftp://topex.ucsd.edu/pub/srtm30_plus/srtm30/data/）不能访问了，导致数据不能下载下来。你可以把路径换成别的路径进行测试，应该没有问题的。最主要的是要把这个二进制的数据的元数据信息指定正确，比如长宽，波段数，数据类型，行偏移量等，其他的信息就是投影坐标等。这些信息设置正确的话，基本上就是OK的，没有什么大的问题。