From c4b96d03dc524a98c93c2c1255a3ce6707787a48 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Mon, 14 Oct 2024 15:50:03 +0200
Subject: [PATCH] netCDF CreateCopy(): fix taking into account NETCDF_DIM_EXTRA
 when source dataset is not georeferenced

Fixes https://lists.osgeo.org/pipermail/gdal-dev/2024-October/059618.html
---
 autotest/gdrivers/netcdf.py    | 24 +++++++++++++++++++++++-
 frmts/netcdf/netcdfdataset.cpp | 11 ++++++-----
 2 files changed, 29 insertions(+), 6 deletions(-)

diff --git a/autotest/gdrivers/netcdf.py b/autotest/gdrivers/netcdf.py
index 17863588251b..d79d303dfa02 100755
--- a/autotest/gdrivers/netcdf.py
+++ b/autotest/gdrivers/netcdf.py
@@ -6542,7 +6542,7 @@ def test_netcdf_create_metadata_with_equal_sign(tmp_path):
 # Test force opening a HDF55 file with netCDF driver
 
 
-def test_netcdf_force_opening_hdf5_file(tmp_vsimem):
+def test_netcdf_force_opening_hdf5_file():
 
     ds = gdal.OpenEx("data/hdf5/groups.h5", allowed_drivers=["netCDF"])
     assert ds.GetDriver().GetDescription() == "netCDF"
@@ -6559,3 +6559,25 @@ def test_netcdf_force_opening_no_match():
 
     drv = gdal.IdentifyDriverEx("data/byte.tif", allowed_drivers=["netCDF"])
     assert drv is None
+
+
+###############################################################################
+
+
+def test_netcdf_extra_dim_no_georef(tmp_path):
+
+    fname = tmp_path / "out.nc"
+
+    src_ds = gdal.Translate("", "../gcore/data/stefan_full_rgba.tif", format="MEM")
+    size_z = 4
+    src_ds.SetMetadataItem("NETCDF_DIM_EXTRA", "{Z}")
+    src_ds.SetMetadataItem("NETCDF_DIM_Z_DEF", f"{{{size_z},4}}")
+    src_ds.SetMetadataItem("NETCDF_DIM_Z_VALUES", "{0,1,2,3}")
+    src_ds.SetMetadataItem("Z#axis", "Z")
+
+    # Create netCDF file
+    gdal.GetDriverByName("netCDF").CreateCopy(fname, src_ds)
+
+    ds = gdal.Open(fname)
+    assert ds.RasterCount == 4
+    assert ds.ReadRaster() == src_ds.ReadRaster()
diff --git a/frmts/netcdf/netcdfdataset.cpp b/frmts/netcdf/netcdfdataset.cpp
index bfdd50e372ee..25c253cca45d 100644
--- a/frmts/netcdf/netcdfdataset.cpp
+++ b/frmts/netcdf/netcdfdataset.cpp
@@ -9604,11 +9604,6 @@ netCDFDataset::CreateCopy(const char *pszFilename, GDALDataset *poSrcDS,
         pScaledProgress =
             GDALCreateScaledProgress(0.1, 0.25, pfnProgress, pProgressData);
         poDS->AddProjectionVars(false, GDALScaledProgress, pScaledProgress);
-        // Save X,Y dim positions.
-        panDimIds[nDim - 1] = poDS->nXDimID;
-        panBandDimPos[0] = nDim - 1;
-        panDimIds[nDim - 2] = poDS->nYDimID;
-        panBandDimPos[1] = nDim - 2;
         GDALDestroyScaledProgress(pScaledProgress);
     }
     else
@@ -9622,6 +9617,12 @@ netCDFDataset::CreateCopy(const char *pszFilename, GDALDataset *poSrcDS,
         }
     }
 
+    // Save X,Y dim positions.
+    panDimIds[nDim - 1] = poDS->nXDimID;
+    panBandDimPos[0] = nDim - 1;
+    panDimIds[nDim - 2] = poDS->nYDimID;
+    panBandDimPos[1] = nDim - 2;
+
     // Write extra dim values - after projection for optimization.
     if (nDim > 2)
     {