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

Rahkonen Jukka jukka.rahkonen at maanmittauslaitos.fi
Tue May 30 06:20:15 PDT 2023


Hi,

I can’t really read the code, but based on what happens, maybe your cutline is outside the image. So you

  *   create an output band
  *   fill it with zeroes
  *   and then you write pixels within the cutline into the output
Because the output stays all zeroes, and there are no errors, then I guess that the pixels within your cutline are all nodata.
Maybe printing the WKT out from “geom” could prove or disprove my guess.

-Jukka Rahkonen-

Lähettäjä: gdal-dev <gdal-dev-bounces at lists.osgeo.org> Puolesta Ahmet Temiz
Lähetetty: tiistai 30. toukokuuta 2023 15.57
Vastaanottaja: Norman Vine <nhv at meganet.net>
Kopio: gdal-dev <gdal-dev at lists.osgeo.org>
Aihe: Re: [gdal-dev] wrong with clipping a tif file with a polygon

Thank you ,

I already worked on that code. But clipping was not succesful. output.tif is blank on qgis.

On Tue, May 30, 2023 at 3:24 PM Norman Vine <nhv at meganet.net<mailto:nhv at meganet.net>> wrote:
perhaps this will help
https://stackoverflow.com/questions/70664043/clip-raster-with-polygon-with-gdal-c

> On May 30, 2023, at 8:12 AM, Ahmet Temiz <ahmettemiz88 at gmail.com<mailto:ahmettemiz88 at gmail.com>> wrote:
>
> 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
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>
> https://lists.osgeo.org/mailman/listinfo/gdal-dev



--
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/c5671cf6/attachment-0001.htm>


More information about the gdal-dev mailing list