Allow real-world co-ordinates to be specified in single-point timeseries.

cset-workflow/includes/deterministic_single_point_surface_time_series.cylc

@@ -11,6 +11,7 @@
             --MODEL_NAME='{{model["name"]}}'
             --LONGITUDE_POINT='{{LONGITUDE_POINT}}'
             --LATITUDE_POINT='{{LATITUDE_POINT}}'
+            --LATLON_IN_TYPE='{{LATLON_IN_TYPE}}'
             --SINGLE_POINT_METHOD='{{SINGLE_POINT_METHOD}}'
         """
         MODEL_IDENTIFIERS = {{model["id"]}}

cset-workflow/meta/diagnostics/rose-meta.conf

@@ -137,25 +137,35 @@ type=python_boolean
 compulsory=true
 trigger=template variables=LATITUDE_POINT: True;
         template variables=LONGITUDE_POINT: True;
+        template variables=LATLON_IN_TYPE: True;
         template variables=SINGLE_POINT_METHOD: True;
 sort-key=0surface6
 
 [template variables=LATITUDE_POINT]
 ns=Diagnostics/Quicklook
-description=Latitude of selected point in the same coordinate system as the data.
+description=Latitude of selected point.
 help=The latitude must exist within the domain. Value should be a float: for example, -1.5.
 type=real
 compulsory=true
 sort-key=0surface7
 
 [template variables=LONGITUDE_POINT]
 ns=Diagnostics/Quicklook
-description=Longitude of selected point in the same coordinate system as the data.
+description=Longitude of selected point.
 help=The longitude must exist within the domain. Value should be a float: for example, 0.8.
 type=real
 compulsory=true
 sort-key=0surface8
 
+[template variables=LATLON_IN_TYPE]
+ns=Diagnostics/Quicklook
+description=Method used to map model data onto selected gridpoints.
+help=This switch indicates whether the selected latitude longitude coordinates are on the real world grid or
+    the rotated grid specified by the input data.
+values="realworld", "rotated"
+compulsory=true
+sort-key=0surface6
+
 [template variables=SINGLE_POINT_METHOD]
 ns=Diagnostics/Quicklook
 description=Method used to map model data onto selected gridpoints.

src/CSET/operators/regrid.py

@@ -202,7 +202,8 @@ def regrid_to_single_point(
     cube: iris.cube.Cube,
     lat_pt: float,
     lon_pt: float,
-    method: str,
+    latlon_in_type: str = "rotated",
+    method: str = "Nearest",
     boundary_margin: int = 8,
     **kwargs,
 ) -> iris.cube.Cube:
@@ -274,6 +275,12 @@ def regrid_to_single_point(
             f"Does not currently support {cube.coord(y_coord).coord_system} regrid method"
         )
 
+    # Transform input coordinates onto rotated grid if requested
+    if latlon_in_type == "realworld":
+        lon_tr, lat_tr = transform_lat_long_points(lon_pt, lat_pt, cube)
+    elif latlon_in_type == "rotated":
+        lon_tr, lat_tr = lon_pt, lat_pt
+
     # Get axis
     lat, lon = cube.coord(y_coord), cube.coord(x_coord)
 
@@ -292,23 +299,23 @@ def regrid_to_single_point(
     )
 
     # Check to see if selected point is outside the domain
-    if (lat_pt < lat_min) or (lat_pt > lat_max):
+    if (lat_tr < lat_min) or (lat_tr > lat_max):
         raise ValueError("Selected point is outside the domain.")
     else:
-        if (lon_pt < lon_min) or (lon_pt > lon_max):
-            if (lon_pt + 360.0 >= lon_min) and (lon_pt + 360.0 <= lon_max):
-                lon_pt += 360.0
-            elif (lon_pt - 360.0 >= lon_min) and (lon_pt - 360.0 <= lon_max):
-                lon_pt -= 360.0
+        if (lon_tr < lon_min) or (lon_tr > lon_max):
+            if (lon_tr + 360.0 >= lon_min) and (lon_tr + 360.0 <= lon_max):
+                lon_tr += 360.0
+            elif (lon_tr - 360.0 >= lon_min) and (lon_tr - 360.0 <= lon_max):
+                lon_tr -= 360.0
             else:
                 raise ValueError("Selected point is outside the domain.")
 
     # Check to see if selected point is near the domain boundaries
     if (
-        (lat_pt < lat_min_bound)
-        or (lat_pt > lat_max_bound)
-        or (lon_pt < lon_min_bound)
-        or (lon_pt > lon_max_bound)
+        (lat_tr < lat_min_bound)
+        or (lat_tr > lat_max_bound)
+        or (lon_tr < lon_min_bound)
+        or (lon_tr > lon_max_bound)
     ):
         warnings.warn(
             f"Selected point is within {boundary_margin} gridlengths of the domain edge, data may be unreliable.",
@@ -319,5 +326,41 @@ def regrid_to_single_point(
     regrid_method = getattr(iris.analysis, method, None)
     if not callable(regrid_method):
         raise NotImplementedError(f"Does not currently support {method} regrid method")
-    cube_rgd = cube.interpolate(((lat, lat_pt), (lon, lon_pt)), regrid_method())
+    cube_rgd = cube.interpolate(((lat, lat_tr), (lon, lon_tr)), regrid_method())
     return cube_rgd
+
+
+def transform_lat_long_points(lon, lat, cube):
+    """Transform a selected point in longitude and latitude.
+
+    Transform the coordinates of a point from the real world
+    grid to the corresponding point on the rotated grid of a cube.
+
+    Parameters
+    ----------
+    cube: Cube
+        An iris cube of data defining the rotated grid to be used in
+        the longitude-latitude transformation.
+    lon: float
+        Selected value of longitude: this should be in the range -180 degrees to
+        180 degrees.
+    lat: float
+        Selected value of latitude: this should be in the range -90 degrees to
+        90 degrees.
+
+    Returns
+    -------
+    lon_rot, lat_rot: float
+        Coordinates of the selected point on the rotated grid specified within
+        the selected cube.
+
+    """
+    import cartopy.crs as ccrs
+
+    rot_pole = cube.coord_system().as_cartopy_crs()
+    true_grid = ccrs.Geodetic()
+    rot_coords = rot_pole.transform_point(lon, lat, true_grid)
+    lon_rot = rot_coords[0]
+    lat_rot = rot_coords[1]
+
+    return lon_rot, lat_rot

src/CSET/recipes/generic_surface_single_point_time_series.yaml

@@ -1,5 +1,5 @@
 category: Time Series
-title: $MODEL_NAME Time series of $VARNAME at $LATITUDE_POINT N, $LONGITUDE_POINT E
+title: $MODEL_NAME Time series of $VARNAME at $LATITUDE_POINT N, $LONGITUDE_POINT E ($LATLON_IN_TYPE)
 description: Plots a time series of the surface $VARNAME at a selected gridpoint.
 
 steps:
@@ -17,6 +17,7 @@ steps:
   - operator: regrid.regrid_to_single_point
     lat_pt: $LATITUDE_POINT
     lon_pt: $LONGITUDE_POINT
+    latlon_in_type: $LATLON_IN_TYPE
     method: $SINGLE_POINT_METHOD
 
   # Make a single NetCDF with all the data inside it.

tests/operators/test_regrid.py

@@ -178,7 +178,7 @@ def test_regrid_to_single_point_east(cube):
     """
     cube_fix = read._longitude_fix_callback(cube, None, None)
     regrid_cube = regrid.regrid_to_single_point(
-        cube_fix, 0.5, -1.5, "Nearest", boundary_margin=1
+        cube_fix, 0.5, -1.5, "rotated", "Nearest", boundary_margin=1
     )
     lat_name, lon_name = _utils.get_cube_yxcoordname(regrid_cube)
     lat_coord = regrid_cube.coord(lat_name)
