Merge pull request #577 from MetOffice/528_extract_single_point_data

Extract single point data

cset-workflow/includes/deterministic_single_point_surface_time_series.cylc

@@ -0,0 +1,16 @@
+{% if SURFACE_SINGLE_POINT_TIME_SERIES %}
+{% for model_field in SURFACE_MODEL_FIELDS %}
+[runtime]
+    [[pre_process_surface_single_point_time_series_{{model_field}}]]
+    inherit = PARALLEL
+        [[[environment]]]
+        CSET_RECIPE_NAME = "generic_surface_single_point_time_series.yaml"
+        CSET_ADDOPTS = "--VARNAME={{model_field}} --LONGITUDE_POINT={{LONGITUDE_POINT}} --LATITUDE_POINT={{LATITUDE_POINT}} --SINGLE_POINT_METHOD={{SINGLE_POINT_METHOD}}"
+
+    [[collate_surface_single_point_time_series_{{model_field}}]]
+    inherit = COLLATE
+        [[[environment]]]
+        CSET_RECIPE_NAME = "generic_surface_single_point_time_series.yaml"
+        CSET_ADDOPTS = "--VARNAME={{model_field}} --LONGITUDE_POINT={{LONGITUDE_POINT}} --LATITUDE_POINT={{LATITUDE_POINT}} --SINGLE_POINT_METHOD={{SINGLE_POINT_METHOD}}"
+{% endfor %}
+{% endif %}

cset-workflow/meta/rose-meta.conf

@@ -702,6 +702,41 @@ type=python_boolean
 compulsory=true
 sort-key=stash2
 
+[template variables=SURFACE_SINGLE_POINT_TIME_SERIES]
+ns=Diagnostics
+description=Plot a time series at a single specified location in a surface field.
+help=Include values of desired longitude and latitude.
+type=python_boolean
+compulsory=true
+trigger=template variables=LATITUDE_POINT: True;
+        template variables=LONGITUDE_POINT: True;
+        template variables=SINGLE_POINT_METHOD: True;
+sort-key=point1
+
+[template variables=LATITUDE_POINT]
+ns=Diagnostics
+description=Latitude of selected point. Note that this could be rotated or not, depending on the data provided.
+help=The latitude must exist within the domain. Value should be a float: for example, -1.5.
+type=real
+compulsory=true
+sort-key=point2
+
+[template variables=LONGITUDE_POINT]
+ns=Diagnostics
+description=Longitude of selected point. Note that this could be rotated or not, depending on the data provided.
+help=The longitude must exist within the domain. Value should be a float: for example, 0.8.
+type=real
+compulsory=true
+sort-key=point2
+
+[template variables=SINGLE_POINT_METHOD]
+ns=Diagnostics
+description=Method used to map model data onto selected gridpoints.
+help=Method used to map model data onto selected gridpoints. These are regrid methods available in Iris.
+values="Nearest", "Linear"
+compulsory=true
+sort-key=point2
+
 [template variables=BASIC_QQ_PLOT]
 ns=Diagnostics
 description=Create a basic quantile-quantile plot for variables specified collapsing over specified coordinates.

docs/source/changelog.rst

@@ -15,6 +15,9 @@ Unreleased
 .. Highlight any user facing changes. E.g:
 .. "* `@gh-user`_ did foo to bar in :pr:`9999`. This enables baz."
 
+* `@JKPShonk`_ and `@cehalliwell`_ added functionality to CSET to allow it to
+   generate time series plots from model data mapped on to a selected
+   longitude/latitude location in :pr:`577`
 * `@Sylviabohnenstengel`_ add pdf functionality for spatial fields in :pr:`594`.
 * `@Sylviabohnenstengel`_ documentation: add info on quick pytesting in :pr:`696`
 * `@Sylviabohnenstengel`_ add constraint operator for lfric full_levels and half_levels
@@ -105,6 +108,9 @@ to come.
 * `@jfrost-mo`_ updated the example rose-suite.conf to reflect what a modern
   version should look like in :pr:`508`
 
+.. _@JKPShonk: https://github.com/JKPShonk
+.. _@cehalliwell: https://github.com/cehalliwell
+
 .. _CSET discussion forum: https://github.com/MetOffice/simulation-systems/discussions/categories/cset-toolkit
 .. _@dasha-shchep: https://github.com/dasha-shchep
 .. _@cjohnson-pi: https://github.com/cjohnson-pi

src/CSET/operators/regrid.py

@@ -14,6 +14,8 @@
 
 """Operators to regrid cubes."""
 
+import warnings
+
 import iris
 import iris.cube
 import numpy as np
@@ -22,6 +24,14 @@
 from CSET.operators._utils import get_cube_yxcoordname
 
 
