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;
 }
1
11

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.

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
    Jan 8 '13 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). Mar 22 '13 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
    Feb 17 '16 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
    Dec 2 '19 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
    May 21 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.

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().

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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