@@ -201,7 +201,7 @@ def test_regrid_to_single_point_west(cube):
     cube.coord("grid_longitude").points = long_coord
     cube_fix = read._longitude_fix_callback(cube, None, None)
     regrid_cube = regrid.regrid_to_single_point(
-        cube_fix, 0.5, -1.5, "Nearest", boundary_margin=1
+        cube_fix, 0.5, -1.5, "rotated", "Nearest", boundary_margin=1
     )
     lat_name, lon_name = _utils.get_cube_yxcoordname(regrid_cube)
     lat_coord = regrid_cube.coord(lat_name)
@@ -220,7 +220,7 @@ def test_regrid_to_single_point_longitude_transform_1(cube):
     back into the standard range (-180 deg to 180 deg).
     """
     regrid_cube = regrid.regrid_to_single_point(
-        cube, 0.5, 358.5, "Nearest", boundary_margin=1
+        cube, 0.5, 358.5, "rotated", "Nearest", boundary_margin=1
     )
     lat_name, lon_name = _utils.get_cube_yxcoordname(regrid_cube)
     lat_coord = regrid_cube.coord(lat_name)
@@ -239,7 +239,7 @@ def test_regrid_to_single_point_longitude_transform_2(cube):
     back into the standard range (-180 deg to 180 deg).
     """
     regrid_cube = regrid.regrid_to_single_point(
-        cube, 0.5, -361.5, "Nearest", boundary_margin=1
+        cube, 0.5, -361.5, "rotated", "Nearest", boundary_margin=1
     )
     lat_name, lon_name = _utils.get_cube_yxcoordname(regrid_cube)
     lat_coord = regrid_cube.coord(lat_name)
