1

我正在尝试使用 GDAL 的多边形来剪辑栅格。目前我收到一个错误，即在初始化 WarpOperation 时存在读取访问冲突。我可以访问我的 Shapefile 并检查功能的数量，所以我认为访问很好。我也可以访问我的栅格数据（GetProjectionRef）。所有文件都在同一个 CRS 中。有没有办法将 GdalWarp 与 Cutline 一起使用？

      const char* inputPath = "input.tif";
      const char* outputPath = "output.tif";
     
      //clipper Polygon
      auto w_read_filenamePoly = "Polygon.shp";
      char* read_filenamePoly = new char[w_read_filenamePoly.length() + 1];
      wcstombs(read_filenamePoly, w_read_filenamePoly.c_str(), w_read_filenamePoly.length() + 1);

      GDALDataset* hSrcDS; 
      GDALDataset* hDstDS;

      GDALAllRegister();
      hSrcDS =(GDALDataset *) GDALOpen(inputPath, GA_Update);
      hDstDS = (GDALDataset*)GDALOpen(outputPath, GA_Update);
     
      const char* proj = hSrcDS->GetProjectionRef();
      const char* proj2 = hDstDS->GetProjectionRef();
      //clipper Layer
      GDALDataset* poDSClipper;
      poDSClipper = (GDALDataset*)GDALOpenEx(read_filenamePoly, GDAL_OF_UPDATE, NULL, NULL, NULL);
      Assert::IsNotNull(poDSClipper);
      delete[]read_filenamePoly;
      OGRLayer* poLayerClipper;
      poLayerClipper = poDSClipper->GetLayerByName("Polygon");
      int numClip = poLayerClipper->GetFeatureCount();

      //setup warp options
      GDALWarpOptions* psWarpOptions = GDALCreateWarpOptions();
      psWarpOptions->hSrcDS = hSrcDS;
      psWarpOptions->hDstDS = hDstDS;
      psWarpOptions->nBandCount = 1;
      psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
      psWarpOptions->panSrcBands[0] = 1;
      psWarpOptions->panDstBands = (int*)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
      psWarpOptions->panDstBands[0] = 1;
      psWarpOptions->pfnProgress = GDALTermProgress;
      psWarpOptions->hCutline = poLayerClipper;
      
      
      // Establish reprojection transformer.
      psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer(hSrcDS,proj, hDstDS, proj2, FALSE, 0.0, 1);
      psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
     
      GDALWarpOperation oOperation;
      oOperation.Initialize(psWarpOptions);
      oOperation.ChunkAndWarpImage(0, 0, GDALGetRasterXSize(hDstDS), GDALGetRasterYSize(hDstDS));
      GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
      GDALDestroyWarpOptions(psWarpOptions);
      GDALClose(hDstDS);
      GDALClose(hSrcDS);
4

1 回答 1

0

你psWarpOptions->hCutline应该是一个多边形，而不是一个图层。

切割线也应该在源像素/线坐标中。

检查TransformCutlineToSourcefrom gdalwarp_lib.cpp，您可能可以简单地从那里获取代码。

这个特殊的 GDAL 操作，当从 C++ 调用时，充满了陷阱 - 这里有很多关于它的未解决问题 - 我正在复制一个完整的工作示例：

使用多边形蒙版（切割线）扭曲（重新投影）光栅图像：

#include <gdal/gdal.h>
#include <gdal/gdal_priv.h>
#include <gdal/gdalwarper.h>
#include <gdal/ogrsf_frmts.h>

int main() {
  const char *inputPath = "input.tif";
  const char *outputPath = "output.tif";

  // clipper Polygon
  // THIS FILE MUST BE IN PIXEL/LINE COORDINATES or otherwise one should
  // copy the function gdalwarp_lib.cpp:TransformCutlineToSource()
  // from GDAL's sources
  // It is expected that it contains a single polygon feature
  const char *read_filenamePoly = "cutline.json";

  GDALDataset *hSrcDS;
  GDALDataset *hDstDS;

  GDALAllRegister();
  auto poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
  hSrcDS = (GDALDataset *)GDALOpen(inputPath, GA_ReadOnly);

  hDstDS = (GDALDataset *)poDriver->CreateCopy(
    outputPath, hSrcDS, 0, nullptr, nullptr, nullptr);
  // Without this step the cutline is useless - because the background
  // will be carried over from the original image
  CPLErr e = hDstDS->GetRasterBand(1)->Fill(0);

  const char *src_srs = hSrcDS->GetProjectionRef();
  const char *dst_srs = hDstDS->GetProjectionRef();

  // clipper Layer
  GDALDataset *poDSClipper;
  poDSClipper = (GDALDataset *)GDALOpenEx(
    read_filenamePoly, GDAL_OF_UPDATE, NULL, NULL, NULL);
  auto poLayerClipper = poDSClipper->GetLayer(0);
  auto geom = poLayerClipper->GetNextFeature()->GetGeometryRef();

  // setup warp options
  GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
  psWarpOptions->hSrcDS = hSrcDS;
  psWarpOptions->hDstDS = hDstDS;
  psWarpOptions->nBandCount = 1;
  psWarpOptions->panSrcBands =
    (int *)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
  psWarpOptions->panSrcBands[0] = 1;
  psWarpOptions->panDstBands =
    (int *)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
  psWarpOptions->panDstBands[0] = 1;
  psWarpOptions->pfnProgress = GDALTermProgress;
  psWarpOptions->hCutline = geom;

  // Establish reprojection transformer.
  psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer(
    hSrcDS, src_srs, hDstDS, dst_srs, TRUE, 1000, 1);
  psWarpOptions->pfnTransformer = GDALGenImgProjTransform;

  GDALWarpOperation oOperation;
  oOperation.Initialize(psWarpOptions);
  oOperation.ChunkAndWarpImage(
    0, 0, GDALGetRasterXSize(hDstDS), GDALGetRasterYSize(hDstDS));
  GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
  GDALDestroyWarpOptions(psWarpOptions);
  GDALClose(hDstDS);
  GDALClose(hSrcDS);
}
于 2022-01-11T11:45:41.523 回答