+class BoundaryWarning(UserWarning):
+    """Selected gridpoint is close to the domain edge.
+
+    In many cases gridpoints near the domain boundary contain non-physical
+    values, so caution is advised when interpreting them.
+    """
+
+
 def regrid_onto_cube(
     toregrid: iris.cube.Cube | iris.cube.CubeList,
     target: iris.cube.Cube,
@@ -184,3 +194,108 @@ def regrid_onto_xyspacing(
         return regridded_cubes[0]
     else:
         return regridded_cubes
+
+
+def regrid_to_single_point(
+    incube: iris.cube.Cube, lat_pt: float, lon_pt: float, method: str, **kwargs
+) -> iris.cube.Cube:
+    """Select data at a single point by longitude and latitude.
+
+    Selection of model grid point is performed by a regrid function, either selecting the
+    nearest gridpoint to the selected longitude and latitude values or using linear
+    interpolation across the surrounding points.
+
+    Parameters
+    ----------
+    incube: Cube
+        An iris cube of the data to regrid. As a minimum, it needs to be 2D with latitude,
+        longitude coordinates.
+    lon_pt: float
+        Selected value of longitude.
+    lat_pt: float
+        Selected value of latitude.
+    method: str
+        Method used to determine the values at the selected longitude and latitude.
+        The recommended approach is to use iris.analysis.Nearest(), which selects the
+        nearest gridpoint. An alternative is iris.analysis.Linear(), which obtains
+        the values at the selected longitude and latitude by linear interpolation.
+
+    Returns
+    -------
+    cube_rgd: Cube
+        An iris cube of the data at the specified point (this may have time
+        and/or height dimensions).
+
+    Raises
+    ------
+    ValueError
+        If a unique x/y coordinate cannot be found; also if, for selecting a single
+        gridpoint, the chosen longitude and latitude point is outside the domain.
+    NotImplementedError
+        If the cubes grid, or the method for regridding, is not yet supported.
+
+    Notes
+    -----
+    The acceptable coordinate names for X and Y coordinates are currently described
+    in X_COORD_NAMES and Y_COORD_NAMES. These cover commonly used coordinate types,
+    though a user can append new ones.
+    Currently rectilinear grids (uniform) are supported.
+    Warnings are raised if the selected gridpoint is within eight gridlengths of the
+    domain boundary as data here is potentially unreliable.
+
+    """
+    # Get x and y coordinate names.
+    y_coord, x_coord = get_cube_yxcoordname(incube)
+
+    # List of supported grids - check if it is compatible
+    # NOTE: The "RotatedGeogCS" option below seems to be required for rotated grids --
+    #  this may need to be added in other places in these Operators.
+    supported_grids = (iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS)
+    if not isinstance(incube.coord(x_coord).coord_system, supported_grids):
+        raise NotImplementedError(
+            f"Does not currently support {incube.coord(x_coord).coord_system} regrid method"
+        )
+    if not isinstance(incube.coord(y_coord).coord_system, supported_grids):
+        raise NotImplementedError(
+            f"Does not currently support {incube.coord(y_coord).coord_system} regrid method"
+        )
+
+    # Get axis
+    lat, lon = incube.coord(y_coord), incube.coord(x_coord)
+
+    # Get bounds
+    lat_min, lon_min = lat.points.min(), lon.points.min()
+    lat_max, lon_max = lat.points.max(), lon.points.max()
+
+    # Get bounds
+    # Boundaries of frame to avoid selecting gridpoint close to domain edge
+    # Currently hardwired to 8 but could be a user input
+    lat_min_bound, lon_min_bound = lat.points[7], lon.points[7]
+    lat_max_bound, lon_max_bound = lat.points[-8], lon.points[-8]
+
+    # Check to see if selected point is outside the domain
+    if (
+        (lat_pt < lat_min)
+        or (lat_pt > lat_max)
+        or (lon_pt < lon_min)
+        or (lon_pt > lon_max)
+    ):
+        raise ValueError("Selected point is outside the domain.")
+    else:
+        if (
+            (lat_pt < lat_min_bound)
+            or (lat_pt > lat_max_bound)
+            or (lon_pt < lon_min_bound)
+            or (lon_pt > lon_max_bound)
+        ):
+            warnings.warn(
+                "Selected point is within 8 gridlengths of the domain edge.",
+                category=BoundaryWarning,
+                stacklevel=2,
+            )
+
+    regrid_method = getattr(iris.analysis, method, None)
+    if not callable(regrid_method):
+        raise NotImplementedError(f"Does not currently support {method} regrid method")
+    cube_rgd = incube.interpolate(((lat, lat_pt), (lon, lon_pt)), regrid_method())
+    return cube_rgd

src/CSET/recipes/generic_surface_single_point_time_series.yaml

@@ -0,0 +1,41 @@
+category: Quick Look
+title: Time series of $VARNAME at $LATITUDE_POINT N, $LONGITUDE_POINT E
+description: Plots a time series of the surface $VARNAME at a selected gridpoint.
+
+# Parallel steps.
+parallel:
+  - operator: read.read_cube
+    constraint:
+      operator: constraints.combine_constraints
+      varname_constraint:
+        operator: constraints.generate_var_constraint
+        varname: $VARNAME
+      pressure_level_constraint:
+        operator: constraints.generate_level_constraint
+        levels: []
+        coordinate: "pressure"
+      validity_time_constraint:
+        operator: constraints.generate_time_constraint
+        time_start: $VALIDITY_TIME
+
+  - operator: regrid.regrid_to_single_point
+    lat_pt: $LATITUDE_POINT
+    lon_pt: $LONGITUDE_POINT
+    method: $SINGLE_POINT_METHOD
+
+  # Save single-point variable to a file per validity time.
+  - operator: write.write_cube_to_nc
+    filename: intermediate/single_point_values
+
+# Collation steps.
+# Reads in intermediate cube and plots it.
+collate:
+  - operator: read.read_cube
+    filename_pattern: intermediate/*.nc
+
+  # Make a single NetCDF with all the data inside it.
+  - operator: write.write_cube_to_nc
+    overwrite: True
+
+  # Plot the data.
+  - operator: plot.plot_line_series

tests/operators/test_regrid.py

@@ -21,6 +21,7 @@
 import pytest
 
 import CSET.operators.regrid as regrid
+from CSET.operators.regrid import BoundaryWarning
 
 
 # Session scope fixtures, so the test data only has to be loaded once.
@@ -148,3 +149,70 @@ def test_regrid_onto_xyspacing_unknown_method(regrid_source_cube):
         regrid.regrid_onto_xyspacing(
             regrid_source_cube, xspacing=0.5, yspacing=0.5, method="nonexistent"
         )
+
+
+@pytest.mark.filterwarnings("ignore:Selected point is within")
+def test_regrid_to_single_point(cube):
+    """Regrid to single point."""
+    # Test extracting a single point.
+    regrid_cube = regrid.regrid_to_single_point(cube, 0.5, 358.5, "Nearest")
+    expected_cube = "<iris 'Cube' of air_temperature / (K) (time: 3)>"
+    assert repr(regrid_cube) == expected_cube
+
+
+@pytest.mark.filterwarnings("ignore:Selected point is within")
+def test_regrid_to_single_point_missing_coord(cube):
+    """Missing coordinate raises error."""
+    # Missing X coordinate.
+    source = cube.copy()
+    source.remove_coord("grid_longitude")
+    with pytest.raises(ValueError):
+        regrid.regrid_to_single_point(source, 0.5, 358.5, "Nearest")
+
+    # Missing Y coordinate.
+    source = cube.copy()
+    source.remove_coord("grid_latitude")
+    with pytest.raises(ValueError):
+        regrid.regrid_to_single_point(source, 0.5, 358.5, "Nearest")
+
+
+def test_regrid_to_single_point_unknown_crs_x(cube):
+    """X coordinate reference system is unrecognised."""
+    # Exchange to unsupported coordinate system.
+    cube.coord("grid_longitude").coord_system = iris.coord_systems.OSGB()
+    with pytest.raises(NotImplementedError):
+        regrid.regrid_to_single_point(cube, 0.5, 358.5, "Nearest")
+
+
+def test_regrid_to_single_point_unknown_crs_y(cube):
+    """Y coordinate reference system is unrecognised."""
+    # Exchange to unsupported coordinate system.
+    cube.coord("grid_latitude").coord_system = iris.coord_systems.OSGB()
+    with pytest.raises(NotImplementedError):
+        regrid.regrid_to_single_point(cube, 0.5, 358.5, "Nearest")
+
+
+def test_regrid_to_single_point_outside_domain(regrid_source_cube):
+    """Error if coordinates are outside the model domain."""
+    with pytest.raises(ValueError):
+        regrid.regrid_to_single_point(regrid_source_cube, 0.5, 178.5, "Nearest")
+
+
+def test_regrid_to_single_point_unknown_method(cube):
+    """Method does not exist."""
+    with pytest.raises(NotImplementedError):
+        regrid.regrid_to_single_point(cube, 0.5, 358.5, method="nonexistent")
+
+
+@pytest.mark.filterwarnings(
+    "ignore:this date/calendar/year zero convention is not supported by CF"
+)
+def test_boundary_warning(regrid_source_cube):
+    """Ensures a warning is raised when chosen point is too close to boundary."""
+    with pytest.warns(
+        BoundaryWarning, match="Selected point is within 8 gridlengths"
+    ) as warning:
+        regrid.regrid_to_single_point(regrid_source_cube, -0.9, 393.0, "Nearest")
+
+    assert len(warning) == 1
+    assert issubclass(warning[-1].category, BoundaryWarning)