Обрезка растрового файла с помощью GDAL warp C ++ API

Aug 20 2020

Мне нужно обрезать растровый файл DEM с помощью шейп-файла. В этом шейп-файле много многоугольников. Мне нужно рассчитать несколько переменных после операции с урожаем.

Алгоритм написан на C ++, поэтому я хочу вызвать GDAL C ++ API. Мне не нужно записывать искаженный файл на мой диск, и я намерен использовать либо VRT, либо MEM после прохождения теста.

Я уже получил координаты каждого многоугольника и передал его как многоугольник, используя аргумент линии пореза, но операция пока не имеет никакого эффекта.

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

Я предполагаю, что мне нужно построить вершину многоугольника для индекса изображения в соответствии с: CUTLINE: Это может содержать геометрию WKT для линии пореза. Он будет преобразован в геометрию с помощью GDALWarpOperation :: Initialize () и назначен полю GDALWarpOptions hCutline. Координаты должны быть выражены в исходных координатах пикселя / линии. Примечание: это отличается от предположений, сделанных для опции -cutline утилиты gdalwarp!https://gdal.org/api/gdalwarp_cpp.html

Ответы

1 Chang Aug 24 2020 at 04:20

После нескольких попыток я смог исправить проблемы после публикации в блоге: http://www.jeepxie.net/article/293437.html

Основная часть заключается в том, что нам нужно преобразовать координаты:

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

и мы не должны использовать одновременно hCutline и cutline.