Fix units and long names of derived ERA5 variables

RELEASE_NOTES.rst

@@ -46,6 +46,11 @@ Upcoming Release
 
 **Bug fixes**
 
+* Fix the ``units`` and ``long_name`` attributes of the ERA5 variables ``height``
+  (was geopotential ``m**2 s**-2``), ``wnd_azimuth`` (was the ``u100`` wind
+  component) and ``solar_altitude``/``solar_azimuth`` (inherited ``latitude``).
+  Only metadata changes, values are unchanged (`#509 <https://github.com/PyPSA/atlite/issues/509>`_).
+
 * Fix ``get_oedb_windturbineconfig`` applying the documented ``turbine_type``
   search parameter to the value of ``name``. Searching by ``turbine_type``
   alone raised ``KeyError: 'name'``, and combining it with ``name`` silently

atlite/datasets/era5.py

@@ -98,7 +98,7 @@ def _add_height(ds: xr.Dataset) -> xr.Dataset:
     z = ds["z"]
     if "time" in z.coords:
         z = z.isel(time=0, drop=True)
-    ds["height"] = z / g0
+    ds["height"] = (z / g0).assign_attrs(units="m", long_name="Height")
     return ds.drop_vars("z")
 
 
@@ -159,7 +159,9 @@ def _process_wind(ds: xr.Dataset, single_precision: bool = False) -> xr.Dataset:
 
     # span the whole circle: 0 is north, π/2 is east, -π is south, 3π/2 is west
     azimuth = np.arctan2(ds["u100"], ds["v100"])
-    azimuth = azimuth.where(azimuth >= 0, azimuth + 2 * np.pi)
+    azimuth = azimuth.where(azimuth >= 0, azimuth + 2 * np.pi).assign_attrs(
+        units="rad", long_name="100 metre wind azimuth"
+    )
     ds["wnd_azimuth"] = azimuth.astype(np.float32) if single_precision else azimuth
 
     ds = ds.drop_vars(["u100", "v100", "u10", "v10", "wnd10m"])

atlite/pv/solar_position.py

@@ -114,17 +114,23 @@ def SolarPosition(ds: xr.Dataset, time_shift: str | pd.Timedelta = "0H") -> xr.D
     alt = arcsin(
         (sin(dec) * sin(lat) + cos(dec) * cos(lat) * cos(h)).clip(min=-1.0, max=1.0)
     ).rename("altitude")
-    alt.attrs["time shift"] = f"{time_shift}"
-    alt.attrs["units"] = "rad"
+    alt.attrs = {
+        "time shift": f"{time_shift}",
+        "units": "rad",
+        "long_name": "solar altitude",
+    }
 
     az = arccos(
         ((sin(dec) * cos(lat) - cos(dec) * sin(lat) * cos(h)) / cos(alt)).clip(
             min=-1.0, max=1.0
         )
     )
     az = az.where(h <= 0, 2 * pi - az).rename("azimuth")
-    az.attrs["time shift"] = f"{time_shift}"
-    az.attrs["units"] = "rad"
+    az.attrs = {
+        "time shift": f"{time_shift}",
+        "units": "rad",
+        "long_name": "solar azimuth",
+    }
 
     vars = {da.name: da for da in [alt, az]}
     return xr.Dataset(vars)

test/test_era5_attrs.py

@@ -0,0 +1,62 @@
+# SPDX-FileCopyrightText: Contributors to atlite <https://github.com/pypsa/atlite>
+#
+# SPDX-License-Identifier: MIT
+
+"""Tests for the units and names of derived ERA5 variables (#509)."""
+
+import numpy as np
+import pandas as pd
+import xarray as xr
+
+from atlite.datasets.era5 import _add_height, _process_influx, _process_wind
+
+
+def _field(value: float, units: str, long_name: str) -> xr.DataArray:
+    return xr.DataArray(
+        np.full((2, 2, 2), value),
+        coords={
+            "time": pd.date_range("2013-01-01", periods=2, freq="h"),
+            "y": [56.0, 56.25],
+            "x": [0.0, 0.25],
+        },
+        dims=("time", "y", "x"),
+        attrs={"units": units, "long_name": long_name},
+    )
+
+
+def test_height_attrs():
+    ds = _add_height(xr.Dataset({"z": _field(100.0, "m**2 s**-2", "Geopotential")}))
+    assert ds["height"].attrs == {"units": "m", "long_name": "Height"}
+
+
+def test_wind_azimuth_attrs():
+    raw = xr.Dataset({
+        name: _field(value, "m s**-1", f"{name} wind component")
+        for name, value in [("u10", 1.0), ("v10", 2.0), ("u100", 3.0), ("v100", 4.0)]
+    })
+    raw["fsr"] = _field(0.1, "m", "Forecast surface roughness")
+    ds = _process_wind(raw)
+    assert ds["wnd_azimuth"].attrs == {
+        "units": "rad",
+        "long_name": "100 metre wind azimuth",
+    }
+
+
+def test_solar_position_attrs():
+    raw = xr.Dataset({
+        name: _field(value, "J m**-2", name)
+        for name, value in [
+            ("ssrd", 3600.0),
+            ("ssr", 1800.0),
+            ("fdir", 1800.0),
+            ("tisr", 7200.0),
+        ]
+    })
+    raw = raw.assign_coords(
+        lon=raw.x.assign_attrs(long_name="longitude"),
+        lat=raw.y.assign_attrs(units="degrees_north", long_name="latitude"),
+    )
+    ds = _process_influx(raw)
+    assert ds["solar_altitude"].attrs["long_name"] == "solar altitude"
+    assert ds["solar_azimuth"].attrs["long_name"] == "solar azimuth"
+    assert ds["solar_altitude"].attrs["units"] == "rad"