3

I need to crop a DEM raster file with a shapefile. This shapefile has many polygons. I need to calculate several variables after the crop operation.

The algorithm is written in C++ so I want to call GDAL C++ API. I do not need to write the warped file out to my disk and I intend to use either VRT or MEM once the test is passed.

I already obtain the coordinates of each polygon and I passed as a polygon using the cutline argument, but the operation has no effect so far.

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

I assume I need to construct the polygon vertex to image index according to: CUTLINE: This may contain the WKT geometry for a cutline. It will be converted into a geometry by GDALWarpOperation::Initialize() and assigned to the GDALWarpOptions hCutline field. The coordinates must be expressed in source pixel/line coordinates. Note: this is different from the assumptions made for the -cutline option of the gdalwarp utility! https://gdal.org/api/gdalwarp_cpp.html

7
  • What error messages are you getting? Which line is the error message being generated on? Commented Aug 20, 2020 at 2:44
  • The code does not give error. The the result tif is the same with the input, so the warp operation does not has effect.
    – Chang
    Commented Aug 20, 2020 at 3:01
  • Is your clip shape and raster input in the same coordinate reference system? If one is undefined or there isn't a geotransform you may see problems, though that scenario should trip a breakpoint. Commented Aug 20, 2020 at 3:26
  • My input files have the same map projection. This is the part I am not sure. The hCutline option is "Optional OGRPolygonH for a masking cutline", it is not a shapefile. But it can be defined from the shapefile. gdal.org/doxygen/…
    – Chang
    Commented Aug 20, 2020 at 4:24
  • OGR_G_AssignSpatialReference( hMultiPolygon, OGR_L_GetSpatialRef(hLayer) ); *phCutlineRet = hMultiPolygon; It appears the spatial reference is required for the geometry object.
    – Chang
    Commented Aug 20, 2020 at 5:11

1 Answer 1

1

After several attempts, I was able to fix the issues following some blog post: http://www.jeepxie.net/article/293437.html

The core part is that we need to convert coordinates:

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

and we should not use both hCutline and cutline at the same time.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.