Repository: pangeo-data/xESMF

Issue:
xesmf Regridder no longer keeps time information from original file.
Hello!

Let's get versions out of the way:

xESMF: 0.8.8
xarray: 2024.11.0
numpy: 2.0.2
python: 3.11.9

## Issue

I am using xesmf to regrid GFS data from the original GRIB2 files to the MRMS 0.01-degree grid. The time variables from the original GFS file are no longer carried through to the regridded file. Here is some sample code:

```
def build_archive(inFile, config, noclobber):

    print(f"Currently processing: {inFile}")
    # First, read/open the file.
    # This is going to be a little different. Due to the oddities of the GFS variables
    # that we're using for hail, the GFS files require some extra TLC.

    gfs_vars = config['variables']
    varlist = [v for v in gfs_vars.values()]
    datasets = []
    for k, v in gfs_vars.items():
        ds = xr.open_dataset(inFile, engine = 'cfgrib', filter_by_keys = {'cfVarName': k}, decode_timedelta=True)
        if k == 'unknown':
            ds = ds.rename_vars(name_dict={k:v})
        datasets.append(ds)

    dsGFS = xr.merge(datasets, compat = 'override').drop_vars('heightAboveGroundLayer')

    print(dsGFS.valid_time.values)
    dsMRMS = xr.open_dataset(config['mrms_file'], engine='cfgrib', decode_timedelta=True)

    # Get time information into datetime objects.
    initTime = dsGFS.time.values.astype('M8[ms]').astype('O')
    validTime = dsGFS.valid_time.values.astype('M8[ms]').astype('O')

    print(initTime)
    print(validTime)
    # Check to see if file currently exists AND noclobber is set.
    if outputPathFull.exists() and noclobber:

        dsGFS.close()
        dsMRMS.close()
        del dsGFS
        del dsMRMS

    else:

        # Read in Regridder file
        regridder = xe.Regridder(dsGFS, dsMRMS, method = 'bilinear', weights = config['regrid_file'], reuse_weights=True)

        # Regrid
        print("Regridding.")
        dsGFSonMRMS = regridder(dsGFS[varlist])

        edict = {'dtype' : 'f4', 'compression' : 'gzip', 'compression_opts' : 4, 'shuffle' : True}
        encoding = {var : edict for var in varlist}
        print(dsGFSonMRMS)
        print(dsGFSonMRMS.valid_time.data)
```

And here is the output from the print statements:

```
Currently processing: /hdd0/archive/hail/gfs/20230101/gfs_select.t00z.pgrb2.0p25.f000.grb2
2023-01-01T00:00:00.000000000
2023-01-01 00:00:00
2023-01-01 00:00:00
Regridding.
<xarray.Dataset> Size: 1GB
Dimensions:                  (latitude: 3500, longitude: 7000)
Coordinates:
    time                     datetime64[ns] 8B 2023-02-01T01:00:00
    step                     timedelta64[ns] 8B 00:00:00
    heightAboveGround        float64 8B 2.0
    valid_time               datetime64[ns] 8B 2023-02-01T01:00:00
    pressureFromGroundLayer  float64 8B 3e+03
    atmosphereSingleLayer    float64 8B 0.0
    surface                  float64 8B 0.0
    isobaricInhPa            float64 8B 500.0
    isothermZero             float64 8B 0.0
    heightAboveSea           float64 8B 500.0
  * latitude                 (latitude) float64 28kB 54.99 54.98 ... 20.02 20.01
  * longitude                (longitude) float64 56kB 230.0 230.0 ... 300.0
Data variables:
    t2m                      (latitude, longitude) float32 98MB 272.7 ... 298.7
    r                        (latitude, longitude) float32 98MB 99.49 ... 71.23
    pwat                     (latitude, longitude) float32 98MB 9.837 ... 33.28
    cape                     (latitude, longitude) float32 98MB 0.0 ... 1.097...
    cin                      (latitude, longitude) float32 98MB 0.4551 ... 0....
    w                        (latitude, longitude) float32 98MB 0.5072 ... 0....
    gh                       (latitude, longitude) float32 98MB 440.0 ... 4.3...
    hpbl                     (latitude, longitude) float32 98MB 31.67 ... 1.0...
    ustm                     (latitude, longitude) float32 98MB 4.766 ... -6.158
    vstm                     (latitude, longitude) float32 98MB -4.588 ... -9...
    hlcy                     (latitude, longitude) float32 98MB 122.0 ... 40.56
Attributes:
    regrid_method:  bilinear
2023-02-01T01:00:00.000000000
```

So as you can see, `valid_time` has now changed. That specific DTG, 2023-02-01T01:00, is the same time/valid_time from the MRMS file I'm using as the basis for the regrid-to grid. I can manually set the `time` and `valid_time` values from the original GFS file after the regridding via:

```
        dsGFSonMRMS['valid_time'] = dsGFS['valid_time']
        dsGFSonMRMS['time'] = dsGFS['time']
```

Is this expected behavior? I have done a similar remapping from GFS to the GOES geostationary grid, and the chain doesn't do this -- it'll keep the correct time information for the GFS. This chain uses older versions. My older versions are:

xESMF: 0.8.5
xarray: 2024.3.0
numpy: 1.26.4
python: 3.11.8

Discussion:

Patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 0f2ebd3..c4e6339 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -26,6 +26,7 @@ repos:
     rev: 7.1.1
     hooks:
       - id: flake8
+
   - repo: https://github.com/PyCQA/isort
     rev: 6.0.0
     hooks:
