Dem gridded_data shifted coordinates (#1682)

Author: pat-schmitt

docs/whats-new.rst

@@ -24,6 +24,10 @@ Bug fixes
   glacier mask when preserving the total values. This is a bad
   bug that is now fixed (:pull:`1661`).
   By `Patrick Schmitt <https://github.com/pat-schmitt>`_
+- When converting a variable of gridded_data to an tiff-file using
+  ``tasks.gridded_data_var_to_geotiff`` the resulting coordinates where
+  shifted half a pixel, this is now fixed (:pull:`1682`).
+  By `Patrick Schmitt <https://github.com/pat-schmitt>`_
 
 
 v1.6.1 (August 27, 2023)

oggm/core/gis.py

@@ -1678,7 +1678,7 @@ def gridded_data_var_to_geotiff(gdir, varname, fname=None):
         # Prepare the profile dict
         crs = ds.pyproj_srs
         var = ds[varname]
-        grid = ds.salem.grid
+        grid = ds.salem.grid.corner_grid
         data = var.data
         data_type = data.dtype.name
         height, width = var.data.shape

oggm/tests/test_prepro.py

@@ -13,6 +13,7 @@
 salem = pytest.importorskip('salem')
 rasterio = pytest.importorskip('rasterio')
 gpd = pytest.importorskip('geopandas')
+rioxr = pytest.importorskip('rioxarray')
 
 # Local imports
 import oggm
@@ -514,10 +515,18 @@ def test_gridded_data_var_to_geotiff(self):
 
         with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
             gridded_topo = ds[target_var]
-            gtiff_ds = salem.open_xr_dataset(gtiff_path)
-            assert ds.salem.grid == gtiff_ds.salem.grid
+            gtiff_ds = rioxr.open_rasterio(gtiff_path)
+            np.allclose(ds.salem.grid.x_coord, gtiff_ds.x)
+            np.allclose(ds.salem.grid.y_coord, gtiff_ds.y)
             assert np.allclose(gridded_topo.data, gtiff_ds.data)
 
+        # compare coordinates of topo.tif with dem.tif
+        demtiff_ds = rioxr.open_rasterio(gdir.get_filepath('dem'))
+        gtiff_ds = rioxr.open_rasterio(gtiff_path)
+        assert np.allclose(demtiff_ds.data, gtiff_ds.data)
+        assert np.allclose(demtiff_ds.x, gtiff_ds.x)
+        assert np.allclose(demtiff_ds.y, gtiff_ds.y)
+
 
 class TestCenterlines(unittest.TestCase):