[gdal-dev] wrong with clipping a tif file with a polygon
Ahmet Temiz
ahmettemiz88 at gmail.com
Tue May 30 05:12:15 PDT 2023
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>
More information about the gdal-dev
mailing list