5

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;
 }
8

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.

  • 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). – matt wilkie Mar 22 '13 at 16:22
  • 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
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().

2

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

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.