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

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

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

[cpp] view plaincopyprint?
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;  

}
调用的代码，如下：

[cpp] view plaincopyprint?
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的，没有什么大的问题。