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 <[email protected]> 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 > [email protected] > https://lists.osgeo.org/mailman/listinfo/gdal-dev _______________________________________________ gdal-dev mailing list [email protected] https://lists.osgeo.org/mailman/listinfo/gdal-dev
