[gdal-dev] wrong with clipping a tif file with a polygon

[Ahmet Temiz]

Hi,
I am trying to clip a tif file with a polygon. But resulting output.tif
file shows nothing.

Where am I doing wrong ?

"  #include <gdal/gdal.h>
#include <gdal/gdal_priv.h>
#include <gdal/gdalwarper.h>
#include <gdal/ogrsf_frmts.h>
#include <gdal/gdalwarper.h>
#include <string> // for string class
using namespace std;

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 = "input.geojson";

  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
  hDstDS->GetRasterBand(1)->Fill(0);

  hDstDS->SetProjection(hSrcDS->GetProjectionRef());

  const char *src_srs = hSrcDS->GetProjectionRef();
  const char *dst_srs = hDstDS->GetProjectionRef();
  printf("%s ",src_srs);
  // 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();
  //char* wkt = nullptr;
  if (geom != nullptr && geom->IsValid())
  {
      printf("\n geom is valid \n");
  }
  // 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);
}
"

-- 
Ahmet Temiz
Jeoloji Müh.
Afet ve Acil Durum Yönetimi Başkanlığı
Deprem  Dairesi Başkanlığı


________________________

Ahmet Temiz
Geological Eng.

Disaster and Emergency Management
of Presidency
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20230530/6d88e7e7/attachment-0001.htm>