Implement relative altitude filter

scripts/make_topography_subset.py

@@ -0,0 +1,14 @@
+import xarray as xr
+
+# Script which extracts only the first time point and topography value from the emep topography
+# file to reduce data for the tests.
+data = xr.open_dataset("/lustre/storeB/project/fou/kl/emep/Auxiliary/topography.nc")
+
+start_time = data["time"][0]
+
+data = data.sel(time=slice(start_time))
+
+data = data["topography"]
+
+
+data.to_netcdf("tests/testdata/datadir_elevation/topography.nc")

setup.cfg

@@ -25,6 +25,8 @@ python_version = >=3.10
 install_requires =
     numpy >= 1.13
     importlib-metadata >= 3.6; python_version < "3.10"
+    xarray
+    netcdf4
 package_dir =
     =src
 packages = pyaro, pyaro.timeseries, pyaro.csvreader

src/pyaro/timeseries/AutoFilterReaderEngine.py

@@ -31,7 +31,7 @@ def supported_filters(cls) -> list[Filter]:
 
         :return: list of filters
         """
-        supported = "variables,stations,countries,bounding_boxes,duplicates,time_bounds,time_resolution,flags,time_variable_station,altitude".split(
+        supported = "variables,stations,countries,bounding_boxes,duplicates,time_bounds,time_resolution,flags,time_variable_station,altitude,relaltitude".split(
             ","
         )
         return [filters.get(name) for name in supported]

src/pyaro/timeseries/Filter.py

@@ -1,3 +1,4 @@
+import logging
 import math
 import abc
 from collections import defaultdict
@@ -9,10 +10,12 @@
 from typing import Any
 
 import numpy as np
+import xarray as xr
 
 from .Data import Data, Flag
 from .Station import Station
 
+logger = logging.getLogger(__name__)
 
 class Filter(abc.ABC):
     """Base-class for all filters used from pyaro-Readers"""
@@ -855,4 +858,83 @@ def filter_stations(self, stations: dict[str, Station]) -> dict[str, Station]:
         if self._max_altitude is not None:
             stations = {n: s for n, s in stations.items() if (not math.isnan(s["altitude"]) and s["altitude"] <= self._max_altitude) }
 
-        return stations
+        return stations
+
+@registered_filter
+class RelativeAltitudeFilter(StationFilter):
+    """
+    Filter class which filters stations based on the relative difference between
+    the station altitude, and the model topography altitude.
+
+    https://github.com/metno/pyaro/issues/39
+    """
+    def __init__(self, topo_file: str | None = None, topo_var: str = "topography", rtol: float = 1):
+        """
+        :param topo_file : A .nc file from which to read model topography data.
+        :param topo_var : Name of variable that stores altitude.
+        :param rtol : Relative toleranse.
+
+        Note:
+        -----
+        Stations will be kept if abs(altobs-altmod) <= rtol*abs(altobs)
+        """
+        self._topo_file = topo_file
+        self._topo_var = topo_var
+        self._rtol = rtol
+
+        self._topography = None
+        if topo_file is not None:
+            self._topography = xr.open_dataset(topo_file)
+        else:
+            logger.warning("No topography data provided (topo_file='%s'). Relative elevation filtering will not be applied.", topo_file)
+
+    def _model_altitude_from_lat_lon(self, lat: float, lon: float) -> float:
+        # TODO: Include a tolerance?
+        data = self._topography.sel({"lat": lat, "lon": lon}, method="nearest")
+
+        # Should not vary in time too much so picking the first one here.
+        altitude = data["topography"][0]
+
+        return float(altitude)
+
+    def _is_close(self, altmod: float, altobs: float) -> bool:
+        """
+        Function to check if two altitudes are within a relative tolerance of each
+        other.
+
+        :param altmod : Model altitude (in meters).
+        :param altobs : Observation / station altitude (in meters).
+
+        :returns :
+            True if the values are close with station altitude as the reference 
+            value.
+        """
+        return abs(altmod-altobs) <= (self._rtol * abs(altobs))
+
+    def init_kwargs(self):
+        return {
+            "topo_file": self._file,
+            "topo_var": self._topo_var,
+            "rtol": self._rtol 
+        }
+
+    def name(self):
+        return "relaltitude"
+
+    def filter_stations(self, stations: dict[str, Station]) -> dict[str, Station]:
+        if self._topography is None:
+            return stations
+
+        new_stations = dict()
+
+        for name, station in stations.items():
+            lat = station["latitude"]
+            lon = station["longitude"]
+
+            altobs = station["altitude"]
+            altmod = self._model_altitude_from_lat_lon(lat, lon)
+
+            if self._is_close(altmod, altobs):
+                new_stations[name] = station
+
+        return new_stations

tests/test_CSVTimeSeriesReader.py

@@ -485,6 +485,82 @@ def test_altitude_filter_3(self):
             }
         ) as ts:
             self.assertEqual(len(ts.stations()), 1)
+
+    def test_relaltitude_filter_1(self):
+        engines = pyaro.list_timeseries_engines()
+        with engines["csv_timeseries"].open(
+            filename=self.elevation_file,
+            filters=[pyaro.timeseries.filters.get("relaltitude", topo_file = "./tests/testdata/datadir_elevation/topography.nc", rtol=0)],
+            columns={
+                "variable": 0,
+                "station": 1,
+                "longitude": 2,
+                "latitude": 3,
+                "value": 4,
+                "units": 5,
+                "start_time": 6,
+                "end_time": 7,
+                "altitude": 9,
+                "country": "NO",
+                "standard_deviation": "NaN",
+                "flag": "0",
+            }
+        ) as ts:
+            # Altitudes in dataset:
+            # Station     | Alt_obs   | Modeobs |
+            # Station 1   | 100       | 12.2554 |
+            # Station 2   | 200       |  4.9016 |
+            # Station 3   | 300       |  4.9016 |
+            # Since rtol = 0, no station should be included.
+            self.assertEqual(len(ts.stations()), 0)
+
+    def test_relaltitude_filter_2(self):
+        engines = pyaro.list_timeseries_engines()
+        with engines["csv_timeseries"].open(
+            filename=self.elevation_file,
+            filters=[pyaro.timeseries.filters.get("relaltitude", topo_file = "./tests/testdata/datadir_elevation/topography.nc", rtol=0.89)],
+            columns={
+                "variable": 0,
+                "station": 1,
+                "longitude": 2,
+                "latitude": 3,
+                "value": 4,
+                "units": 5,
+                "start_time": 6,
+                "end_time": 7,
+                "altitude": 9,
+                "country": "NO",
+                "standard_deviation": "NaN",
+                "flag": "0",
+            }
+        ) as ts:
+            # At rtol = 0.89, only the first station should be included.
+            self.assertEqual(len(ts.stations()), 1)
+
+    def test_relaltitude_filter_3(self):
+        engines = pyaro.list_timeseries_engines()
+        with engines["csv_timeseries"].open(
+            filename=self.elevation_file,
+            filters=[pyaro.timeseries.filters.get("relaltitude", topo_file = "./tests/testdata/datadir_elevation/topography.nc", rtol=1)],
+            columns={
+                "variable": 0,
+                "station": 1,
+                "longitude": 2,
+                "latitude": 3,
+                "value": 4,
+                "units": 5,
+                "start_time": 6,
+                "end_time": 7,
+                "altitude": 9,
+                "country": "NO",
+                "standard_deviation": "NaN",
+                "flag": "0",
+            }
+        ) as ts:
+            # Since rtol=1, all stations should be included.
+            self.assertEqual(len(ts.stations()), 3)
+
+
 
 
 if __name__ == "__main__":