NCDF: use of GeoTransform attribute if lon/lat dimensions are available

[rhijmans]

What is the bug?

When reading this NCDF file GDAL uses the GeoTransform parameters that are specified in the grid_mapping variable, instead of deriving them from the lon/lat dimensions.

In this case these are not the same because the file was edited with CDO after editing it with GDAL. From the driver documentation. It is not clear to me if this is by design (it seems to me that it is not).

This was reported as an issue with R/terra, see the discussion here.

The problem goes away when removing the GeoTransform attribute. It should be possible to avoid the problem when writing the NCDF file with GDAL by using "WRITE_GDAL_TAGS=NO" but that is easily missed.

If the current behavior is as intended, perhaps an option to use the GeoTransform only as a last resort could be added?

Steps to reproduce the issue

# terra::describe is a wrapper to gdalinfo
terra::describe("SWE_snodas_20250401_sub.nc")[63:66] 
#[1] "Upper Left  (-124.7333333,  52.8750000) (124d44' 0.00\"W, 52d52'30.00\"N)"
#[2] "Lower Left  (-124.7333333,  44.8750000) (124d44' 0.00\"W, 44d52'30.00\"N)"
#[3] "Upper Right (-116.0000000,  52.8750000) (116d 0' 0.00\"W, 52d52'30.00\"N)"
#[4] "Lower Right (-116.0000000,  44.8750000) (116d 0' 0.00\"W, 44d52'30.00\"N)"

# remove the GeoTransform attribute
nc <- ncdf4::nc_open("SWE_snodas_20250401_sub.nc", write=TRUE)
ncdf4::ncatt_put( nc, "crs", "GeoTransform",  "")
ncdf4::nc_close(nc)

terra::describe("SWE_snodas_20250401_sub.nc")[63:66]
#[1] "Upper Left  (-124.7333333,  49.5000000) (124d44' 0.00\"W, 49d30' 0.00\"N)"
#[2] "Lower Left  (-124.7333333,  41.5000000) (124d44' 0.00\"W, 41d30' 0.00\"N)"
#[3] "Upper Right (-116.0000000,  49.5000000) (116d 0' 0.00\"W, 49d30' 0.00\"N)"
#[4] "Lower Right (-116.0000000,  41.5000000) (116d 0' 0.00\"W, 41d30' 0.00\"N)"

Versions and provenance

Windows 10;
GDAL 3.11.0 via RTools

Additional context

No response

[dbaston]

Hi @rhijmans , this behavior is intended as of GDAL 3.11. For some background, see #11993

Still, GDAL should emit a warning where the stored GeoTransform is inconsistent with the one calculated from dimension variables. Are you seeing anything?

I can see an argument for preferring the stored GeoTransform only when it is reasonably consistent with the one computed from dimension variables.

[rhijmans]

Thanks @dbaston. The warning was suppressed, I now see:

gdal(warn=1)
x <- rast(f)
Warning message:
GeoTransform read from attribute of crs variable differs from value calculated from dimension variables (max diff = 3.375). Using value from attribute. (GDAL 1) 

The GeoTransform is not part of the CF convention, and other tools do not modify or remove it when changing the dimensions. For that reason, I think it would be sensible to use the dimension variables if there is a discrepancy like this.