Recorte de archivos ráster con GDAL warp C ++ API

Aug 20 2020

Necesito recortar un archivo ráster DEM con un shapefile. Este shapefile tiene muchos polígonos. Necesito calcular varias variables después de la operación de cultivo.

El algoritmo está escrito en C ++, así que quiero llamar a la API de GDAL C ++. No necesito escribir el archivo deformado en mi disco y tengo la intención de usar VRT o MEM una vez que se pase la prueba.

Ya obtengo las coordenadas de cada polígono y las pasé como polígono usando el argumento de línea de corte, pero la operación no ha tenido ningún efecto hasta ahora.

            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 );

Supongo que necesito construir el vértice del polígono al índice de imagen de acuerdo con: CUTLINE: Esto puede contener la geometría WKT para una línea de corte. Se convertirá en una geometría por GDALWarpOperation :: Initialize () y se asignará al campo GDALWarpOptions hCutline. Las coordenadas deben expresarse en coordenadas de línea / píxel de origen. Nota: ¡esto es diferente de las suposiciones hechas para la opción -cutline de la utilidad gdalwarp!https://gdal.org/api/gdalwarp_cpp.html

Respuestas

1 Chang Aug 24 2020 at 04:20

Después de varios intentos, pude solucionar los problemas después de una publicación de blog: http://www.jeepxie.net/article/293437.html

La parte central es que necesitamos convertir coordenadas:

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

y no debemos usar hCutline y cutline al mismo tiempo.