8

I have read pixel values of an raster image using GDAL libraries in visual studio 2010(vc++).

Now I have to crop the image (subset) according to the grid given in shape file.

I have got a code for cropping image on internet as given below, but its having some error in it. Please review it and help .

GDALDataset* crop(GDALDataset* poDS, int xoff, int yoff, int xsize, int ysize)
{
        VRTDataset *poVDS = (VRTDataset*)VRTCreate(xsize, ysize);

        // Copy dataset info
        const char* pszProjection = poDS->GetProjectionRef();
        if (pszProjection != NULL && strlen(pszProjection) > 0)
                poVDS->SetProjection(pszProjection);

        double adfGeoTransform[6];
        if (poDS->GetGeoTransform(adfGeoTransform) == CE_None)
        {
                 // Adapt the geotransform matrix to the subarea
                    adfGeoTransform[0] += xoff * adfGeoTransform[1]
                        + yoff * adfGeoTransform[2];
                adfGeoTransform[3] += xoff * adfGeoTransform[4]
                        + yoff * adfGeoTransform[5];

                poVDS->SetGeoTransform(adfGeoTransform);
        }

        poVDS->SetMetadata(poDS->GetMetadata());

        // Here I also copy metadata from my own domain
        char **papszMD;
         papszMD = poDS->GetMetadata(MD_DOMAIN_MSAT);
               if (papszMD != NULL)
                 poVDS->SetMetadata(papszMD, MD_DOMAIN_MSAT);

         for (int i = 0; i < poDS->GetRasterCount(); ++i)
         {
                 GDALRasterBand* poSrcBand = poDS->GetRasterBand(i + 1);
                 GDALDataType eBandType = poSrcBand->GetRasterDataType();
                 poVDS->AddBand(eBandType, NULL);
                 VRTSourcedRasterBand* poVRTBand = (VRTSourcedRasterBand*)poVDS- >GetRasterBand(i + 1);
                 poVRTBand->AddSimpleSource(poSrcBand,
                                 xoff, yoff,
                                 xsize, ysize,
                                 0, 0, xsize, ysize);
                 poVRTBand->CopyCommonInfoFrom(poSrcBand);

                 // Again, I copy my own metadata
                 papszMD = poSrcBand->GetMetadata(MD_DOMAIN_MSAT);
                 if (papszMD != NULL)
                         poVRTBand->SetMetadata(papszMD, MD_DOMAIN_MSAT);
         } 

         return poVDS;
 }
Improve this question
1

3 Answers 3

Reset to default
12

Unprojected:

gdal_translate -srcwin xstart ystart xstop ystop input.raster output.raster

Projected:

gdal_translate -projwin ulx uly lrx lry input.raster output.raster

If you are using gdal on your system, gdal_translate is installed.

Improve this answer
5
  • Hello jay laura: thank you for the reply but the lines you have suggested works for python i guess. i only uses c++ code. So i dont think gdal_translate work in c++ coding .
    – Hardik
    Commented Jan 8, 2013 at 17:01
  • 1
    @hardik, gdal_translate is written in C++ and operates independently of the gdal python bindings, here's the code: trac.osgeo.org/gdal/browser/trunk/gdal/apps/gdal_translate.cpp. It's part of the core utility set and a compiled binary is generally installed by default (unless you're not installing and only building from source).
    – matt wilkie
    Commented Mar 22, 2013 at 16:22
  • 1
    The unprojected command seems inconsistent with the documentation gdal.org/gdal_translate.html, which says -srcwin xoff yoff xsize ysize.
    – user55937
    Commented Feb 17, 2016 at 18:44
  • Can you explain what you mean by "projected" and "unprojected". My raster has a projection, but I don't want to reproject it — which form is correct?
    – colllin
    Commented Dec 2, 2019 at 18:59
  • @colllin "projected" means you can specify the coordinates using the coordinate format of the projection (e.g. for example if using EPSG:4326 you can provide 90/180 lat/lon coordinates) and "unprojected" means you specify coordinates as x/y pixel values of the raster. For projected, you can specify the SRS using -projwin_srs. It does not actually reproject anything in the output.
    – rooby
    Commented May 21, 2021 at 8:08
4

See the GDALRasterBand Class. Carefully using the nXOff and nXSize (etc.) parameters of GDALRasterBand::RasterIO will read a "cropped" part of the raster.

Improve this answer
3

Thank you Mike. I have solved the problem of cropping the image using VC++ with gdal libraries. I have created VRTDataset of my desired size of raster to be cropped and then save it using CreateCopy().

Improve this answer

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.