使用 gdal - 编程之家

如何解决使用 gdal

在一些 C++ 代码中，我使用 gdal 读取了 geotiff 文件。我想要做的是设置从 GeoTIFFs 坐标系到 wgs84 (lat,long) 的坐标转换。

这在某些情况下有效，但对于其他一些 GeoTIFF，轴会被交换。我希望 gdal 库能够处理所有的转换工作，而无需我处理这个问题，但是对于我的代码，这并不能按预期工作。

这是相关代码行的摘录。我希望这些勾勒出相关的计算，即使这不是 MWE。我做了一些错误检查，但为了简洁起见删除了它。

    // data structures used below
    double GT_[6];

    OGRSpatialReference srcref_;
    OGRSpatialReference wgs84ref_;

    OGRCoordinateTransformation *transform_ = nullptr;
    OGRCoordinateTransformation *backtransform_ = nullptr;

    // some initializations,and open file
    GDALAllRegister();
    wgs84ref_.SetWellKNownGeogCS("epsg:4326");

    data_set_ = static_cast<GDALDataset *>(
        GDALOpen(filename_.c_str(),GA_ReadOnly));

    data_set_->GetGeoTransform(GT_)

    // get projection and set up transformatuon matrices to/from wgs84
    const char *pt =  data_set_->GetProjectionRef();
    auto res = srcref_.importFromWkt(const_cast<char **>(&pt));
    transform_ = OGRCreateCoordinateTransformation(&srcref_,&wgs84ref_);
    backtransform_ = OGRCreateCoordinateTransformation(&wgs84ref_,&srcref_);

后来我逐行读取数据，确定“本地”坐标如下：

    inline double geox(int x,int y) const {
        return GT_[0] + x*GT_[1] + y*GT_[2];
    }
    inline double geoy(int x,int y) const {
        return GT_[3] + x*GT_[4] + y*GT_[5];
    }


    GDALRasterBand *band = data_set_->GetRasterBand(b);
    int nx = band->GetXSize();
    int ny = band->GetYSize();

    for (int nyoffset = 0; nyoffset < ny; nyoffset++) {
  
        auto res = band->Rasterio(GF_Read,nyoffset,nx,1,scanline,GEOTIFF_TYPE,0);

        for (int x=0; x < nx; x++) {
            lat[x] = gx[x] = geox(x,nyoffset);
            lon[x] = gy[x] = geoy(x,nyoffset);
        }
        transform_->Transform(nx,&lat[0],&lon[0])
         
        // do some processing involving (lat,lon)
    }

出现的问题是某些文件的纬度和经度被交换了。这似乎与坐标传递给 Transform 的顺序有关，或者我如何从 GT_ 计算这些。当我在“问题 TIFF”上运行 gdalinfo 时，我看到：

Coordinate System is:
GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],...
Data axis to CRS axis mapping: 2,3

我如何正确看待这个“轴映射”？我希望 gdal 有办法为我处理这一切，但我不知道如何修复上面的代码。从 GT_ 计算 gx 和 gx 的想法是错误的，还是我需要明确考虑轴映射？

（我链接到 libgdal.so.28，gdalinfo 说 GDAL 3.2.1,released 2020/12/29）

更新：我现在添加了代码来“修复”轴映射：

    // figure out axis mapping
    auto spatialref = data_set_->GetSpatialRef();
    auto axismapping = spatialref->GetDataAxisToSRSAxisMapping();

    if (axismapping.size() < 2)
        throw std::string("Cannot find axis mapping information");

    if (axismapping[0] == 1 && axismapping[1] == 2) 
        axes_swapped_ = false;
    else if (axismapping[0] == 2 && axismapping[1] == 1) 
        axes_swapped_ = true;  
    else 
        throw std::string("I don't understand this axis mapping");

及以后：

     // fill gx and gy vectors and transform these into wgs84
     if ( !axes_swapped_ ) {
         for (int x=0; x < nx; x++) {
             lat[x] = gx[x] = geox(x,nyoffset);
             lon[x] = gy[x] = geoy(x,nyoffset);
         }
     } else {
         for (int x=0; x < nx; x++) {
             lon[x] = gx[x] = geox(x,nyoffset);
             lat[x] = gy[x] = geoy(x,nyoffset);
         }
     }
     if (!transform_->Transform(nx,&lon[0]))
         throw std::string("transform Failed");

这似乎有效，但我想知道我是否考虑了所有情况，以及 gdal 是否提供了一些功能来为我完成这一切。