Cropping images with GDAL

eng pdo sw tips

I am working to get a better integration between Meteosatlib and GDAL.

A nice aspect of GDAL is that it allows to create read/write drivers around two functions: Open and CreateCopy.

Open opens a datset read only, and to implement that all you have to do is to implement read access to your data using the GDALDataset interface.

To implement CreateCopy, all you have to do is to create a new file reading information from a GDALDataset; then call Open on it. This means that there is no need to support incremental updates, and that all the data required to create a new file is readily available. This simplifies matters a lot.

GDAL provides some interesting image manipulation functions, that can work over just these two calls. The way it does it is by exploiting the concept of a virtual dataset, which wraps an existing dataset but changes some parameters on the fly.

This means that you can wrap a read only dataset with a virtual dataset that transforms it somehow, and then pass the virtual dataset to CreateCopy to save the transformed image in a format of your choice.

On top of that, for many transformations you do not need to create your own virtual datasets, but you can use the functions provided by the VRT GDAL driver.

One can learn a lot on how to use VRTDataset by reading the source code for gdal_translate. By doing that I could come up with this code for cropping an image, which is an interesting VRTDataset example.

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

This function wraps a dataset with a virtual dataset that crops it. Just pass the resulting dataset to GDALCreateCopy to save it in the format that you need.