Warping a Raster File to Match Another in C++ with GDAL

In C++, how to warp a raster file to the same spatial reference, resolution and extent with another raster file and output an vrt file

To warp a raster file to match the spatial reference, resolution, and extent of another raster file in C++ and output a VRT (Virtual Raster) file, you can use the GDAL (Geospatial Data Abstraction Library) library. GDAL provides a powerful set of tools for raster and vector data manipulation.

Here's a step-by-step guide to achieve this:

1. Install GDAL: Make sure GDAL is installed on your system. You can download it from the GDAL website or install it using a package manager like apt on Linux or brew on macOS.

2. Include GDAL in Your Project: Ensure that your C++ project is set up to include the GDAL headers and link against the GDAL library.

3. Write the C++ Code: Below is an example of how you can warp a raster file to match another raster file and output a VRT file.

#include <gdal_priv.h>
#include <cpl_conv.h> // for CPLMalloc()

int main(int argc, char **argv) {
    if (argc < 4) {
        printf("Usage: %s <source_raster> <reference_raster> <output_vrt>\n", argv[0]);
        return 1;
    }

    const char *srcFilename = argv[1];
    const char *refFilename = argv[2];
    const char *dstFilename = argv[3];

    GDALAllRegister();

    // Open the source raster
    GDALDataset *srcDataset = (GDALDataset *) GDALOpen(srcFilename, GA_ReadOnly);
    if (srcDataset == nullptr) {
        printf("Failed to open source raster: %s\n", srcFilename);
        return 1;
    }

    // Open the reference raster
    GDALDataset *refDataset = (GDALDataset *) GDALOpen(refFilename, GA_ReadOnly);
    if (refDataset == nullptr) {
        printf("Failed to open reference raster: %s\n", refFilename);
        GDALClose(srcDataset);
        return 1;
    }

    // Get the spatial reference, resolution, and extent from the reference raster
    double refGeoTransform[6];
    refDataset->GetGeoTransform(refGeoTransform);
    const char *refProjection = refDataset->GetProjectionRef();
    int refXSize = refDataset->GetRasterXSize();
    int refYSize = refDataset->GetRasterYSize();

    // Set up the warp options
    GDALWarpOptions *warpOptions = GDALCreateWarpOptions();
    warpOptions->hSrcDS = srcDataset;
    warpOptions->hDstDS = nullptr; // We will create the VRT dataset later
    warpOptions->nBandCount = srcDataset->GetRasterCount();
    warpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) * warpOptions->nBandCount);
    warpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) * warpOptions->nBandCount);
    for (int i = 0; i < warpOptions->nBandCount; i++) {
        warpOptions->panSrcBands[i] = i + 1;
        warpOptions->panDstBands[i] = i + 1;
    }

    // Create the VRT dataset
    GDALDriver *vrtDriver = GetGDALDriverManager()->GetDriverByName("VRT");
    GDALDataset *vrtDataset = vrtDriver->Create(dstFilename, refXSize, refYSize, 0, GDT_Unknown, nullptr);
    vrtDataset->SetGeoTransform(refGeoTransform);
    vrtDataset->SetProjection(refProjection);

    // Set the destination dataset in the warp options
    warpOptions->hDstDS = vrtDataset;

    // Perform the warp
    GDALWarpOperation warpOperation;
    warpOperation.Initialize(warpOptions);
    CPLErr err = warpOperation.ChunkAndWarpImage(0, 0, refXSize, refYSize);
    if (err != CE_None) {
        printf("Warping failed.\n");
    }

    // Clean up
    GDALDestroyWarpOptions(warpOptions);
    GDALClose(srcDataset);
    GDALClose(refDataset);
    GDALClose(vrtDataset);

    return (err == CE_None) ? 0 : 1;
}

Explanation:

1. GDAL Initialization: GDALAllRegister() registers all GDAL drivers.
2. Open Source and Reference Rasters: GDALOpen() opens the source and reference raster files.
3. Get Reference Raster Properties: GetGeoTransform(), GetProjectionRef(), GetRasterXSize(), and GetRasterYSize() retrieve the spatial reference, resolution, and extent of the reference raster.
4. Set Up Warp Options: GDALCreateWarpOptions() creates warp options, and we set the source and destination bands.
5. Create VRT Dataset: GetGDALDriverManager()->GetDriverByName("VRT") gets the VRT driver, and Create() creates the VRT dataset with the same properties as the reference raster.
6. Perform the Warp: GDALWarpOperation performs the warp operation.
7. Clean Up: GDALDestroyWarpOptions(), GDALClose() clean up resources.

Compilation:

Make sure to link against the GDAL library when compiling your code. For example:

g++ -o warp_raster warp_raster.cpp -lgdal

This example assumes you have basic knowledge of C++ and GDAL. Adjust paths and error handling as needed for your specific use case.