[GIS] How to subset raster image using gdal

cgdalraster

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

Best Answer

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.

Related Question