OSGeo
gdalwarp VRT NoData issue
What is the bug?
Gdalwarp in c++ when using VRT seems to ignore nans despite telling gdal explicitly
Steps to reproduce the issue
I have a simple 61x61 dataset generated randomly :
NoData value is nan, when I project it using this command line :
gdalwarp -s_srs "+proj=utm +zone=15 +north +datum=WGS84 +units=m +no_defs" -t_srs "+proj=utm +zone=15 +north +datum=WGS84 +units=m +no_defs" -tr 60 60 -r bilinear -srcnodata nan -ot Float64 -overwrite -of GTiff input_generated.tif output_gdal.tif
everything works fine as you can see :
bilinear interpolation didn't take into account nan values due to "-srcnodata nan"
now when I try to do it programmaticaly with a VRT :
- opening the dataset
// Create the in-memory dataset
GDALDriver* memDriver = GetGDALDriverManager()->GetDriverByName("MEM");
dataset.reset(memDriver->Create("", nCols, nRows, 1, GDT_Float64, NULL));
// Get the raster band
GDALRasterBand* band = dataset->GetRasterBand(1);
// Set NoData value for NaN handling
band->SetNoDataValue(std::numeric_limits<double>::quiet_NaN());
// Copy data from xarray to GDAL dataset
// Note: GDAL expects data in row-major order, same as xtensor default
CPLErr result = band->RasterIO(GF_Write, 0, 0, nCols, nRows, const_cast<double*>(data.data()),
nCols, nRows, GDT_Float64, 0, 0);
return dataset;
Then projection using the required NoData options :
projectedDataset.reset(
memDriver->Create("", nCols, nRows, inputDataset->GetRasterCount(), GDT_Float64, NULL));
if (projectedDataset == NULL) {
LogHandler::handleLog(
"I-008", vector<string>{"Failed to open projected dataset for " + datasetType, "open"});
return;
}
// Set the projection to UTM
OGRSpatialReference projectedSRS;
projectedSRS.SetWellKnownGeogCS("WGS84");
projectedSRS.SetUTM(ProjectionInfo::get().getUtmZone(), ProjectionInfo::get().getNorthing());
projectedDataset->SetSpatialRef(&projectedSRS);
auto geoTransform = ProjectionInfo::get().getGeoTransform();
// Set the geo transform
projectedDataset->SetGeoTransform(geoTransform.data());
std::unique_ptr<GDALWarpOptions, decltype(&GDALDestroyWarpOptions)> warpOptions{
nullptr, GDALDestroyWarpOptions};
// get Nodata value
GDALRasterBand *srcBand = inputDataset->GetRasterBand(1);
int hasNoData;
double noDataValue = srcBand->GetNoDataValue(&hasNoData);
cout << "NoData value: " << noDataValue << endl;
cout << "hasNoData prep value: " << hasNoData << endl;
int nBands = inputDataset->GetRasterCount();
cout << "nBands: " << nBands << endl;
// Setup warp options
warpOptions.reset(GDALCreateWarpOptions());
warpOptions->hSrcDS = inputDataset;
warpOptions->hDstDS = projectedDataset.get();
// CRITICAL: Set these options to handle nodata properly during bilinear interpolation
warpOptions->papszWarpOptions =
CSLSetNameValue(warpOptions->papszWarpOptions, "INIT_DEST", "NO_DATA");
warpOptions->papszWarpOptions =
CSLSetNameValue(warpOptions->papszWarpOptions, "UNIFIED_SRC_NODATA", "YES");
warpOptions->papszWarpOptions =
CSLSetNameValue(warpOptions->papszWarpOptions, "SAMPLE_GRID", "YES");
// Only set nodata if the source actually has nodata defined
if (hasNoData == 1) {
// Allocate arrays for source and destination nodata values
warpOptions->padfSrcNoDataReal = (double *)CPLMalloc(sizeof(double) * nBands);
warpOptions->padfDstNoDataReal = (double *)CPLMalloc(sizeof(double) * nBands);
cout << "warp options allocated" << endl;
// Set nodata values for each band
for (int i = 0; i < nBands; i++) {
// Also set nodata on the destination bands
GDALRasterBand *dstBand = projectedDataset->GetRasterBand(i + 1);
dstBand->SetNoDataValue(noDataValue);
cout << "setting nodata value for output band " << i << " to " << noDataValue << endl;
cout << "setting nodata value for band " << i << endl;
warpOptions->padfSrcNoDataReal[i] = noDataValue;
warpOptions->padfDstNoDataReal[i] = noDataValue;
}
} else {
// If no nodata is defined, you might want to set a specific value
// or leave these as NULL to avoid nodata handling issues
cout << "No nodata value defined in source dataset" << endl;
}
warpOptions->pfnProgress = GDALTermProgress;
// Establish reprojection transformer
warpOptions->pTransformerArg = GDALCreateGenImgProjTransformer(
inputDataset, GDALGetProjectionRef(inputDataset), projectedDataset.get(),
GDALGetProjectionRef(projectedDataset.get()), FALSE, 0.0, 1);
warpOptions->pfnTransformer = GDALGenImgProjTransform;
warpOptions->eResampleAlg = resampleAlg;
warpOptions->pfnProgress = GDALDummyProgress; // Add this to suppress progress messages
warpOptions->pProgressArg = nullptr;
// Initialize and execute the warp operation
try {
LogCapture capture("GDAL");
capture.start();
GDALWarpOperation operation;
operation.Initialize(warpOptions.get());
operation.ChunkAndWarpImage(0, 0, GDALGetRasterXSize(projectedDataset.get()),
GDALGetRasterYSize(projectedDataset.get()));
capture.stop();
} catch (const std::exception &e) {
throw GenericException("I-016", {"GDAL ", string(e.what())});
}
The prints inform me that noData is indeed nan and that it was correctly set inside input and output dataset, now the output :
we can see that nan values propagated, so they were not correctly interpreted as nans, despite NoData being shown as nan when opening them inside qgis.
What am I doing wrong? Or is it a bug? ( Gdal downloaded from conan center, version 3.8.3 )
Versions and provenance
3.8.3
Additional context
No response
it's not actually a problem of "MEM" driver, I have the same issue when opening the dataset from a .tif file and writing to a .tif dataset.
Another interesting information : if I switch nan values to 9999 and set NoData value as 9999 : gdalwarp tool still works, but c++ warp interprets 9999 as values and doesn't work.
Apart from warpOptions->nBandCount = 1 which is missing, nothing obvious strikes me as wrong. I'd suggest you use the GDALWarp() higher level API. Otherwise please include a fully self-contained C/C++ reproducer (ie something that just needs to be compiled)
Thanks for your answer, somehow warpOptions->nBandCount = 1 made the program idle forever...
I asked LLMS for a simpler version using GDALReprojectImage :
void GenerateUtils::testProjectionWithGDALWarp() {
int nRows = 61, nCols = 61;
double noData = 9999.0;
xt::xarray<double> data = xt::random::randn<double>({nRows, nCols}) * 10.0 + 283.0;
data = xt::where(xt::random::rand<double>({nRows, nCols}) > 0.5, data, noData);
string inputPath = "test_input.tif";
RasterUtils::writeSingleBandXarrayToGeoTiff<double>(data, inputPath, true, noData);
GDALDataset* srcDS = (GDALDataset*)GDALOpen(inputPath.c_str(), GA_ReadOnly);
OGRSpatialReference srcSRS, dstSRS;
srcSRS.SetWellKnownGeogCS("WGS84");
srcSRS.SetUTM(15, 1);
srcDS->SetSpatialRef(&srcSRS);
double srcGeoTransform[6] = {699960, 1800.0, 0.0, 5700000, 0.0, -1800.0};
srcDS->SetGeoTransform(srcGeoTransform);
dstSRS.SetWellKnownGeogCS("WGS84");
dstSRS.SetUTM(15, 1);
auto geoTransform = ProjectionInfo::get().getGeoTransform();
// Define output dimensions BEFORE creating VRT
int outRows = 1830;
int outCols = 1830;
char* dstWKT = nullptr;
dstSRS.exportToWkt(&dstWKT);
GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* finalDS =
driver->Create("test_output.tif", outCols, outRows, 1, GDT_Float64, nullptr);
finalDS->SetSpatialRef(&dstSRS);
finalDS->SetGeoTransform(geoTransform.data());
finalDS->GetRasterBand(1)->SetNoDataValue(noData);
GDALReprojectImage(srcDS, nullptr, finalDS, dstWKT, GRA_Bilinear, 0.0, 0.0, nullptr, nullptr,
nullptr);
xt::xarray<double> result = xt::zeros<double>({outRows, outCols});
finalDS->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, outCols, outRows, result.data(), outCols,
outRows, GDT_Float64, 0, 0);
cout << "Result shape: " << xt::adapt(result.shape()) << endl;
GDALClose(srcDS);
GDALClose(finalDS);
CPLFree(dstWKT);
exit(0);
}
and it worked, posting here so other people can see the solution.
I thought it was OK, but the issue came from somewhere else, it came from the fact that the input raster had multiple bands :
I made this reproductible small code for testing
#include "gdal.h"
#include "gdal_priv.h"
#include "gdalwarper.h"
#include "xtensor/xadapt.hpp"
#include "xtensor/xarray.hpp"
#include "xtensor/xio.hpp"
#include "xtensor/xmanipulation.hpp"
#include "xtensor/xmath.hpp"
#include "xtensor/xrandom.hpp"
#include "xtensor/xview.hpp"
void GenerateUtils::testProjectionWithGDALWarpGeogToUtm(int nBands) {
int nRows = 61, nCols = 61;
double noData = std::numeric_limits<double>::quiet_NaN();
xt::xarray<double> data[nBands];
// Create test data and double it
for (int i = 0; i < nBands; i++) {
data[i] = xt::random::randn<double>({nRows, nCols}) * 10.0 + 283.0;
data[i] = xt::where(xt::random::rand<double>({nRows, nCols}) > 0.5, data[i], noData);
}
// Basic geographic projection setup
const char* srcWKT =
"GEOGCS[\"WGS84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS "
"84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]]";
double srcGeoTransform[6] = {-91.0, 0.05, 0.0, 51.5, 0.0, -0.025};
// Create input dataset
GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* srcDS = driver->Create(("test_input_" + to_string(nBands) + "_bands.tif").c_str(),
nCols, nRows, nBands, GDT_Float64, nullptr);
OGRSpatialReference srcSRS;
char* pSrcWKT = CPLStrdup(srcWKT);
srcSRS.importFromWkt(&pSrcWKT);
srcDS->SetSpatialRef(&srcSRS);
srcDS->SetGeoTransform(srcGeoTransform);
for (int i = 1; i <= srcDS->GetRasterCount(); i++) {
srcDS->GetRasterBand(i)->SetNoDataValue(noData);
srcDS->GetRasterBand(i)->RasterIO(GF_Write, 0, 0, nCols, nRows, data[i - 1].data(), nCols,
nRows, GDT_Float64, 0, 0);
}
// Setup UTM destination
OGRSpatialReference dstSRS;
dstSRS.SetWellKnownGeogCS("WGS84");
dstSRS.SetUTM(15, 1);
// Create output dataset
int outRows = 915, outCols = 915;
auto geoTransform = ProjectionInfo::get().getGeoTransform();
GDALDataset* dstDS = driver->Create(("test_output_" + to_string(nBands) + "_bands.tif").c_str(),
outCols, outRows, nBands, GDT_Float64, nullptr);
dstDS->SetSpatialRef(&dstSRS);
dstDS->SetGeoTransform(geoTransform.data());
// Initialize destination bands with NoData
for (int i = 1; i <= dstDS->GetRasterCount(); i++) {
dstDS->GetRasterBand(i)->SetNoDataValue(noData);
// Initialize the entire band with NoData values
double* buffer = new double[outCols * outRows];
std::fill_n(buffer, outCols * outRows, noData);
dstDS->GetRasterBand(i)->RasterIO(GF_Write, 0, 0, outCols, outRows, buffer, outCols,
outRows, GDT_Float64, 0, 0);
delete[] buffer;
}
// Setup warp options
char* dstWKT = nullptr;
dstSRS.exportToWkt(&dstWKT);
char** papszWarpOptions = NULL;
papszWarpOptions = CSLSetNameValue(papszWarpOptions, "INIT_DEST", "NO_DATA");
papszWarpOptions = CSLSetNameValue(papszWarpOptions, "WRITE_FLUSH", "YES");
papszWarpOptions = CSLSetNameValue(papszWarpOptions, "SKIP_NOSOURCE", "YES");
papszWarpOptions = CSLSetNameValue(papszWarpOptions, "UNIFIED_SRC_NODATA", "NO");
// Create and configure GDALWarpOptions
GDALWarpOptions* warpOptions = GDALCreateWarpOptions();
warpOptions->papszWarpOptions = papszWarpOptions;
warpOptions->hSrcDS = srcDS;
warpOptions->hDstDS = dstDS;
warpOptions->eResampleAlg = GRA_Bilinear;
// Set band count and allocate band arrays
int bandCount = srcDS->GetRasterCount();
warpOptions->nBandCount = bandCount;
// Allocate band mapping arrays
warpOptions->panSrcBands = (int*)CPLMalloc(sizeof(int) * bandCount);
warpOptions->panDstBands = (int*)CPLMalloc(sizeof(int) * bandCount);
// Set band mapping (1-based indexing for GDAL)
for (int i = 0; i < bandCount; i++) {
warpOptions->panSrcBands[i] = i + 1;
warpOptions->panDstBands[i] = i + 1;
}
// Allocate and set NoData values for each band
warpOptions->padfSrcNoDataReal = (double*)CPLMalloc(sizeof(double) * bandCount);
warpOptions->padfDstNoDataReal = (double*)CPLMalloc(sizeof(double) * bandCount);
warpOptions->padfSrcNoDataImag = nullptr;
warpOptions->padfDstNoDataImag = nullptr;
// Set NoData values for each band
for (int i = 0; i < bandCount; i++) {
warpOptions->padfSrcNoDataReal[i] = noData;
warpOptions->padfDstNoDataReal[i] = noData;
}
// Perform reprojection
CPLErr err = GDALReprojectImage(srcDS, nullptr, dstDS, dstWKT, GRA_Bilinear, 0.0, 0.0, nullptr,
nullptr, warpOptions);
// Force write to disk
GDALFlushCache(dstDS);
// Cleanup
GDALDestroyWarpOptions(warpOptions);
GDALClose(srcDS);
GDALClose(dstDS);
}
this function generate a random array of values ( nRows * nCols), for each band if you call it with 1 band it will work perfectly :
but as soon as you have 2 or more bands :
it seems like papszWarpOptions = CSLSetNameValue(papszWarpOptions, "UNIFIED_SRC_NODATA", "NO"); is completely ignored
it seems like
papszWarpOptions = CSLSetNameValue(papszWarpOptions, "UNIFIED_SRC_NODATA", "NO");is completely ignored
From the image, it isn't obvious at all to me what the bug is. That may be an issue with your viewer or how you set it up ? If you want to make a reproducer to analyze, you should create very small rasters, like input_band1 = [1, 2, 1, 3] input_band2 = [1, 1, 2, 4] and nodata_band1 = nodata_band2 = 1