0

I am trying to clip a raster using a polygon an GDAL. The following code is running without an error but the ouputimage is all zeros. All files are in the same CRS. Is there a way to use GdalWarp with 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();

  //get my geometry as OGRPolygon
  std::vector<OGRFeature*> vecClipFeature;
  OGRFeature* feat;
  while ((feat = poLayerClipper->GetNextFeature()) != NULL)
  {
    vecClipFeature.push_back(feat);
  }
  OGRGeometry* geom = vecClipFeature.at(0)->GetGeometryRef();
  OGRPolygon *clipPoly = geom->toPolygon();

  //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 = clipPoly;
  
  
  // 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);
1

1 Answer 1

0

This is a duplicate of a question on SO. There are many pitfalls when using the C++ interface.

#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);
}

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.