Ritaglio di file raster con l'API GDAL warp C ++

Question

Devo ritagliare un file raster DEM con uno shapefile. Questo shapefile ha molti poligoni. Devo calcolare diverse variabili dopo l'operazione di ritaglio.

L'algoritmo è scritto in C ++, quindi voglio chiamare l'API GDAL C ++. Non ho bisogno di scrivere il file deformato sul mio disco e intendo utilizzare VRT o MEM una volta superato il test.

Ho già ottenuto le coordinate di ogni poligono e ho passato come poligono utilizzando l'argomento cutline, ma l'operazione non ha alcun effetto fino ad ora.

            OGRLinearRing poExteriorRing;
            for (i = 0; i < nPt; i++) 
              {
                dX = (*iIterator).vVertex.at(i).dX;
                dY = (*iIterator).vVertex.at(i).dY;
                double dDummy1 = (dX - dX_origin) / dResolution_elevation;
                long lColumn_index = long(round(dDummy1));
                double dDummy2 = (dY_origin - dY) / dResolution_elevation;
                long lRow_index = long(round(dDummy2));     
                poExteriorRing.addPoint(lColumn_index, lRow_index);
              }
            poExteriorRing.closeRings();
            OGRPolygon polygon;
            polygon.addRing(&poExteriorRing);

            GDALDataset  *poVRTDS;
            poDS_elevation->GetGeoTransform(adfGeoTransform);
            poVRTDS = poDriver->CreateCopy( "./test.tiff", poDS_elevation, FALSE, NULL, NULL, NULL );

            GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
            ...
            psWarpOptions->hCutline = &polygon;                
            GDALWarpOperation oOperation;
            oOperation.Initialize( psWarpOptions );
            oOperation.ChunkAndWarpImage( 0, 0,
                            GDALGetRasterXSize( poDS_elevation ),
                            GDALGetRasterYSize( poDS_elevation ) );
            GDALClose( poVRTDS );

Presumo di dover costruire il vertice del poligono sull'indice dell'immagine in base a: CUTLINE: Questo può contenere la geometria WKT per una linea di taglio. Verrà convertito in una geometria da GDALWarpOperation :: Initialize () e assegnato al campo GDALWarpOptions hCutline. Le coordinate devono essere espresse in coordinate pixel / linea di origine. Nota: questo è diverso dalle ipotesi fatte per l'opzione -cutline dell'utilità gdalwarp!https://gdal.org/api/gdalwarp_cpp.html

Answer 1

Dopo diversi tentativi, sono stato in grado di risolvere i problemi seguendo alcuni post sul blog: http://www.jeepxie.net/article/293437.html

La parte fondamentale è che dobbiamo convertire le coordinate:

eError = TransformCutlineToSource( hWrkSrcDS, hCutline,
                                  &(psWO->papszWarpOptions),
                                  psOptions->papszTO );

e non dovremmo usare contemporaneamente hCutline e cutline.