@@ -250,19 +250,49 @@ def test_regrid_to_single_point_longitude_transform_2(cube):
     assert lon_coord.points[0] == -1.5
 
 
+def test_regrid_to_single_point_realworld(cube):
+    """Test extracting a single point.
+
+    Test that, if a real world coordinate is specified, the code is
+    mapping this point correctly onto the rotated grid.
+    """
+    cube_fix = read._longitude_fix_callback(cube, None, None)
+    regrid_cube = regrid.regrid_to_single_point(
+        cube_fix, 52.98, -5.0, "realworld", "Nearest", boundary_margin=1
+    )
+    expected_array = "array(288.59375, dtype=float32)"
+    assert repr(regrid_cube[0].data) == expected_array
+
+
+def test_regrid_to_single_point_rotated(cube):
+    """Test extracting a single point.
+
+    Test that, if a rotated coordinate is specified, the answer
+    matches the corresponding coordinate on the real world grid.
+    """
+    cube_fix = read._longitude_fix_callback(cube, None, None)
+    regrid_cube = regrid.regrid_to_single_point(
+        cube_fix, 0.5, -1.5, "rotated", "Nearest", boundary_margin=1
+    )
+    expected_array = "array(288.59375, dtype=float32)"
+    assert repr(regrid_cube[0].data) == expected_array
+
+
 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", boundary_margin=1)
+        regrid.regrid_to_single_point(
+            source, 0.5, 358.5, "rotated", "Nearest", boundary_margin=1
+        )
 
     # 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")
+        regrid.regrid_to_single_point(source, 0.5, 358.5, "rotated", "Nearest")
 
 
 def test_longitude_fix_callback_missing_coord(cube):
@@ -285,34 +315,38 @@ def test_regrid_to_single_point_unknown_crs_x(cube):
     # 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")
+        regrid.regrid_to_single_point(cube, 0.5, 358.5, "rotated", "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")
+        regrid.regrid_to_single_point(cube, 0.5, 358.5, "rotated", "Nearest")
 
 
 def test_regrid_to_single_point_outside_domain_longitude(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")
+        regrid.regrid_to_single_point(
+            regrid_source_cube, 0.5, 178.5, "rotated", "Nearest"
+        )
 
 
 def test_regrid_to_single_point_outside_domain_latitude(regrid_source_cube):
     """Error if coordinates are outside the model domain."""
     with pytest.raises(ValueError):
-        regrid.regrid_to_single_point(regrid_source_cube, 80.5, 358.5, "Nearest")
+        regrid.regrid_to_single_point(
+            regrid_source_cube, 80.5, 358.5, "rotated", "Nearest"
+        )
 
 
 @pytest.mark.filterwarnings("ignore:Selected point is within")
 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")
+        regrid.regrid_to_single_point(cube, 0.5, 358.5, "rotated", method="nonexistent")
 
 
 @pytest.mark.filterwarnings(
@@ -323,7 +357,9 @@ def test_boundary_warning(regrid_source_cube):
     with pytest.warns(
         regrid.BoundaryWarning, match="Selected point is within 8 gridlengths"
     ) as warning:
-        regrid.regrid_to_single_point(regrid_source_cube, -0.9, 393.0, "Nearest")
+        regrid.regrid_to_single_point(
+            regrid_source_cube, -0.9, 393.0, "rotated", "Nearest"
+        )
 
     assert len(warning) == 1
     assert issubclass(warning[-1].category, regrid.BoundaryWarning)