diff --git a/CHANGES.rst b/CHANGES.rst
index a55028e..2c82fdc 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -4,6 +4,7 @@ What's new
 0.8.9 (unreleased)
 ------------------
 * Destroy grids explicitly once weights are computed. Do not store them in `grid_in` and  `grid_out` attributes. This fixes segmentation faults introduced by the memory fix of last version. By `Pascal Bourgault <https://github.com/aulemahal>`_.
+* Do not add scalar coordinates of the target grid to the regridded output (:issue:`417`, :pull:`418`). `xe.Regridder.out_coords` is now a dataset instead of a dictionary. By `Pascal Bourgault <https://github.com/aulemahal>`_.
 
 0.8.8 (2024-11-01)
 ------------------
diff --git a/pyproject.toml b/pyproject.toml
index 98908b0..32e3251 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -61,7 +61,7 @@ tag_regex = "^(?P<prefix>v)?(?P<version>[^\\+]+)(?P<suffix>.*)?$"
 [tool.black]
 line-length = 100
 target-version = [
-  'py310',
+  'py311',
 ]
 skip-string-normalization = true
 
diff --git a/setup.cfg b/setup.cfg
index e6774ae..4e2cd68 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -7,4 +7,4 @@ ignore =
 max-line-length = 100
 max-complexity = 18
 select = B,C,E,F,W,T4,B9
-extend-ignore = E203,E501,E402,W605
+extend-ignore = E203,E501,E402,W503,W605
diff --git a/xesmf/frontend.py b/xesmf/frontend.py
index 989662a..a5bd611 100644
--- a/xesmf/frontend.py
+++ b/xesmf/frontend.py
@@ -962,18 +962,23 @@ class Regridder(BaseRegridder):
             self.out_horiz_dims = (lat_out.dims[0], lon_out.dims[0])
 
         if isinstance(ds_out, Dataset):
-            self.out_coords = {
-                name: crd
-                for name, crd in ds_out.coords.items()
-                if set(self.out_horiz_dims).issuperset(crd.dims)
-            }
+            out_coords = ds_out.coords.to_dataset()
             grid_mapping = {
                 var.attrs['grid_mapping']
                 for var in ds_out.data_vars.values()
                 if 'grid_mapping' in var.attrs
             }
-            if grid_mapping:
-                self.out_coords.update({gm: ds_out[gm] for gm in grid_mapping if gm in ds_out})
+            #  to keep : grid_mappings    and    non-scalar coords that have the spatial dims
+            self.out_coords = out_coords.drop_vars(
+                [
+                    name
+                    for name, crd in out_coords.coords.items()
+                    if not (
+                        (name in grid_mapping)
+                        or (len(crd.dims) > 0 and set(self.out_horiz_dims).issuperset(crd.dims))
+                    )
+                ]
+            )
         else:
             self.out_coords = {lat_out.name: lat_out, lon_out.name: lon_out}
 
@@ -1055,10 +1060,14 @@ class Regridder(BaseRegridder):
         chunks = out_chunks | in_chunks
 
         # Rename coords to avoid issues in xr.map_blocks
-        for coord in list(self.out_coords.keys()):
-            # If coords and dims are the same, renaming has already been done.
-            if coord not in self.out_horiz_dims:
-                ds_out = ds_out.rename({coord: coord + '_out'})
+        # If coords and dims are the same, renaming has already been done.
+        ds_out = ds_out.rename(
+            {
+                coord: coord + '_out'
+                for coord in self.out_coords.coords.keys()
+                if coord not in self.out_horiz_dims
+            }
+        )
 
         weights_dims = ('y_out', 'x_out', 'y_in', 'x_in')
         templ = sps.zeros((self.shape_out + self.shape_in))
@@ -1102,7 +1111,7 @@ class Regridder(BaseRegridder):
             # rename dimension name to match output grid
             out = out.rename({nd: od for nd, od in zip(new_dims, self.out_horiz_dims)})
 
-        out = out.assign_coords(**self.out_coords)
+        out = out.assign_coords(self.out_coords.coords)
         out.attrs['regrid_method'] = self.method
 
         if self.sequence_out:


Test patch:
diff --git a/xesmf/tests/test_frontend.py b/xesmf/tests/test_frontend.py
index 4a5a595..1a394a7 100644
--- a/xesmf/tests/test_frontend.py
+++ b/xesmf/tests/test_frontend.py
@@ -732,17 +732,22 @@ def test_regrid_dataset_extracoords():
         x=np.arange(24),
         y=np.arange(20),  # coords to be transfered
         latitude_longitude=xr.DataArray(),  # grid_mapping
-        bogus=ds_out.lev * ds_out.lon,  # coord not to be transfered
+        bogus=ds_out.lev * ds_out.lon,  # coords not to be transfered
+        scalar1=1,  #
+        scalar2=1,  #
     )
     ds_out2['data_ref'].attrs['grid_mapping'] = 'latitude_longitude'
     ds_out2['data4D_ref'].attrs['grid_mapping'] = 'latitude_longitude'
 
+    ds_in2 = ds_in.assign_coords(scalar2=5)
     regridder = xe.Regridder(ds_in, ds_out2, 'conservative')
-    ds_result = regridder(ds_in)
+    ds_result = regridder(ds_in2)
 
     assert 'x' in ds_result.coords
     assert 'y' in ds_result.coords
     assert 'bogus' not in ds_result.coords
+    assert 'scalar1' not in ds_result.coords
+    assert ds_result.scalar2 == 5
     assert 'latitude_longitude' in ds_result.coords