使用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的，没有什么大的问题。