From e47273e37066dd067fe885e6113107b89f64ec20 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Thu, 9 Jul 2026 10:38:52 +0100 Subject: [PATCH 01/29] Placed sketch solution for rotating co-ordinate systems --- lib/ants/cli/ancil_create_shapefile.py | 63 ++++++++++++++++++++++++-- 1 file changed, 59 insertions(+), 4 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index 7081626..b9a4301 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -13,7 +13,9 @@ import argparse import json +import cartopy.crs as ccrs import numpy as np +from ants.io.load import load_landsea_mask from osgeo import ogr from shapely.geometry import Polygon @@ -39,7 +41,49 @@ def _load_polygon_from_json(json_file): return polygon -def main(json_file, output): +def _load_data(target_lsm_path): + """ + Loads in the target lsm which defines the rotated pole coordinate system + which is used as the new reference coordinate system for the shapefile. + + Parameters + ---------- + target_lsm_path : str + Path to the lsm file. + + Returns + ------- + + """ + + target_lsm = load_landsea_mask(target_lsm_path) + + return target_lsm + + +def _rotate_domains(target_lsm, lon, lat): + """ + Testing could consist of rotating through obvious angles + One set where the say the new pole is just rotated to the equator + + Two approaches are available: + 1. create an iris cube with the lat lon pairs + 2. just create a reference crs + """ + source_crs = ccrs.PlateCarree() + rotated_points = ( + target_lsm.coord_system() + .as_cartopy_crs() + .transform_points( + source_crs, # source_crs + lon, + lat, + ) + ) + return rotated_points + + +def main(json_file, output, target_lsm_path): """ Loads in a provided json file that defines pairs of longitude, latitude points to create a polygon from. That polygon is then used to create a @@ -51,12 +95,21 @@ def main(json_file, output): Path to json file output : str Location to store generated shape file - + target_lsm_path : str + File path for a land sea mask that provides the new rotated pole + coordinates to which the data be mapped. """ # Load a json and make a polygon polygon = _load_polygon_from_json(json_file) + # TODO replace with lat lon pairs from polygon file + lon, lat = [], [] + lsm = _load_data(target_lsm_path) + rotated_points = _rotate_domains(lsm, lon, lat) + # TODO assign rotated points to polygon (this is a placeholder for commit) + rotated_points = rotated_points + 1 + # Now convert it to a shapefile with OGR driver = ogr.GetDriverByName("Esri Shapefile") datasource = driver.CreateDataSource(output) @@ -82,7 +135,9 @@ def main(json_file, output): def _get_parser(): - parser = argparse.ArgumentParser() + parser = argparse.ArgumentParser( + target_lsm=True, + ) parser.add_argument( "json_file", help="Path to json file defining polygon to generate." ) @@ -93,7 +148,7 @@ def _get_parser(): def cli_interface(): parser = _get_parser() args = parser.parse_args() - main(args.json_file, args.output) + main(args.json_file, args.output, args.target_lsm) if __name__ == "__main__": From ca534fb9aa735a609e8630b92d2d86ed3d0c5d2b Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Mon, 13 Jul 2026 11:35:54 +0100 Subject: [PATCH 02/29] Allow points to be transformed to a rotated pole. --- lib/ants/cli/ancil_create_shapefile.py | 185 +++++++++++++----- .../rose-app.conf | 6 + 2 files changed, 142 insertions(+), 49 deletions(-) create mode 100644 rose-stem/app/ancil_create_shapefile_with_target_lsm/rose-app.conf diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index b9a4301..d5d90d4 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -9,78 +9,163 @@ Creates and saves a shapefile from a list of pairs of longitude, latitude points defining a single polygon in a specified polygon file. + +Rotated pole domains can be specified using the land sea mask argument, +where the longitude, latitude pairs are rotated to the new pole +It is assumed that the shapefile is on a standard unrotated lat-lon grid. """ import argparse import json +import logging +import warnings +import ants import cartopy.crs as ccrs +import iris.coord_systems import numpy as np from ants.io.load import load_landsea_mask from osgeo import ogr from shapely.geometry import Polygon +_LOGGER = logging.getLogger(__name__) + -def _load_polygon_from_json(json_file): +def _validate_coord_system(target_lsm): """ - Load a json file containing a list of pairs of longitude, latitude points - to create a polygon from. + Check that target_lsm has a rotated pole co-ordinate system. Parameters ---------- - json_file : str - Path to json file + target_lsm : :class:`iris.cube.Cube` + The lsm cube specifiying the rotated pole co-ordinate system. - Returns - ------- - : :class:`~shapely.geometry.Polygon` + Raises + ------ + ValueError : + If target_lsm does not have a rotated pole co-ordinate system. """ - with open(json_file, "r") as polygon_json: - polygon = json.load(polygon_json) - polygon = np.array(polygon) - polygon = Polygon(polygon) - return polygon + + is_pole_coords = isinstance( + target_lsm.coord_system(), iris.coord_systems.RotatedGeogCS + ) + if not is_pole_coords: + raise ValueError( + f"""target_lsm.coord_system() {target_lsm.coord_system()} is not + an instance of {iris.coord_systems.RotatedGeogCS}. The landsea mask + should specify a valid rotated pole co-ordinate system.""" + ) -def _load_data(target_lsm_path): +def _check_pole_orientation(target_lsm): """ - Loads in the target lsm which defines the rotated pole coordinate system - which is used as the new reference coordinate system for the shapefile. + Check if target_lsm is already orientated with the true north pole. Parameters ---------- - target_lsm_path : str - Path to the lsm file. + target_lsm : :class:`iris.cube.Cube` + The lsm cube specifiying the rotated pole co-ordinate system. + + Raises + ------ + warning : + If target_lsm is already coincident with the true north pole. + """ + + # TODO : I don't know what to specify here + sphere = ccrs.Globe(semimajor_axis=6371000.0, semiminor_axis=6371000.0) + + ref_crs = ccrs.RotatedGeodetic( + pole_latitude=90.0, + pole_longitude=0, + globe=sphere, + ) + + # get the cube proj4 parameters + target_crs = target_lsm.coord_system().as_cartopy_crs() + target_proj4 = target_crs.proj4_params + + if ref_crs == target_crs: + warnings.warn( + f"""target_lsm has a geodetic co-ordinate system with pole located at + grid_latitude={target_proj4['o_lat_p']}, + grid_longitude={target_proj4['o_lon_p']}. + No rotation will be carried out.""" + ) + + +def _transform_coordinates(target_lsm, lons, lats): + """ + Transforms the longitude, latitude points in the source co-ordinate + system to the rotated pole co-ordinate system defined by target_lsm. + The source co-ordinate is assumed to be unrotated geodetic. + + Parameters + ---------- + target_lsm : :class:`iris.cube.Cube` + The lsm cube specifiying the rotated pole co-ordinate system. + lons : :class:`np.ndarry` + Longitude points to rotate of length m. + lats : :class:`np.ndarry` + Latitude points to rotate of length m. Returns ------- + rotated_points : :class:`np.ndarry` + A (m,2) numpy array of rotated longitude and latitude pairs. """ - target_lsm = load_landsea_mask(target_lsm_path) + sphere = ccrs.Globe(semimajor_axis=6371000.0, semiminor_axis=6371000.0) + + source_crs = ccrs.Geodetic(globe=sphere) + target_coord = target_lsm.coord_system() - return target_lsm + rotated_points = target_coord.as_cartopy_crs().transform_points( + source_crs, lons, lats # longitude # latitude + )[:, :2] + + _LOGGER.info( + "Input json file transformed to new pole rotated co-ordinate system at" + "grid_latitude=%s, grid_longitude=%s.", + target_coord.grid_north_pole_latitude, + target_coord.grid_north_pole_longitude, + ) + + return rotated_points -def _rotate_domains(target_lsm, lon, lat): +def _load_polygon_from_json(json_file, target_lsm_path): """ - Testing could consist of rotating through obvious angles - One set where the say the new pole is just rotated to the equator + Load a json file containing a list of pairs of longitude, latitude points + to create a polygon from. - Two approaches are available: - 1. create an iris cube with the lat lon pairs - 2. just create a reference crs + Parameters + ---------- + json_file : str + Path to json file + target_lsm_path : str + File path for a land sea mask that provides the new rotated pole + coordinates to which the longitude, latitude pairs will be mapped. + + Returns + ------- + : :class:`~shapely.geometry.Polygon` """ - source_crs = ccrs.PlateCarree() - rotated_points = ( - target_lsm.coord_system() - .as_cartopy_crs() - .transform_points( - source_crs, # source_crs - lon, - lat, - ) - ) - return rotated_points + with open(json_file, "r") as polygon_json: + points = json.load(polygon_json) + points = np.array(points) + + if target_lsm_path is not None: + target_lsm = load_landsea_mask(target_lsm_path) + _validate_coord_system(target_lsm) + _check_pole_orientation(target_lsm) + + lon, lat = points[:, 0], points[:, -1] + points = _transform_coordinates(target_lsm, lon, lat) + + polygon = Polygon(points) + + return polygon def main(json_file, output, target_lsm_path): @@ -89,6 +174,10 @@ def main(json_file, output, target_lsm_path): points to create a polygon from. That polygon is then used to create a shape file that is saved to the specified output location. + If a target_lsm_path is provided, the points will first be transformed + from a geodetic co-ordinate system to a rotated pole co-ordinate system + specified by the lsm file. + Parameters ---------- json_file : str @@ -97,18 +186,11 @@ def main(json_file, output, target_lsm_path): Location to store generated shape file target_lsm_path : str File path for a land sea mask that provides the new rotated pole - coordinates to which the data be mapped. + coordinates to which the longitude, latitude pairs will be transformed. """ # Load a json and make a polygon - polygon = _load_polygon_from_json(json_file) - - # TODO replace with lat lon pairs from polygon file - lon, lat = [], [] - lsm = _load_data(target_lsm_path) - rotated_points = _rotate_domains(lsm, lon, lat) - # TODO assign rotated points to polygon (this is a placeholder for commit) - rotated_points = rotated_points + 1 + polygon = _load_polygon_from_json(json_file, target_lsm_path) # Now convert it to a shapefile with OGR driver = ogr.GetDriverByName("Esri Shapefile") @@ -135,12 +217,17 @@ def main(json_file, output, target_lsm_path): def _get_parser(): - parser = argparse.ArgumentParser( - target_lsm=True, - ) + parser = argparse.ArgumentParser() parser.add_argument( "json_file", help="Path to json file defining polygon to generate." ) + parser.add_argument( + "--target-lsm", + type=ants.config.filepath_readable, + required=False, + help="Path to the land sea mask containing the rotated pole" + " co-ordinate system.", + ) parser.add_argument("output", help="File to save shape file to.") return parser diff --git a/rose-stem/app/ancil_create_shapefile_with_target_lsm/rose-app.conf b/rose-stem/app/ancil_create_shapefile_with_target_lsm/rose-app.conf new file mode 100644 index 0000000..4f37bbf --- /dev/null +++ b/rose-stem/app/ancil_create_shapefile_with_target_lsm/rose-app.conf @@ -0,0 +1,6 @@ +[command] +default=ants-launch ancil_create_shapefile.py ${SHAPEFILE_JSON} ${OUTPUT} + +[env] +OUTPUT=${ROSE_DATA}/ite.shp +SHAPEFILE_JSON=${TEST_SOURCES_DIR}/ancil_create_shapefile/ite.json From 2d8879eef5f1f153c519aeefb172107f28d6b51f Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Tue, 14 Jul 2026 11:10:14 +0100 Subject: [PATCH 03/29] Testing with optional argument in rose-stem. --- lib/ants/cli/ancil_create_shapefile.py | 97 ++++++++----------- .../app/ancil_create_shapefile/rose-app.conf | 3 +- .../rose-app.conf | 6 -- 3 files changed, 45 insertions(+), 61 deletions(-) delete mode 100644 rose-stem/app/ancil_create_shapefile_with_target_lsm/rose-app.conf diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index d5d90d4..b45dfa7 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -12,7 +12,7 @@ Rotated pole domains can be specified using the land sea mask argument, where the longitude, latitude pairs are rotated to the new pole -It is assumed that the shapefile is on a standard unrotated lat-lon grid. +It is assumed that the shapefile is on a standard unrotated lon-lat grid. """ import argparse import json @@ -32,17 +32,17 @@ def _validate_coord_system(target_lsm): """ - Check that target_lsm has a rotated pole co-ordinate system. + Check that target_lsm has a rotated pole coordinate system. Parameters ---------- target_lsm : :class:`iris.cube.Cube` - The lsm cube specifiying the rotated pole co-ordinate system. + The lsm cube specifiying the rotated pole coordinate system. Raises ------ ValueError : - If target_lsm does not have a rotated pole co-ordinate system. + If target_lsm does not have a rotated pole coordinate system. """ is_pole_coords = isinstance( @@ -52,28 +52,43 @@ def _validate_coord_system(target_lsm): raise ValueError( f"""target_lsm.coord_system() {target_lsm.coord_system()} is not an instance of {iris.coord_systems.RotatedGeogCS}. The landsea mask - should specify a valid rotated pole co-ordinate system.""" + should specify a valid rotated pole coordinate system.""" ) -def _check_pole_orientation(target_lsm): +def _transform_coordinates(target_lsm, points): """ - Check if target_lsm is already orientated with the true north pole. + Transforms the longitude, latitude points in the source coordinate + system to the rotated pole coordinate system defined by target_lsm. + + he source coordinate system is assumed to be unrotated geodetic + defined on a sphere. Does nothing if the target coordinate system + is unrotated. Parameters ---------- target_lsm : :class:`iris.cube.Cube` - The lsm cube specifiying the rotated pole co-ordinate system. + The lsm cube specifiying the rotated pole coordinate system. + points : :class:`np.ndarry` + An (m,2) sized numpy array of m longitude and m latitude points. + + Returns + ------- + rotated_points : :class:`np.ndarry` + An (m,2) sized numpy array of rotated longitude and latitude points. Raises ------ - warning : - If target_lsm is already coincident with the true north pole. + Warns : UserWarning + If target_lsm is coincident with the true north pole and is + defined on a sphere. """ - # TODO : I don't know what to specify here sphere = ccrs.Globe(semimajor_axis=6371000.0, semiminor_axis=6371000.0) + source_crs = ccrs.Geodetic(globe=sphere) + target_coord = target_lsm.coord_system() + ref_crs = ccrs.RotatedGeodetic( pole_latitude=90.0, pole_longitude=0, @@ -84,51 +99,27 @@ def _check_pole_orientation(target_lsm): target_crs = target_lsm.coord_system().as_cartopy_crs() target_proj4 = target_crs.proj4_params + # Checks that the poles and also underlying ellipses are the same if ref_crs == target_crs: warnings.warn( - f"""target_lsm has a geodetic co-ordinate system with pole located at + f"""target_lsm has a geodetic coordinate system with pole located at grid_latitude={target_proj4['o_lat_p']}, grid_longitude={target_proj4['o_lon_p']}. - No rotation will be carried out.""" + No transformation will be carried out.""" ) - - -def _transform_coordinates(target_lsm, lons, lats): - """ - Transforms the longitude, latitude points in the source co-ordinate - system to the rotated pole co-ordinate system defined by target_lsm. - The source co-ordinate is assumed to be unrotated geodetic. - - Parameters - ---------- - target_lsm : :class:`iris.cube.Cube` - The lsm cube specifiying the rotated pole co-ordinate system. - lons : :class:`np.ndarry` - Longitude points to rotate of length m. - lats : :class:`np.ndarry` - Latitude points to rotate of length m. - - Returns - ------- - rotated_points : :class:`np.ndarry` - A (m,2) numpy array of rotated longitude and latitude pairs. - - """ - - sphere = ccrs.Globe(semimajor_axis=6371000.0, semiminor_axis=6371000.0) - - source_crs = ccrs.Geodetic(globe=sphere) - target_coord = target_lsm.coord_system() + rotated_points = np.copy(points) + return rotated_points rotated_points = target_coord.as_cartopy_crs().transform_points( - source_crs, lons, lats # longitude # latitude + source_crs, points[:, 0], points[:, 1] )[:, :2] _LOGGER.info( - "Input json file transformed to new pole rotated co-ordinate system at" - "grid_latitude=%s, grid_longitude=%s.", - target_coord.grid_north_pole_latitude, + "Input json file transformed to new pole rotated coordinate system at" + "pole longitude=%s, pole latitude=%s, central rotated longitude=%s.", target_coord.grid_north_pole_longitude, + target_coord.grid_north_pole_latitude, + target_coord.north_pole_grid_longitude, ) return rotated_points @@ -145,7 +136,7 @@ def _load_polygon_from_json(json_file, target_lsm_path): Path to json file target_lsm_path : str File path for a land sea mask that provides the new rotated pole - coordinates to which the longitude, latitude pairs will be mapped. + coordinates to which the longitude, latitude pairs will be transformed. Returns ------- @@ -158,10 +149,7 @@ def _load_polygon_from_json(json_file, target_lsm_path): if target_lsm_path is not None: target_lsm = load_landsea_mask(target_lsm_path) _validate_coord_system(target_lsm) - _check_pole_orientation(target_lsm) - - lon, lat = points[:, 0], points[:, -1] - points = _transform_coordinates(target_lsm, lon, lat) + points = _transform_coordinates(target_lsm, points) polygon = Polygon(points) @@ -174,9 +162,10 @@ def main(json_file, output, target_lsm_path): points to create a polygon from. That polygon is then used to create a shape file that is saved to the specified output location. - If a target_lsm_path is provided, the points will first be transformed - from a geodetic co-ordinate system to a rotated pole co-ordinate system - specified by the lsm file. + If target_lsm_path is provided, the points are first transformed from an + unrotated geodetic coordinate system to a rotated pole coordinate system + specified by the lsm. It is assumed that the points specified in the json + file are on an unrotated geodetic grid. Parameters ---------- @@ -226,7 +215,7 @@ def _get_parser(): type=ants.config.filepath_readable, required=False, help="Path to the land sea mask containing the rotated pole" - " co-ordinate system.", + " coordinate system.", ) parser.add_argument("output", help="File to save shape file to.") return parser diff --git a/rose-stem/app/ancil_create_shapefile/rose-app.conf b/rose-stem/app/ancil_create_shapefile/rose-app.conf index 4f37bbf..05f94bd 100644 --- a/rose-stem/app/ancil_create_shapefile/rose-app.conf +++ b/rose-stem/app/ancil_create_shapefile/rose-app.conf @@ -1,6 +1,7 @@ [command] -default=ants-launch ancil_create_shapefile.py ${SHAPEFILE_JSON} ${OUTPUT} +default=ants-launch ancil_create_shapefile.py ${SHAPEFILE_JSON} ${OUTPUT} --target-lsm ${LSM} [env] OUTPUT=${ROSE_DATA}/ite.shp SHAPEFILE_JSON=${TEST_SOURCES_DIR}/ancil_create_shapefile/ite.json +LSM=/home/users/alasdair.roy/Documents/cartopy_tutorial/target_lsm.nc diff --git a/rose-stem/app/ancil_create_shapefile_with_target_lsm/rose-app.conf b/rose-stem/app/ancil_create_shapefile_with_target_lsm/rose-app.conf deleted file mode 100644 index 4f37bbf..0000000 --- a/rose-stem/app/ancil_create_shapefile_with_target_lsm/rose-app.conf +++ /dev/null @@ -1,6 +0,0 @@ -[command] -default=ants-launch ancil_create_shapefile.py ${SHAPEFILE_JSON} ${OUTPUT} - -[env] -OUTPUT=${ROSE_DATA}/ite.shp -SHAPEFILE_JSON=${TEST_SOURCES_DIR}/ancil_create_shapefile/ite.json From 36a38190d091b90eb1ba86702971129571cf9869 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Wed, 15 Jul 2026 11:49:02 +0100 Subject: [PATCH 04/29] Added source cube as argument, added some tests for shapefile creation. --- lib/ants/cli/ancil_create_shapefile.py | 102 +++++++----- lib/ants/tests/cli/__init__.py | 4 + .../tests/cli/test_ancil_create_shapefile.py | 151 ++++++++++++++++++ 3 files changed, 220 insertions(+), 37 deletions(-) create mode 100644 lib/ants/tests/cli/__init__.py create mode 100644 lib/ants/tests/cli/test_ancil_create_shapefile.py diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index b45dfa7..2a6852c 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -13,6 +13,8 @@ Rotated pole domains can be specified using the land sea mask argument, where the longitude, latitude pairs are rotated to the new pole It is assumed that the shapefile is on a standard unrotated lon-lat grid. + +# TODO I need to decide if the central rotated longitude has to be set to 180 """ import argparse import json @@ -23,7 +25,7 @@ import cartopy.crs as ccrs import iris.coord_systems import numpy as np -from ants.io.load import load_landsea_mask +from ants.io.load import load_cube, load_landsea_mask from osgeo import ogr from shapely.geometry import Polygon @@ -50,18 +52,39 @@ def _validate_coord_system(target_lsm): ) if not is_pole_coords: raise ValueError( - f"""target_lsm.coord_system() {target_lsm.coord_system()} is not - an instance of {iris.coord_systems.RotatedGeogCS}. The landsea mask - should specify a valid rotated pole coordinate system.""" + f"target_lsm.coord_system() {target_lsm.coord_system()} is not" + f" an instance of {iris.coord_systems.RotatedGeogCS}." + f" The landsea mask should specify a valid rotated pole coordinate" + f" system." ) -def _transform_coordinates(target_lsm, points): +def _validate_orientation(target_lsm): + + ref_crs = ccrs.RotatedGeodetic(pole_latitude=90.0, pole_longitude=0.0) + + # get the cube proj4 parameters + target_crs = target_lsm.coord_system().as_cartopy_crs() + target_proj4 = target_crs.proj4_params + + # Checks that the poles and also underlying ellipses are the same + check_crs = ref_crs == target_crs + if check_crs: + warnings.warn( + f"""target_lsm has a geodetic coordinate system with pole located at + grid_latitude={target_proj4['o_lat_p']}, + grid_longitude={target_proj4['o_lon_p']}. + No transformation will be carried out.""" + ) + return check_crs + + +def _transform_coordinates(target_lsm, source_cube, points): """ Transforms the longitude, latitude points in the source coordinate system to the rotated pole coordinate system defined by target_lsm. - he source coordinate system is assumed to be unrotated geodetic + The source coordinate system is assumed to be unrotated geodetic defined on a sphere. Does nothing if the target coordinate system is unrotated. @@ -69,6 +92,9 @@ def _transform_coordinates(target_lsm, points): ---------- target_lsm : :class:`iris.cube.Cube` The lsm cube specifiying the rotated pole coordinate system. + source_cube : :class:`iris.cube.Cube` + An iris cube specifying the co-ordinate system of the + input json file. points : :class:`np.ndarry` An (m,2) sized numpy array of m longitude and m latitude points. @@ -84,36 +110,22 @@ def _transform_coordinates(target_lsm, points): defined on a sphere. """ - sphere = ccrs.Globe(semimajor_axis=6371000.0, semiminor_axis=6371000.0) - - source_crs = ccrs.Geodetic(globe=sphere) - target_coord = target_lsm.coord_system() + # sphere = ccrs.Globe(semimajor_axis=6371000.0, semiminor_axis=6371000.0) - ref_crs = ccrs.RotatedGeodetic( - pole_latitude=90.0, - pole_longitude=0, - globe=sphere, - ) + # source_crs = ccrs.Geodetic(globe=sphere) - # get the cube proj4 parameters - target_crs = target_lsm.coord_system().as_cartopy_crs() - target_proj4 = target_crs.proj4_params - - # Checks that the poles and also underlying ellipses are the same - if ref_crs == target_crs: - warnings.warn( - f"""target_lsm has a geodetic coordinate system with pole located at - grid_latitude={target_proj4['o_lat_p']}, - grid_longitude={target_proj4['o_lon_p']}. - No transformation will be carried out.""" - ) - rotated_points = np.copy(points) - return rotated_points + source_crs = source_cube.coord_system().as_cartopy_crs() + target_coord = target_lsm.coord_system() + target_crs = target_coord.as_cartopy_crs() - rotated_points = target_coord.as_cartopy_crs().transform_points( + rotated_points = target_crs.transform_points( source_crs, points[:, 0], points[:, 1] )[:, :2] + # Enforce that the points are all measured in the positive direction + negative_rotated_lons = np.where(rotated_points[:, 0] < 0) + rotated_points[negative_rotated_lons, 0] += 360.0 + _LOGGER.info( "Input json file transformed to new pole rotated coordinate system at" "pole longitude=%s, pole latitude=%s, central rotated longitude=%s.", @@ -125,7 +137,7 @@ def _transform_coordinates(target_lsm, points): return rotated_points -def _load_polygon_from_json(json_file, target_lsm_path): +def _load_polygon_from_json(json_file, target_lsm_path, source_cube_path): """ Load a json file containing a list of pairs of longitude, latitude points to create a polygon from. @@ -137,6 +149,9 @@ def _load_polygon_from_json(json_file, target_lsm_path): target_lsm_path : str File path for a land sea mask that provides the new rotated pole coordinates to which the longitude, latitude pairs will be transformed. + source_cube_path : str + File path to an iris cube specifying the co-ordinate system of the + input json file. Returns ------- @@ -146,17 +161,19 @@ def _load_polygon_from_json(json_file, target_lsm_path): points = json.load(polygon_json) points = np.array(points) - if target_lsm_path is not None: + if target_lsm_path is not None and source_cube_path is not None: target_lsm = load_landsea_mask(target_lsm_path) - _validate_coord_system(target_lsm) - points = _transform_coordinates(target_lsm, points) + source_cube = load_cube(source_cube_path) + + if not _validate_coord_system(target_lsm): + points = _transform_coordinates(target_lsm, source_cube, points) polygon = Polygon(points) return polygon -def main(json_file, output, target_lsm_path): +def main(json_file, output, target_lsm_path, source_cube_path): """ Loads in a provided json file that defines pairs of longitude, latitude points to create a polygon from. That polygon is then used to create a @@ -176,10 +193,13 @@ def main(json_file, output, target_lsm_path): target_lsm_path : str File path for a land sea mask that provides the new rotated pole coordinates to which the longitude, latitude pairs will be transformed. + source_cube_path : str + File path to an iris cube specifying the co-ordinate system of the + input json file. """ # Load a json and make a polygon - polygon = _load_polygon_from_json(json_file, target_lsm_path) + polygon = _load_polygon_from_json(json_file, target_lsm_path, source_cube_path) # Now convert it to a shapefile with OGR driver = ogr.GetDriverByName("Esri Shapefile") @@ -217,6 +237,12 @@ def _get_parser(): help="Path to the land sea mask containing the rotated pole" " coordinate system.", ) + parser.add_argument( + "--source-cube", + type=ants.config.filepath_readable, + required=False, + help="Path to the source co-ordinate system", + ) parser.add_argument("output", help="File to save shape file to.") return parser @@ -224,7 +250,9 @@ def _get_parser(): def cli_interface(): parser = _get_parser() args = parser.parse_args() - main(args.json_file, args.output, args.target_lsm) + if args.search_method is not None and args.source_cube is None: + parser.error("--source-cube is required when --target-lsm is specified") + main(args.json_file, args.output, args.target_lsm, args.source_cube) if __name__ == "__main__": diff --git a/lib/ants/tests/cli/__init__.py b/lib/ants/tests/cli/__init__.py new file mode 100644 index 0000000..89193a5 --- /dev/null +++ b/lib/ants/tests/cli/__init__.py @@ -0,0 +1,4 @@ +# (C) Crown Copyright, Met Office. All rights reserved. +# +# This file is part of ANTS and is released under the BSD 3-Clause license. +# See LICENSE.txt in the root of the repository for full licensing details. diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py new file mode 100644 index 0000000..b34a432 --- /dev/null +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -0,0 +1,151 @@ +# (C) Crown Copyright, Met Office. All rights reserved. +# +# This file is part of ANTS and is released under the BSD 3-Clause license. +# See LICENSE.txt in the root of the repository for full licensing details. + +import re +import unittest.mock as mock + +import ants.tests +import iris +import numpy as np +from ants.cli.ancil_create_shapefile import ( + _transform_coordinates, + _validate_coord_system, +) +from ants.tests.stock import geodetic + + +def _create_cube(north_pole_lon, north_pole_lat, ellipsoid): + # stock geodetic automatically returns an identity rotated geodetic + # back to geodetic + crs = iris.coord_systems.RotatedGeogCS( + grid_north_pole_latitude=north_pole_lat, + grid_north_pole_longitude=north_pole_lon, + ellipsoid=ellipsoid, + ) + lons, lats = [0.0], [0.0] + lon_coord = iris.coords.DimCoord( + lons, + standard_name="grid_longitude", + units="degrees", + coord_system=crs, + ) + + lat_coord = iris.coords.DimCoord( + lats, + standard_name="grid_latitude", + units="degrees", + coord_system=crs, + ) + + data = np.zeros((len(lats), len(lons))) + + cube = iris.cube.Cube(data, dim_coords_and_dims=[(lat_coord, 0), (lon_coord, 1)]) + return cube + + +class Common(object): + def setUp(self): + # unrotated source cube + self.source_cube = geodetic((2, 2)) + # spherical input coordinate systems are simple to test + self.sphere_crs = iris.coord_systems.GeogCS(6371229.0, 6371229.0) + + self.sphere_source = geodetic((2, 2), crs=self.sphere_crs) + + self.sphere_identity_target_lsm = _create_cube(0.0, 90.0, self.sphere_crs) + self.sphere_equator_target_lsm = geodetic( + (2, 2), north_pole_lat=0.0, north_pole_lon=0.0, crs=self.sphere_crs + ) + + +class Test__validate_coord_system(Common, ants.tests.TestCase): + + def test_unrotated_target_lsm(self): + target_lsm = geodetic((2, 2)) + error_msg = re.escape( + f"target_lsm.coord_system() {target_lsm.coord_system()} is not" + f" an instance of {iris.coord_systems.RotatedGeogCS}." + f" The landsea mask should specify a valid rotated pole coordinate" + f" system." + ) + + with self.assertRaisesRegex(ValueError, error_msg): + _validate_coord_system(target_lsm) + + +# add test that checks results are made positive with a mock +class Test__transform_coordinates(Common, ants.tests.TestCase): + def setUp(self): + super().setUp() + self.vary_latitudes = np.array( + [[0.0, 0.0, 0.0, 0.0, 0.0], [90.0, 45.0, 0.0, -45.0, -90.0]] + ).T + self.vary_longitudes = np.array( + [ + [-135.0, -90.0, -45.0, 0.0, 45.0, 90.0, 135.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + ] + ).T + + def test_identity_rotation_sphere(self): + + true_lats = np.array([90, 45, 0, -45, -90]) + # The rotated geodesic system flips the longitude + # Rotated geodesic has some odd properties, like the central + # longitude being set to 180.0 + true_lons = np.array([np.nan, 180.0, 180.0, 180.0, 180.0]) + + expected_rotation = np.array([true_lons, true_lats]).T + # we cannot check the longitude at the pole (any longitude is valid) + check_mask = ~np.isnan(expected_rotation) + + rotated_coords = _transform_coordinates( + self.sphere_identity_target_lsm, self.sphere_source, self.vary_latitudes + ) + + assert ants.utils.ndarray.allclose( + expected_rotation[check_mask], rotated_coords[check_mask] + ) + + def test_rotation_to_equator_sphere(self): + true_lats = np.array([-45.0, 0.0, 45.0, 90.0, 45.0, 0.0, -45.0]) + # Again the results for lon depend on where the central lon is + true_lons = np.array([90.0, 90.0, 90.0, np.nan, 270.0, 270.0, 270.0]) + expected_rotation = np.array([true_lons, true_lats]).T + check_mask = ~np.isnan(expected_rotation) + rotated_coords = _transform_coordinates( + self.sphere_equator_target_lsm, self.sphere_source, self.vary_longitudes + ) + assert ants.utils.ndarray.allclose( + expected_rotation[check_mask], rotated_coords[check_mask] + ) + + def test_negative_longitudes_converted(self): + + true_coords = np.array([[90.0, 160.0, 326.7], [0.0, 0.0, 0.0]]).T + # presumably I can't get a -400.0 for example + neg_lons = np.array([[-270.0, -200.0, -33.3], [0.0, 0.0, 0.0]]).T + + target_lsm, source_cube = mock.Mock(), mock.Mock() + source_crs, target_crs = mock.Mock(), mock.Mock() + target_coord = mock.Mock() + + source_cube.coord_system.return_value.as_cartopy_crs.return_value = source_crs + target_lsm.coord_system.return_value = target_coord + + target_coord.as_cartopy_crs.return_value = target_crs + + target_crs.transform_points.return_value = neg_lons + + result = _transform_coordinates(target_lsm, source_cube, neg_lons) + assert ants.utils.ndarray.allclose(true_coords, result) + + +class test_main: + pass + + +class test_cli: + pass From a6a60eb4481684ccc71dc0ae9c5489ec2e2574b7 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Thu, 16 Jul 2026 16:32:57 +0100 Subject: [PATCH 05/29] Added more tests --- lib/ants/cli/ancil_create_shapefile.py | 62 ++++-- .../tests/cli/test_ancil_create_shapefile.py | 200 ++++++++++++++++-- 2 files changed, 225 insertions(+), 37 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index 2a6852c..5818d5a 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -22,7 +22,6 @@ import warnings import ants -import cartopy.crs as ccrs import iris.coord_systems import numpy as np from ants.io.load import load_cube, load_landsea_mask @@ -60,21 +59,31 @@ def _validate_coord_system(target_lsm): def _validate_orientation(target_lsm): + """ + Check if target_lsm is a non-rotated pole. - ref_crs = ccrs.RotatedGeodetic(pole_latitude=90.0, pole_longitude=0.0) + Parameters + ---------- + target_lsm : :class:`iris.cube.Cube` + The lsm cube specifiying the rotated pole coordinate system. - # get the cube proj4 parameters - target_crs = target_lsm.coord_system().as_cartopy_crs() - target_proj4 = target_crs.proj4_params + Warns + ------ + UserWarning : + If target_lsm is already at latitude=90.0, longitude=0.0. + """ + + target_crs = target_lsm.coord_system() + + grid_lon = target_crs.grid_north_pole_longitude + grid_lat = target_crs.grid_north_pole_latitude - # Checks that the poles and also underlying ellipses are the same - check_crs = ref_crs == target_crs + check_crs = ants.utils.ndarray.allclose([grid_lon, grid_lat], [0.0, 90.0]) if check_crs: warnings.warn( - f"""target_lsm has a geodetic coordinate system with pole located at - grid_latitude={target_proj4['o_lat_p']}, - grid_longitude={target_proj4['o_lon_p']}. - No transformation will be carried out.""" + "target_lsm has a geodetic coordinate system with pole located" + f"at grid_longitude={grid_lon}, grid_latitude={grid_lat}." + "No transformation will be carried out." ) return check_crs @@ -102,12 +111,6 @@ def _transform_coordinates(target_lsm, source_cube, points): ------- rotated_points : :class:`np.ndarry` An (m,2) sized numpy array of rotated longitude and latitude points. - - Raises - ------ - Warns : UserWarning - If target_lsm is coincident with the true north pole and is - defined on a sphere. """ # sphere = ccrs.Globe(semimajor_axis=6371000.0, semiminor_axis=6371000.0) @@ -161,10 +164,14 @@ def _load_polygon_from_json(json_file, target_lsm_path, source_cube_path): points = json.load(polygon_json) points = np.array(points) - if target_lsm_path is not None and source_cube_path is not None: - target_lsm = load_landsea_mask(target_lsm_path) + if target_lsm_path is not None and source_cube_path is None: + crs = iris.coord_systems.GeogCS(6371229.0) + source_cube = ants.utils.cube.CubeBuilder(crs, (2, 2))._cube + elif target_lsm_path is not None and source_cube_path is not None: source_cube = load_cube(source_cube_path) + if target_lsm_path is not None: + target_lsm = load_landsea_mask(target_lsm_path) if not _validate_coord_system(target_lsm): points = _transform_coordinates(target_lsm, source_cube, points) @@ -230,6 +237,7 @@ def _get_parser(): parser.add_argument( "json_file", help="Path to json file defining polygon to generate." ) + parser.add_argument("output", help="File to save shape file to.") parser.add_argument( "--target-lsm", type=ants.config.filepath_readable, @@ -241,17 +249,25 @@ def _get_parser(): "--source-cube", type=ants.config.filepath_readable, required=False, - help="Path to the source co-ordinate system", + help="Path to an iris cube which specifies the co-ordinate" + "system of the json file", ) - parser.add_argument("output", help="File to save shape file to.") return parser +def _validate_args(args): + + if args.source_cube is not None and args.target_lsm is None: + raise ValueError( + "If --source-cube is passed then --target-lsm must" "also be given." + ) + + def cli_interface(): parser = _get_parser() args = parser.parse_args() - if args.search_method is not None and args.source_cube is None: - parser.error("--source-cube is required when --target-lsm is specified") + + _validate_args(args) main(args.json_file, args.output, args.target_lsm, args.source_cube) diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index b34a432..71324ac 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -3,15 +3,20 @@ # This file is part of ANTS and is released under the BSD 3-Clause license. # See LICENSE.txt in the root of the repository for full licensing details. +import argparse import re -import unittest.mock as mock +from unittest import mock +from unittest.mock import Mock import ants.tests import iris import numpy as np from ants.cli.ancil_create_shapefile import ( + _load_polygon_from_json, _transform_coordinates, + _validate_args, _validate_coord_system, + _validate_orientation, ) from ants.tests.stock import geodetic @@ -50,14 +55,15 @@ def setUp(self): # unrotated source cube self.source_cube = geodetic((2, 2)) # spherical input coordinate systems are simple to test - self.sphere_crs = iris.coord_systems.GeogCS(6371229.0, 6371229.0) - + self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) self.sphere_source = geodetic((2, 2), crs=self.sphere_crs) - self.sphere_identity_target_lsm = _create_cube(0.0, 90.0, self.sphere_crs) self.sphere_equator_target_lsm = geodetic( (2, 2), north_pole_lat=0.0, north_pole_lon=0.0, crs=self.sphere_crs ) + self.sphere_rotate_lon = geodetic( + (2, 2), north_pole_lat=90.0, north_pole_lon=90.0, crs=self.sphere_crs + ) class Test__validate_coord_system(Common, ants.tests.TestCase): @@ -75,7 +81,23 @@ def test_unrotated_target_lsm(self): _validate_coord_system(target_lsm) -# add test that checks results are made positive with a mock +class Test__validate_orientation(Common, ants.tests.TestCase): + def setUp(self): + super().setUp() + self.warning_msg = ( + "target_lsm has a geodetic coordinate system with pole located" + "at grid_longitude=0.0, grid_latitude=90.0." + "No transformation will be carried out." + ) + + def test_warning_raised(self): + with self.assertWarnsRegex(UserWarning, self.warning_msg): + _validate_orientation(self.sphere_identity_target_lsm) + + def test_no_warning_raised(self): + self.assertFalse(_validate_orientation(self.sphere_rotate_lon)) + + class Test__transform_coordinates(Common, ants.tests.TestCase): def setUp(self): super().setUp() @@ -109,6 +131,24 @@ def test_identity_rotation_sphere(self): expected_rotation[check_mask], rotated_coords[check_mask] ) + def test_longitudinal_rotation(self): + # these tests can be counter-intuitive as the prime-meridian + # is by default at 180.0 so a rotation of 90.0 causes the prime meridian + # to be at 270.0 (it is always 180 further round from the longitude given) + points = np.array([[0.0, 45.0, 90.0, 135.0], [0.0, 0.0, 0.0, 0.0]]).T + + true_lats = np.array([0.0, 0.0, np.nan, 0.0]) + true_lons = np.array([90.0, 135.0, 180.0, 225.0]) + expected_rotation = np.array([true_lons, true_lats]).T + check_mask = ~np.isnan(expected_rotation) + + rotated_coords = _transform_coordinates( + self.sphere_rotate_lon, self.sphere_source, points + ) + assert ants.utils.ndarray.allclose( + expected_rotation[check_mask], rotated_coords[check_mask] + ) + def test_rotation_to_equator_sphere(self): true_lats = np.array([-45.0, 0.0, 45.0, 90.0, 45.0, 0.0, -45.0]) # Again the results for lon depend on where the central lon is @@ -128,9 +168,9 @@ def test_negative_longitudes_converted(self): # presumably I can't get a -400.0 for example neg_lons = np.array([[-270.0, -200.0, -33.3], [0.0, 0.0, 0.0]]).T - target_lsm, source_cube = mock.Mock(), mock.Mock() - source_crs, target_crs = mock.Mock(), mock.Mock() - target_coord = mock.Mock() + target_lsm, source_cube = Mock(), Mock() + source_crs, target_crs = Mock(), Mock() + target_coord = Mock() source_cube.coord_system.return_value.as_cartopy_crs.return_value = source_crs target_lsm.coord_system.return_value = target_coord @@ -143,9 +183,141 @@ def test_negative_longitudes_converted(self): assert ants.utils.ndarray.allclose(true_coords, result) -class test_main: - pass - - -class test_cli: - pass +class Test__load_polygon_from_json(Common, ants.tests.TestCase): + def setUp(self): + super().setUp() + self.json_values = [[1, 2], [3, 4], [5, 6], [7, 8]] + + def patch_loader(func): + @mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") + @mock.patch("ants.cli.ancil_create_shapefile.load_cube") + @mock.patch("ants.cli.ancil_create_shapefile._validate_coord_system") + @mock.patch("ants.cli.ancil_create_shapefile._transform_coordinates") + @mock.patch("builtins.open", new_callable=mock.mock_open) + @mock.patch("ants.cli.ancil_create_shapefile.json.load") + def wrapper( + self, + mock_json, + mock_open, + mock_transform, + mock_validate, + mock_load, + mock_load_lsm, + *args, + **kwargs, + ): + return func( + self, + mock_json, + mock_open, + mock_transform, + mock_validate, + mock_load, + mock_load_lsm, + *args, + **kwargs, + ) + + return wrapper + + @patch_loader + def test_load_call( + self, + mock_json, + mock_open, + mock_transform, + mock_validate, + mock_load, + mock_load_lsm, + ): + + mock_json.return_value = self.json_values + _ = _load_polygon_from_json("json/path", None, None) + + mock_json.assert_called_once() + mock_load_lsm.assert_not_called() + mock_load.assert_not_called() + mock_validate.assert_not_called() + mock_transform.assert_not_called() + + @patch_loader + def test_transform_not_called( + self, + mock_json, + mock_open, + mock_transform, + mock_validate, + mock_load, + mock_load_lsm, + ): + + mock_json.return_value = self.json_values + mock_validate.return_value = True + mock_load_lsm.return_value = Mock() + mock_load.return_value = Mock() + + _ = _load_polygon_from_json("json/path", "lsm/path", "source_cube/path") + + mock_load_lsm.assert_called_once() + mock_load.assert_called_once() + mock_transform.assert_not_called() + + @patch_loader + def test_transform_called( + self, + mock_json, + mock_open, + mock_transform, + mock_validate, + mock_load, + mock_load_lsm, + ): + + mock_json.return_value = self.json_values + mock_validate.return_value = False + mock_load_lsm.return_value = Mock() + mock_load.return_value = Mock() + mock_transform.return_value = self.json_values + + _ = _load_polygon_from_json("json/path", "lsm/path", "source_cube/path") + + mock_load_lsm.assert_called_once() + mock_load.assert_called_once() + mock_transform.assert_called_once() + + """ + I need to figure out how to get this patch to work + @patch_loader + def test_load_call(self, mock_json, mock_open, mock_transform, + mock_validate, mock_load, mock_load_lsm): + + ants_cube = ants.utils.cube.CubeBuilder(self.sphere_crs,(2,2))._cube + mock_json.return_value = self.json_values + mock_validate.return_value = False + mock_load_lsm.return_value = Mock() + mock_transform.return_value = self.json_values + + _ = _load_polygon_from_json("json/path", "lsm/path", None) + + mock_json.assert_called_once() + mock_load_lsm.assert_called_once() + mock_load.assert_not_called() + mock_validate.assert_called_once() + mock_transform.assert_called_once_with(mock_load_lsm.return_value, + ants_cube, + np.array(self.json_values)) + """ + + +class Test__validate_args(ants.tests.TestCase): + + def test_error_raised(self): + args = argparse.Namespace( + json_file="json/path", + output="output/path", + target_lsm=None, + source_cube="source/path", + ) + error_msg = "If --source-cube is passed then --target-lsm must" "also be given." + with self.assertRaisesRegex(ValueError, error_msg): + _validate_args(args) From 64af015384edeb22fa79b1a24eb17ccaa164f986 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Fri, 17 Jul 2026 11:18:44 +0100 Subject: [PATCH 06/29] Changed source cube to optional and cleaned tests --- lib/ants/cli/ancil_create_shapefile.py | 8 +- .../tests/cli/test_ancil_create_shapefile.py | 141 ++++++------------ 2 files changed, 50 insertions(+), 99 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index 5818d5a..f7eba9f 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -13,8 +13,6 @@ Rotated pole domains can be specified using the land sea mask argument, where the longitude, latitude pairs are rotated to the new pole It is assumed that the shapefile is on a standard unrotated lon-lat grid. - -# TODO I need to decide if the central rotated longitude has to be set to 180 """ import argparse import json @@ -25,6 +23,7 @@ import iris.coord_systems import numpy as np from ants.io.load import load_cube, load_landsea_mask +from ants.utils.cube import CubeBuilder from osgeo import ogr from shapely.geometry import Polygon @@ -166,13 +165,14 @@ def _load_polygon_from_json(json_file, target_lsm_path, source_cube_path): if target_lsm_path is not None and source_cube_path is None: crs = iris.coord_systems.GeogCS(6371229.0) - source_cube = ants.utils.cube.CubeBuilder(crs, (2, 2))._cube + source_cube = CubeBuilder(crs, (2, 2))._cube elif target_lsm_path is not None and source_cube_path is not None: source_cube = load_cube(source_cube_path) if target_lsm_path is not None: target_lsm = load_landsea_mask(target_lsm_path) - if not _validate_coord_system(target_lsm): + _validate_coord_system(target_lsm) + if not _validate_orientation(target_lsm): points = _transform_coordinates(target_lsm, source_cube, points) polygon = Polygon(points) diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index 71324ac..a0586b9 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -19,54 +19,10 @@ _validate_orientation, ) from ants.tests.stock import geodetic +from ants.utils.cube import CubeBuilder -def _create_cube(north_pole_lon, north_pole_lat, ellipsoid): - # stock geodetic automatically returns an identity rotated geodetic - # back to geodetic - crs = iris.coord_systems.RotatedGeogCS( - grid_north_pole_latitude=north_pole_lat, - grid_north_pole_longitude=north_pole_lon, - ellipsoid=ellipsoid, - ) - lons, lats = [0.0], [0.0] - lon_coord = iris.coords.DimCoord( - lons, - standard_name="grid_longitude", - units="degrees", - coord_system=crs, - ) - - lat_coord = iris.coords.DimCoord( - lats, - standard_name="grid_latitude", - units="degrees", - coord_system=crs, - ) - - data = np.zeros((len(lats), len(lons))) - - cube = iris.cube.Cube(data, dim_coords_and_dims=[(lat_coord, 0), (lon_coord, 1)]) - return cube - - -class Common(object): - def setUp(self): - # unrotated source cube - self.source_cube = geodetic((2, 2)) - # spherical input coordinate systems are simple to test - self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) - self.sphere_source = geodetic((2, 2), crs=self.sphere_crs) - self.sphere_identity_target_lsm = _create_cube(0.0, 90.0, self.sphere_crs) - self.sphere_equator_target_lsm = geodetic( - (2, 2), north_pole_lat=0.0, north_pole_lon=0.0, crs=self.sphere_crs - ) - self.sphere_rotate_lon = geodetic( - (2, 2), north_pole_lat=90.0, north_pole_lon=90.0, crs=self.sphere_crs - ) - - -class Test__validate_coord_system(Common, ants.tests.TestCase): +class Test__validate_coord_system(ants.tests.TestCase): def test_unrotated_target_lsm(self): target_lsm = geodetic((2, 2)) @@ -81,9 +37,16 @@ def test_unrotated_target_lsm(self): _validate_coord_system(target_lsm) -class Test__validate_orientation(Common, ants.tests.TestCase): +class Test__validate_orientation(ants.tests.TestCase): def setUp(self): - super().setUp() + self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) + self.sphere_identity_crs = iris.coord_systems.RotatedGeogCS(90.0, 0.0) + self.sphere_identity_target_lsm = CubeBuilder( + self.sphere_identity_crs, (2, 2) + )._cube + self.sphere_rotate_lon = geodetic( + (2, 2), north_pole_lat=90.0, north_pole_lon=90.0, crs=self.sphere_crs + ) self.warning_msg = ( "target_lsm has a geodetic coordinate system with pole located" "at grid_longitude=0.0, grid_latitude=90.0." @@ -98,9 +61,20 @@ def test_no_warning_raised(self): self.assertFalse(_validate_orientation(self.sphere_rotate_lon)) -class Test__transform_coordinates(Common, ants.tests.TestCase): +class Test__transform_coordinates(ants.tests.TestCase): def setUp(self): - super().setUp() + self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) + self.sphere_identity_crs = iris.coord_systems.RotatedGeogCS(90.0, 0.0) + self.sphere_identity_target_lsm = CubeBuilder( + self.sphere_identity_crs, (2, 2) + )._cube + self.sphere_source = geodetic((2, 2), crs=self.sphere_crs) + self.sphere_equator_target_lsm = geodetic( + (2, 2), north_pole_lat=0.0, north_pole_lon=0.0, crs=self.sphere_crs + ) + self.sphere_rotate_lon = geodetic( + (2, 2), north_pole_lat=90.0, north_pole_lon=90.0, crs=self.sphere_crs + ) self.vary_latitudes = np.array( [[0.0, 0.0, 0.0, 0.0, 0.0], [90.0, 45.0, 0.0, -45.0, -90.0]] ).T @@ -126,8 +100,7 @@ def test_identity_rotation_sphere(self): rotated_coords = _transform_coordinates( self.sphere_identity_target_lsm, self.sphere_source, self.vary_latitudes ) - - assert ants.utils.ndarray.allclose( + self.assertArrayAlmostEqual( expected_rotation[check_mask], rotated_coords[check_mask] ) @@ -145,7 +118,7 @@ def test_longitudinal_rotation(self): rotated_coords = _transform_coordinates( self.sphere_rotate_lon, self.sphere_source, points ) - assert ants.utils.ndarray.allclose( + self.assertArrayAlmostEqual( expected_rotation[check_mask], rotated_coords[check_mask] ) @@ -158,7 +131,7 @@ def test_rotation_to_equator_sphere(self): rotated_coords = _transform_coordinates( self.sphere_equator_target_lsm, self.sphere_source, self.vary_longitudes ) - assert ants.utils.ndarray.allclose( + self.assertArrayAlmostEqual( expected_rotation[check_mask], rotated_coords[check_mask] ) @@ -180,13 +153,17 @@ def test_negative_longitudes_converted(self): target_crs.transform_points.return_value = neg_lons result = _transform_coordinates(target_lsm, source_cube, neg_lons) - assert ants.utils.ndarray.allclose(true_coords, result) + self.assertArrayEqual(true_coords, result) -class Test__load_polygon_from_json(Common, ants.tests.TestCase): +class Test__load_polygon_from_json(ants.tests.TestCase): def setUp(self): - super().setUp() self.json_values = [[1, 2], [3, 4], [5, 6], [7, 8]] + self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) + self.sphere_source = geodetic((2, 2), crs=self.sphere_crs) + self.sphere_equator_target_lsm = geodetic( + (2, 2), north_pole_lat=0.0, north_pole_lon=0.0, crs=self.sphere_crs + ) def patch_loader(func): @mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") @@ -200,7 +177,7 @@ def wrapper( mock_json, mock_open, mock_transform, - mock_validate, + mock_validate_coord, mock_load, mock_load_lsm, *args, @@ -211,7 +188,7 @@ def wrapper( mock_json, mock_open, mock_transform, - mock_validate, + mock_validate_coord, mock_load, mock_load_lsm, *args, @@ -226,7 +203,7 @@ def test_load_call( mock_json, mock_open, mock_transform, - mock_validate, + mock_validate_coord, mock_load, mock_load_lsm, ): @@ -237,76 +214,50 @@ def test_load_call( mock_json.assert_called_once() mock_load_lsm.assert_not_called() mock_load.assert_not_called() - mock_validate.assert_not_called() mock_transform.assert_not_called() @patch_loader - def test_transform_not_called( + def test_transform_called( self, mock_json, mock_open, mock_transform, - mock_validate, + mock_validate_coord, mock_load, mock_load_lsm, ): mock_json.return_value = self.json_values - mock_validate.return_value = True - mock_load_lsm.return_value = Mock() + mock_load_lsm.return_value = self.sphere_equator_target_lsm mock_load.return_value = Mock() + mock_transform.return_value = self.json_values _ = _load_polygon_from_json("json/path", "lsm/path", "source_cube/path") mock_load_lsm.assert_called_once() mock_load.assert_called_once() - mock_transform.assert_not_called() + mock_transform.assert_called_once() @patch_loader - def test_transform_called( + def test_default_source( self, mock_json, mock_open, mock_transform, - mock_validate, + mock_validate_coord, mock_load, mock_load_lsm, ): mock_json.return_value = self.json_values - mock_validate.return_value = False - mock_load_lsm.return_value = Mock() + mock_load_lsm.return_value = self.sphere_equator_target_lsm mock_load.return_value = Mock() mock_transform.return_value = self.json_values - _ = _load_polygon_from_json("json/path", "lsm/path", "source_cube/path") - - mock_load_lsm.assert_called_once() - mock_load.assert_called_once() - mock_transform.assert_called_once() - - """ - I need to figure out how to get this patch to work - @patch_loader - def test_load_call(self, mock_json, mock_open, mock_transform, - mock_validate, mock_load, mock_load_lsm): - - ants_cube = ants.utils.cube.CubeBuilder(self.sphere_crs,(2,2))._cube - mock_json.return_value = self.json_values - mock_validate.return_value = False - mock_load_lsm.return_value = Mock() - mock_transform.return_value = self.json_values - _ = _load_polygon_from_json("json/path", "lsm/path", None) - mock_json.assert_called_once() mock_load_lsm.assert_called_once() - mock_load.assert_not_called() - mock_validate.assert_called_once() - mock_transform.assert_called_once_with(mock_load_lsm.return_value, - ants_cube, - np.array(self.json_values)) - """ + mock_transform.assert_called_once() class Test__validate_args(ants.tests.TestCase): From 5eed8ddc0351b3e8cfe04ced9dd39ae6a6ed3d52 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Fri, 17 Jul 2026 13:15:17 +0100 Subject: [PATCH 07/29] Remove some comments that are not needed --- lib/ants/cli/ancil_create_shapefile.py | 48 +++++++++++-------- .../tests/cli/test_ancil_create_shapefile.py | 11 ++++- 2 files changed, 37 insertions(+), 22 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index f7eba9f..d43a176 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -12,7 +12,9 @@ Rotated pole domains can be specified using the land sea mask argument, where the longitude, latitude pairs are rotated to the new pole -It is assumed that the shapefile is on a standard unrotated lon-lat grid. +Unless a source cube is provided with the source coordinate reference system, +it assumed that the point defined in the json file are on a standard spherical +unrotated geodetic coordinate reference system. """ import argparse import json @@ -76,8 +78,11 @@ def _validate_orientation(target_lsm): grid_lon = target_crs.grid_north_pole_longitude grid_lat = target_crs.grid_north_pole_latitude + pole_lon_rotation = target_crs.north_pole_grid_longitude - check_crs = ants.utils.ndarray.allclose([grid_lon, grid_lat], [0.0, 90.0]) + check_crs = ants.utils.ndarray.allclose( + [grid_lon, grid_lat, pole_lon_rotation], [0.0, 90.0, 0.0] + ) if check_crs: warnings.warn( "target_lsm has a geodetic coordinate system with pole located" @@ -87,14 +92,25 @@ def _validate_orientation(target_lsm): return check_crs +def _validate_ellipse(target_lsm, source_cube): + # what is the right thing to do here + pass + + +def _validate_args(target_lsm_path, source_cube_path): + if source_cube_path is not None and target_lsm_path is None: + raise ValueError( + "If --source-cube is passed then --target-lsm must" "also be given." + ) + + def _transform_coordinates(target_lsm, source_cube, points): """ Transforms the longitude, latitude points in the source coordinate system to the rotated pole coordinate system defined by target_lsm. The source coordinate system is assumed to be unrotated geodetic - defined on a sphere. Does nothing if the target coordinate system - is unrotated. + defined on a sphere, unless otherwise specified. Parameters ---------- @@ -112,10 +128,6 @@ def _transform_coordinates(target_lsm, source_cube, points): An (m,2) sized numpy array of rotated longitude and latitude points. """ - # sphere = ccrs.Globe(semimajor_axis=6371000.0, semiminor_axis=6371000.0) - - # source_crs = ccrs.Geodetic(globe=sphere) - source_crs = source_cube.coord_system().as_cartopy_crs() target_coord = target_lsm.coord_system() target_crs = target_coord.as_cartopy_crs() @@ -124,7 +136,7 @@ def _transform_coordinates(target_lsm, source_cube, points): source_crs, points[:, 0], points[:, 1] )[:, :2] - # Enforce that the points are all measured in the positive direction + # Enforce that longitude is measured in the positive direction negative_rotated_lons = np.where(rotated_points[:, 0] < 0) rotated_points[negative_rotated_lons, 0] += 360.0 @@ -144,6 +156,10 @@ def _load_polygon_from_json(json_file, target_lsm_path, source_cube_path): Load a json file containing a list of pairs of longitude, latitude points to create a polygon from. + If a valid rotated pole is specifed in the target lsm, first transforms + the points from the source coordinate system to the rotated pole before + creating a polygon. + Parameters ---------- json_file : str @@ -188,8 +204,8 @@ def main(json_file, output, target_lsm_path, source_cube_path): If target_lsm_path is provided, the points are first transformed from an unrotated geodetic coordinate system to a rotated pole coordinate system - specified by the lsm. It is assumed that the points specified in the json - file are on an unrotated geodetic grid. + specified by the lsm. It is assumed that the points in the json + file are on an unrotated geodetic grid, unless a crs is provided. Parameters ---------- @@ -255,19 +271,11 @@ def _get_parser(): return parser -def _validate_args(args): - - if args.source_cube is not None and args.target_lsm is None: - raise ValueError( - "If --source-cube is passed then --target-lsm must" "also be given." - ) - - def cli_interface(): parser = _get_parser() args = parser.parse_args() - _validate_args(args) + _validate_args(args.target_lsm, args.source_cube) main(args.json_file, args.output, args.target_lsm, args.source_cube) diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index a0586b9..d496aba 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -41,9 +41,13 @@ class Test__validate_orientation(ants.tests.TestCase): def setUp(self): self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) self.sphere_identity_crs = iris.coord_systems.RotatedGeogCS(90.0, 0.0) + self.sphere_rotated_crs = iris.coord_systems.RotatedGeogCS(90.0, 0.0, 180.0) self.sphere_identity_target_lsm = CubeBuilder( self.sphere_identity_crs, (2, 2) )._cube + self.sphere_rotated_target_lsm = CubeBuilder( + self.sphere_rotated_crs, (2, 2) + )._cube self.sphere_rotate_lon = geodetic( (2, 2), north_pole_lat=90.0, north_pole_lon=90.0, crs=self.sphere_crs ) @@ -60,6 +64,10 @@ def test_warning_raised(self): def test_no_warning_raised(self): self.assertFalse(_validate_orientation(self.sphere_rotate_lon)) + # an optional rotation can be applied afterward + def test_no_warning_rotated(self): + self.assertFalse(_validate_orientation(self.sphere_rotated_target_lsm)) + class Test__transform_coordinates(ants.tests.TestCase): def setUp(self): @@ -138,7 +146,6 @@ def test_rotation_to_equator_sphere(self): def test_negative_longitudes_converted(self): true_coords = np.array([[90.0, 160.0, 326.7], [0.0, 0.0, 0.0]]).T - # presumably I can't get a -400.0 for example neg_lons = np.array([[-270.0, -200.0, -33.3], [0.0, 0.0, 0.0]]).T target_lsm, source_cube = Mock(), Mock() @@ -271,4 +278,4 @@ def test_error_raised(self): ) error_msg = "If --source-cube is passed then --target-lsm must" "also be given." with self.assertRaisesRegex(ValueError, error_msg): - _validate_args(args) + _validate_args(args.target_lsm, args.source_cube) From 4448d869cfdc1dd0c159a4632375b9ac8aa1a4ac Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Fri, 17 Jul 2026 15:06:11 +0100 Subject: [PATCH 08/29] Improved comments in tests. --- lib/ants/cli/ancil_create_shapefile.py | 5 --- .../tests/cli/test_ancil_create_shapefile.py | 40 ++++++++++++------- 2 files changed, 25 insertions(+), 20 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index d43a176..6f13c54 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -92,11 +92,6 @@ def _validate_orientation(target_lsm): return check_crs -def _validate_ellipse(target_lsm, source_cube): - # what is the right thing to do here - pass - - def _validate_args(target_lsm_path, source_cube_path): if source_cube_path is not None and target_lsm_path is None: raise ValueError( diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index d496aba..f7c4d36 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -40,8 +40,12 @@ def test_unrotated_target_lsm(self): class Test__validate_orientation(ants.tests.TestCase): def setUp(self): self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) - self.sphere_identity_crs = iris.coord_systems.RotatedGeogCS(90.0, 0.0) - self.sphere_rotated_crs = iris.coord_systems.RotatedGeogCS(90.0, 0.0, 180.0) + self.sphere_identity_crs = iris.coord_systems.RotatedGeogCS( + 90.0, 0.0, ellipsoid=self.sphere_crs + ) + self.sphere_rotated_crs = iris.coord_systems.RotatedGeogCS( + 90.0, 0.0, 180.0, ellipsoid=self.sphere_crs + ) self.sphere_identity_target_lsm = CubeBuilder( self.sphere_identity_crs, (2, 2) )._cube @@ -64,15 +68,17 @@ def test_warning_raised(self): def test_no_warning_raised(self): self.assertFalse(_validate_orientation(self.sphere_rotate_lon)) - # an optional rotation can be applied afterward - def test_no_warning_rotated(self): + def test_no_warning_raised_rotated(self): + # The pole can be optionally rotated in longitude after relocation. self.assertFalse(_validate_orientation(self.sphere_rotated_target_lsm)) class Test__transform_coordinates(ants.tests.TestCase): def setUp(self): self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) - self.sphere_identity_crs = iris.coord_systems.RotatedGeogCS(90.0, 0.0) + self.sphere_identity_crs = iris.coord_systems.RotatedGeogCS( + 90.0, 0.0, ellipsoid=self.sphere_crs + ) self.sphere_identity_target_lsm = CubeBuilder( self.sphere_identity_crs, (2, 2) )._cube @@ -94,15 +100,15 @@ def setUp(self): ).T def test_identity_rotation_sphere(self): + # Rotation to pole at (0.0, 90.0). + # Rotated geodesic by convention defines the new pole 180 degrees. + # from passed longitude. true_lats = np.array([90, 45, 0, -45, -90]) - # The rotated geodesic system flips the longitude - # Rotated geodesic has some odd properties, like the central - # longitude being set to 180.0 true_lons = np.array([np.nan, 180.0, 180.0, 180.0, 180.0]) - expected_rotation = np.array([true_lons, true_lats]).T - # we cannot check the longitude at the pole (any longitude is valid) + + # We cannot check the longitude at the pole (any longitude is valid). check_mask = ~np.isnan(expected_rotation) rotated_coords = _transform_coordinates( @@ -113,14 +119,14 @@ def test_identity_rotation_sphere(self): ) def test_longitudinal_rotation(self): - # these tests can be counter-intuitive as the prime-meridian - # is by default at 180.0 so a rotation of 90.0 causes the prime meridian - # to be at 270.0 (it is always 180 further round from the longitude given) - points = np.array([[0.0, 45.0, 90.0, 135.0], [0.0, 0.0, 0.0, 0.0]]).T + # Rotation to pole (90.0, 90.0). + # The new pole is located at 270.0 after a 90.0 rotation. + points = np.array([[0.0, 45.0, 90.0, 135.0], [0.0, 0.0, 0.0, 0.0]]).T true_lats = np.array([0.0, 0.0, np.nan, 0.0]) true_lons = np.array([90.0, 135.0, 180.0, 225.0]) expected_rotation = np.array([true_lons, true_lats]).T + check_mask = ~np.isnan(expected_rotation) rotated_coords = _transform_coordinates( @@ -131,8 +137,10 @@ def test_longitudinal_rotation(self): ) def test_rotation_to_equator_sphere(self): + # Rotation to pole at (0.0, 0.0). + # Latitudes become longitudes and vice versa. + true_lats = np.array([-45.0, 0.0, 45.0, 90.0, 45.0, 0.0, -45.0]) - # Again the results for lon depend on where the central lon is true_lons = np.array([90.0, 90.0, 90.0, np.nan, 270.0, 270.0, 270.0]) expected_rotation = np.array([true_lons, true_lats]).T check_mask = ~np.isnan(expected_rotation) @@ -144,6 +152,8 @@ def test_rotation_to_equator_sphere(self): ) def test_negative_longitudes_converted(self): + # Check that negative longitudes are converted to positive values. + # No transform performed here. true_coords = np.array([[90.0, 160.0, 326.7], [0.0, 0.0, 0.0]]).T neg_lons = np.array([[-270.0, -200.0, -33.3], [0.0, 0.0, 0.0]]).T From 5d1a68aa38243ad4b6a1fc22881da5cb2561ac4f Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Fri, 17 Jul 2026 16:48:48 +0100 Subject: [PATCH 09/29] Fixed issue causing polygon to invert --- lib/ants/cli/ancil_create_shapefile.py | 4 --- .../tests/cli/test_ancil_create_shapefile.py | 28 +++---------------- 2 files changed, 4 insertions(+), 28 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index 6f13c54..dbef930 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -131,10 +131,6 @@ def _transform_coordinates(target_lsm, source_cube, points): source_crs, points[:, 0], points[:, 1] )[:, :2] - # Enforce that longitude is measured in the positive direction - negative_rotated_lons = np.where(rotated_points[:, 0] < 0) - rotated_points[negative_rotated_lons, 0] += 360.0 - _LOGGER.info( "Input json file transformed to new pole rotated coordinate system at" "pole longitude=%s, pole latitude=%s, central rotated longitude=%s.", diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index f7c4d36..d617084 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -105,7 +105,7 @@ def test_identity_rotation_sphere(self): # from passed longitude. true_lats = np.array([90, 45, 0, -45, -90]) - true_lons = np.array([np.nan, 180.0, 180.0, 180.0, 180.0]) + true_lons = np.array([np.nan, -180.0, -180.0, -180.0, -180.0]) expected_rotation = np.array([true_lons, true_lats]).T # We cannot check the longitude at the pole (any longitude is valid). @@ -124,7 +124,7 @@ def test_longitudinal_rotation(self): points = np.array([[0.0, 45.0, 90.0, 135.0], [0.0, 0.0, 0.0, 0.0]]).T true_lats = np.array([0.0, 0.0, np.nan, 0.0]) - true_lons = np.array([90.0, 135.0, 180.0, 225.0]) + true_lons = np.array([90.0, 135.0, -180.0, -135.0]) expected_rotation = np.array([true_lons, true_lats]).T check_mask = ~np.isnan(expected_rotation) @@ -132,6 +132,7 @@ def test_longitudinal_rotation(self): rotated_coords = _transform_coordinates( self.sphere_rotate_lon, self.sphere_source, points ) + self.assertArrayAlmostEqual( expected_rotation[check_mask], rotated_coords[check_mask] ) @@ -141,7 +142,7 @@ def test_rotation_to_equator_sphere(self): # Latitudes become longitudes and vice versa. true_lats = np.array([-45.0, 0.0, 45.0, 90.0, 45.0, 0.0, -45.0]) - true_lons = np.array([90.0, 90.0, 90.0, np.nan, 270.0, 270.0, 270.0]) + true_lons = np.array([90.0, 90.0, 90.0, np.nan, -90.0, -90.0, -90.0]) expected_rotation = np.array([true_lons, true_lats]).T check_mask = ~np.isnan(expected_rotation) rotated_coords = _transform_coordinates( @@ -151,27 +152,6 @@ def test_rotation_to_equator_sphere(self): expected_rotation[check_mask], rotated_coords[check_mask] ) - def test_negative_longitudes_converted(self): - # Check that negative longitudes are converted to positive values. - # No transform performed here. - - true_coords = np.array([[90.0, 160.0, 326.7], [0.0, 0.0, 0.0]]).T - neg_lons = np.array([[-270.0, -200.0, -33.3], [0.0, 0.0, 0.0]]).T - - target_lsm, source_cube = Mock(), Mock() - source_crs, target_crs = Mock(), Mock() - target_coord = Mock() - - source_cube.coord_system.return_value.as_cartopy_crs.return_value = source_crs - target_lsm.coord_system.return_value = target_coord - - target_coord.as_cartopy_crs.return_value = target_crs - - target_crs.transform_points.return_value = neg_lons - - result = _transform_coordinates(target_lsm, source_cube, neg_lons) - self.assertArrayEqual(true_coords, result) - class Test__load_polygon_from_json(ants.tests.TestCase): def setUp(self): From af32edfa75a2ad9d1c61f4acaad3a64d31397894 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Mon, 20 Jul 2026 13:10:37 +0100 Subject: [PATCH 10/29] Moved functionality out of load polygon to ease testing. --- lib/ants/cli/ancil_create_shapefile.py | 82 +++++-- .../tests/cli/test_ancil_create_shapefile.py | 201 +++++++++--------- 2 files changed, 167 insertions(+), 116 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index dbef930..5f4b542 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -32,7 +32,7 @@ _LOGGER = logging.getLogger(__name__) -def _validate_coord_system(target_lsm): +def _check_coord_system_type(target_lsm): """ Check that target_lsm has a rotated pole coordinate system. @@ -61,17 +61,23 @@ def _validate_coord_system(target_lsm): def _validate_orientation(target_lsm): """ - Check if target_lsm is a non-rotated pole. + Check if target_lsm has a valid rotated pole coordinate system. Parameters ---------- target_lsm : :class:`iris.cube.Cube` - The lsm cube specifiying the rotated pole coordinate system. + The lsm cube specifying the rotated pole coordinate system. + + Returns + ------- + valid_crs : bool + True if target_lsm is a rotated pole, + False otherwise. Warns ------ - UserWarning : - If target_lsm is already at latitude=90.0, longitude=0.0. + UserWarning + If target_lsm is already at latitude=90.0 and longitude=0.0. """ target_crs = target_lsm.coord_system() @@ -80,22 +86,24 @@ def _validate_orientation(target_lsm): grid_lat = target_crs.grid_north_pole_latitude pole_lon_rotation = target_crs.north_pole_grid_longitude - check_crs = ants.utils.ndarray.allclose( - [grid_lon, grid_lat, pole_lon_rotation], [0.0, 90.0, 0.0] + valid_crs = not ( + ants.utils.ndarray.allclose( + [grid_lon, grid_lat, pole_lon_rotation], [0.0, 90.0, 0.0] + ) ) - if check_crs: + if not valid_crs: warnings.warn( "target_lsm has a geodetic coordinate system with pole located" f"at grid_longitude={grid_lon}, grid_latitude={grid_lat}." "No transformation will be carried out." ) - return check_crs + return valid_crs def _validate_args(target_lsm_path, source_cube_path): if source_cube_path is not None and target_lsm_path is None: raise ValueError( - "If --source-cube is passed then --target-lsm must" "also be given." + "If --source-cube is passed then --target-lsm must also be given." ) @@ -120,7 +128,7 @@ def _transform_coordinates(target_lsm, source_cube, points): Returns ------- rotated_points : :class:`np.ndarry` - An (m,2) sized numpy array of rotated longitude and latitude points. + An (m,2) sized numpy array of transformed longitude and latitude points. """ source_crs = source_cube.coord_system().as_cartopy_crs() @@ -142,6 +150,47 @@ def _transform_coordinates(target_lsm, source_cube, points): return rotated_points +def _transform_if_required(target_lsm_path, source_cube_path, points): + """ + Load the target lsm and create a source cube if not provided. The points + are then transformed to the specified rotated pole coordinate system + if the provided coordinate system is a rotated pole and is not + coincident with the north pole. + + Parameters + ---------- + target_lsm_path : str + File path for a land sea mask that provides the new rotated pole + coordinates to which the longitude, latitude pairs will be transformed. + source_cube_path : str + File path to an iris cube specifying the co-ordinate system of the + input json file. + points : :class:`np.ndarry` + An (m,2) sized numpy array of longitude and latitude points. + + Returns + ------- + rotated_points : :class:`np.ndarry` + An (m,2) sized numpy array of transformed longitude and latitude points. + """ + + target_lsm = load_landsea_mask(target_lsm_path) + _check_coord_system_type(target_lsm) + + if source_cube_path is None: + crs = iris.coord_systems.GeogCS(6371229.0) + source_cube = CubeBuilder(crs, (2, 2))._cube + else: + source_cube = load_cube(source_cube_path) + + if _validate_orientation(target_lsm): + rotated_points = _transform_coordinates(target_lsm, source_cube, points) + else: + rotated_points = np.copy(points) + + return rotated_points + + def _load_polygon_from_json(json_file, target_lsm_path, source_cube_path): """ Load a json file containing a list of pairs of longitude, latitude points @@ -170,17 +219,8 @@ def _load_polygon_from_json(json_file, target_lsm_path, source_cube_path): points = json.load(polygon_json) points = np.array(points) - if target_lsm_path is not None and source_cube_path is None: - crs = iris.coord_systems.GeogCS(6371229.0) - source_cube = CubeBuilder(crs, (2, 2))._cube - elif target_lsm_path is not None and source_cube_path is not None: - source_cube = load_cube(source_cube_path) - if target_lsm_path is not None: - target_lsm = load_landsea_mask(target_lsm_path) - _validate_coord_system(target_lsm) - if not _validate_orientation(target_lsm): - points = _transform_coordinates(target_lsm, source_cube, points) + points = _transform_if_required(target_lsm_path, source_cube_path, points) polygon = Polygon(points) diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index d617084..64476fa 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -12,17 +12,18 @@ import iris import numpy as np from ants.cli.ancil_create_shapefile import ( + _check_coord_system_type, _load_polygon_from_json, _transform_coordinates, + _transform_if_required, _validate_args, - _validate_coord_system, _validate_orientation, ) from ants.tests.stock import geodetic from ants.utils.cube import CubeBuilder -class Test__validate_coord_system(ants.tests.TestCase): +class Test__check_coord_system_type(ants.tests.TestCase): def test_unrotated_target_lsm(self): target_lsm = geodetic((2, 2)) @@ -34,7 +35,7 @@ def test_unrotated_target_lsm(self): ) with self.assertRaisesRegex(ValueError, error_msg): - _validate_coord_system(target_lsm) + _check_coord_system_type(target_lsm) class Test__validate_orientation(ants.tests.TestCase): @@ -66,11 +67,11 @@ def test_warning_raised(self): _validate_orientation(self.sphere_identity_target_lsm) def test_no_warning_raised(self): - self.assertFalse(_validate_orientation(self.sphere_rotate_lon)) + self.assertTrue(_validate_orientation(self.sphere_rotate_lon)) def test_no_warning_raised_rotated(self): # The pole can be optionally rotated in longitude after relocation. - self.assertFalse(_validate_orientation(self.sphere_rotated_target_lsm)) + self.assertTrue(_validate_orientation(self.sphere_rotated_target_lsm)) class Test__transform_coordinates(ants.tests.TestCase): @@ -153,119 +154,129 @@ def test_rotation_to_equator_sphere(self): ) -class Test__load_polygon_from_json(ants.tests.TestCase): +class Test__validate_args(ants.tests.TestCase): + + def test_error_raised(self): + args = argparse.Namespace( + json_file="json/path", + output="output/path", + target_lsm=None, + source_cube="source/path", + ) + error_msg = "If --source-cube is passed then --target-lsm must also be given." + with self.assertRaisesRegex(ValueError, error_msg): + _validate_args(args.target_lsm, args.source_cube) + + +class Test__transform_if_required(ants.tests.TestCase): def setUp(self): - self.json_values = [[1, 2], [3, 4], [5, 6], [7, 8]] + self.points = np.array([[1, 2], [3, 4], [5, 6], [7, 8]]) self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) self.sphere_source = geodetic((2, 2), crs=self.sphere_crs) self.sphere_equator_target_lsm = geodetic( (2, 2), north_pole_lat=0.0, north_pole_lon=0.0, crs=self.sphere_crs ) + self.sphere_identity_crs = iris.coord_systems.RotatedGeogCS( + 90.0, 0.0, ellipsoid=self.sphere_crs + ) + self.sphere_identity_target_lsm = CubeBuilder( + self.sphere_identity_crs, (2, 2) + )._cube - def patch_loader(func): - @mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") - @mock.patch("ants.cli.ancil_create_shapefile.load_cube") - @mock.patch("ants.cli.ancil_create_shapefile._validate_coord_system") - @mock.patch("ants.cli.ancil_create_shapefile._transform_coordinates") - @mock.patch("builtins.open", new_callable=mock.mock_open) - @mock.patch("ants.cli.ancil_create_shapefile.json.load") - def wrapper( - self, - mock_json, - mock_open, - mock_transform, - mock_validate_coord, - mock_load, - mock_load_lsm, - *args, - **kwargs, - ): - return func( - self, - mock_json, - mock_open, - mock_transform, - mock_validate_coord, - mock_load, - mock_load_lsm, - *args, - **kwargs, - ) - - return wrapper - - @patch_loader - def test_load_call( - self, - mock_json, - mock_open, - mock_transform, - mock_validate_coord, - mock_load, - mock_load_lsm, + @mock.patch("ants.cli.ancil_create_shapefile._validate_orientation") + @mock.patch("ants.cli.ancil_create_shapefile._transform_coordinates") + @mock.patch("ants.cli.ancil_create_shapefile.CubeBuilder") + @mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") + def test_transform_coords_not_called( + self, mock_lsm, mock_source, mock_transform_coord, mock_validate ): - mock_json.return_value = self.json_values - _ = _load_polygon_from_json("json/path", None, None) + mock_lsm.return_value = self.sphere_identity_target_lsm + mock_validate.return_value = False - mock_json.assert_called_once() - mock_load_lsm.assert_not_called() - mock_load.assert_not_called() - mock_transform.assert_not_called() + rotated_points = _transform_if_required("lsm/path", None, self.points) + + mock_lsm.assert_called_once() + mock_source.assert_called_once() + mock_transform_coord.assert_not_called() + self.assertArrayEqual(rotated_points, self.points) + + @mock.patch("ants.cli.ancil_create_shapefile._transform_coordinates") + @mock.patch("ants.cli.ancil_create_shapefile.CubeBuilder") + @mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") + def test_transform_called(self, mock_lsm, mock_source, mock_transform_coord): + + mock_lsm.return_value = self.sphere_equator_target_lsm + + _ = _transform_if_required("lsm/path", None, self.points) + + mock_lsm.assert_called_once() + mock_source.assert_called_once() + mock_transform_coord.assert_called_once() - @patch_loader - def test_transform_called( - self, - mock_json, - mock_open, - mock_transform, - mock_validate_coord, - mock_load, - mock_load_lsm, + @mock.patch("ants.cli.ancil_create_shapefile.load_cube") + @mock.patch("ants.cli.ancil_create_shapefile._transform_coordinates") + @mock.patch("ants.cli.ancil_create_shapefile.CubeBuilder") + @mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") + def test_source_cube_loaded( + self, mock_lsm, mock_source, mock_transform_coord, mock_load_cube ): + mock_lsm.return_value = self.sphere_equator_target_lsm + mock_load_cube.return_value = Mock() + + _ = _transform_if_required("lsm/path", "source/path", self.points) + + mock_lsm.assert_called_once() + mock_source.assert_not_called() + mock_load_cube.assert_called_once() + mock_transform_coord.assert_called_once() + + +class Test__load_polygon_from_json(ants.tests.TestCase): + def setUp(self): + self.json_values = [[1, 2], [3, 4], [5, 6], [7, 8]] + + @mock.patch("ants.cli.ancil_create_shapefile._transform_if_required") + @mock.patch("builtins.open", new_callable=mock.mock_open) + @mock.patch("ants.cli.ancil_create_shapefile.json.load") + def test_json_load(self, mock_json, mock_open, mock_transform): + mock_json.return_value = self.json_values - mock_load_lsm.return_value = self.sphere_equator_target_lsm - mock_load.return_value = Mock() - mock_transform.return_value = self.json_values - _ = _load_polygon_from_json("json/path", "lsm/path", "source_cube/path") + _ = _load_polygon_from_json("json/path", None, None) - mock_load_lsm.assert_called_once() - mock_load.assert_called_once() - mock_transform.assert_called_once() + mock_json.assert_called_once() + mock_transform.assert_not_called() - @patch_loader - def test_default_source( - self, - mock_json, - mock_open, - mock_transform, - mock_validate_coord, - mock_load, - mock_load_lsm, - ): + @mock.patch("ants.cli.ancil_create_shapefile._transform_if_required") + @mock.patch("builtins.open", new_callable=mock.mock_open) + @mock.patch("ants.cli.ancil_create_shapefile.json.load") + def test_transform_called(self, mock_json, mock_open, mock_transform): mock_json.return_value = self.json_values - mock_load_lsm.return_value = self.sphere_equator_target_lsm - mock_load.return_value = Mock() - mock_transform.return_value = self.json_values _ = _load_polygon_from_json("json/path", "lsm/path", None) - mock_load_lsm.assert_called_once() + target_path, source_path, points = mock_transform.call_args.args + mock_json.assert_called_once() mock_transform.assert_called_once() + self.assertArrayEqual(points, self.json_values) + self.assertEqual("lsm/path", target_path) + self.assertEqual(None, source_path) + @mock.patch("ants.cli.ancil_create_shapefile._transform_if_required") + @mock.patch("builtins.open", new_callable=mock.mock_open) + @mock.patch("ants.cli.ancil_create_shapefile.json.load") + def test_transform_called_with_source(self, mock_json, mock_open, mock_transform): -class Test__validate_args(ants.tests.TestCase): + mock_json.return_value = self.json_values - def test_error_raised(self): - args = argparse.Namespace( - json_file="json/path", - output="output/path", - target_lsm=None, - source_cube="source/path", - ) - error_msg = "If --source-cube is passed then --target-lsm must" "also be given." - with self.assertRaisesRegex(ValueError, error_msg): - _validate_args(args.target_lsm, args.source_cube) + _ = _load_polygon_from_json("json/path", "lsm/path", "source/path") + + target_path, source_path, points = mock_transform.call_args.args + mock_json.assert_called_once() + mock_transform.assert_called_once() + self.assertArrayEqual(points, self.json_values) + self.assertEqual("lsm/path", target_path) + self.assertEqual("source/path", source_path) From fa4bbebed6971db40c97a6ff705de9123fa12307 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Mon, 20 Jul 2026 13:45:53 +0100 Subject: [PATCH 11/29] Minor reformatting of docstrings. --- lib/ants/cli/ancil_create_shapefile.py | 29 ++++++++++--------- .../tests/cli/test_ancil_create_shapefile.py | 2 -- .../app/ancil_create_shapefile/rose-app.conf | 3 +- 3 files changed, 16 insertions(+), 18 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index 5f4b542..c903bee 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -11,9 +11,9 @@ points defining a single polygon in a specified polygon file. Rotated pole domains can be specified using the land sea mask argument, -where the longitude, latitude pairs are rotated to the new pole +where the longitude, latitude pairs are rotated to the new pole. Unless a source cube is provided with the source coordinate reference system, -it assumed that the point defined in the json file are on a standard spherical +it assumed that the points defined in the json file are on a standard spherical unrotated geodetic coordinate reference system. """ import argparse @@ -61,12 +61,13 @@ def _check_coord_system_type(target_lsm): def _validate_orientation(target_lsm): """ - Check if target_lsm has a valid rotated pole coordinate system. + Check if the rotated pole coordinate system specified by the target lsm + is coincident with the north pole. Parameters ---------- target_lsm : :class:`iris.cube.Cube` - The lsm cube specifying the rotated pole coordinate system. + The target lsm cube specifying the rotated pole coordinate system. Returns ------- @@ -120,15 +121,15 @@ def _transform_coordinates(target_lsm, source_cube, points): target_lsm : :class:`iris.cube.Cube` The lsm cube specifiying the rotated pole coordinate system. source_cube : :class:`iris.cube.Cube` - An iris cube specifying the co-ordinate system of the + An iris cube specifying the coordinate system of the input json file. points : :class:`np.ndarry` - An (m,2) sized numpy array of m longitude and m latitude points. + An (m,2) shaped numpy array of m longitude and m latitude points. Returns ------- rotated_points : :class:`np.ndarry` - An (m,2) sized numpy array of transformed longitude and latitude points. + An (m,2) shaped numpy array of transformed longitude and latitude points. """ source_crs = source_cube.coord_system().as_cartopy_crs() @@ -163,15 +164,15 @@ def _transform_if_required(target_lsm_path, source_cube_path, points): File path for a land sea mask that provides the new rotated pole coordinates to which the longitude, latitude pairs will be transformed. source_cube_path : str - File path to an iris cube specifying the co-ordinate system of the + File path to an iris cube specifying the coordinate system of the input json file. points : :class:`np.ndarry` - An (m,2) sized numpy array of longitude and latitude points. + An (m,2) shaped numpy array of longitude and latitude points. Returns ------- rotated_points : :class:`np.ndarry` - An (m,2) sized numpy array of transformed longitude and latitude points. + An (m,2) shaped numpy array of transformed longitude and latitude points. """ target_lsm = load_landsea_mask(target_lsm_path) @@ -208,7 +209,7 @@ def _load_polygon_from_json(json_file, target_lsm_path, source_cube_path): File path for a land sea mask that provides the new rotated pole coordinates to which the longitude, latitude pairs will be transformed. source_cube_path : str - File path to an iris cube specifying the co-ordinate system of the + File path to an iris cube specifying the coordinate system of the input json file. Returns @@ -236,7 +237,7 @@ def main(json_file, output, target_lsm_path, source_cube_path): If target_lsm_path is provided, the points are first transformed from an unrotated geodetic coordinate system to a rotated pole coordinate system specified by the lsm. It is assumed that the points in the json - file are on an unrotated geodetic grid, unless a crs is provided. + file are on an unrotated geodetic grid, unless otherwise specified. Parameters ---------- @@ -248,7 +249,7 @@ def main(json_file, output, target_lsm_path, source_cube_path): File path for a land sea mask that provides the new rotated pole coordinates to which the longitude, latitude pairs will be transformed. source_cube_path : str - File path to an iris cube specifying the co-ordinate system of the + File path to an iris cube specifying the cosordinate system of the input json file. """ @@ -297,7 +298,7 @@ def _get_parser(): type=ants.config.filepath_readable, required=False, help="Path to an iris cube which specifies the co-ordinate" - "system of the json file", + " system of the json file.", ) return parser diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index 64476fa..57befb2 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -6,7 +6,6 @@ import argparse import re from unittest import mock -from unittest.mock import Mock import ants.tests import iris @@ -223,7 +222,6 @@ def test_source_cube_loaded( ): mock_lsm.return_value = self.sphere_equator_target_lsm - mock_load_cube.return_value = Mock() _ = _transform_if_required("lsm/path", "source/path", self.points) diff --git a/rose-stem/app/ancil_create_shapefile/rose-app.conf b/rose-stem/app/ancil_create_shapefile/rose-app.conf index 05f94bd..4f37bbf 100644 --- a/rose-stem/app/ancil_create_shapefile/rose-app.conf +++ b/rose-stem/app/ancil_create_shapefile/rose-app.conf @@ -1,7 +1,6 @@ [command] -default=ants-launch ancil_create_shapefile.py ${SHAPEFILE_JSON} ${OUTPUT} --target-lsm ${LSM} +default=ants-launch ancil_create_shapefile.py ${SHAPEFILE_JSON} ${OUTPUT} [env] OUTPUT=${ROSE_DATA}/ite.shp SHAPEFILE_JSON=${TEST_SOURCES_DIR}/ancil_create_shapefile/ite.json -LSM=/home/users/alasdair.roy/Documents/cartopy_tutorial/target_lsm.nc From 791118e457d1b925ecc93260951c86947b9953ab Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Mon, 20 Jul 2026 14:40:20 +0100 Subject: [PATCH 12/29] Added docstrings to tests. --- .../tests/cli/test_ancil_create_shapefile.py | 49 +++++++++++++++---- 1 file changed, 39 insertions(+), 10 deletions(-) diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index 57befb2..0163ea6 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -25,6 +25,9 @@ class Test__check_coord_system_type(ants.tests.TestCase): def test_unrotated_target_lsm(self): + """Test that an error message is raised if the coordinate system in + target_lsm is not a rotated pole.""" + target_lsm = geodetic((2, 2)) error_msg = re.escape( f"target_lsm.coord_system() {target_lsm.coord_system()} is not" @@ -62,14 +65,22 @@ def setUp(self): ) def test_warning_raised(self): + """Test that a warning is raised with a coordinate system at latitude=90.0, + longitude=0.0 is passed. Check the function returns False.""" + with self.assertWarnsRegex(UserWarning, self.warning_msg): - _validate_orientation(self.sphere_identity_target_lsm) + valid_crs = _validate_orientation(self.sphere_identity_target_lsm) + self.assertFalse(valid_crs) + + def test_true_returned(self): + """Test that True is returned when target_lsm has a valid rotated pole.""" - def test_no_warning_raised(self): self.assertTrue(_validate_orientation(self.sphere_rotate_lon)) def test_no_warning_raised_rotated(self): - # The pole can be optionally rotated in longitude after relocation. + """Test that passing a non-zero central rotated longitude is still + accounted for when checking if the coordinate system is rotated.""" + self.assertTrue(_validate_orientation(self.sphere_rotated_target_lsm)) @@ -100,9 +111,10 @@ def setUp(self): ).T def test_identity_rotation_sphere(self): - # Rotation to pole at (0.0, 90.0). - # Rotated geodesic by convention defines the new pole 180 degrees. - # from passed longitude. + """Test that a rotation to a pole at latitude=90.0, longitude=0.0 + returns the same points, except a 180.0 degree rotation in longitude. By + convention the new pole is defined 180.0 rotated from the provided + longitude.""" true_lats = np.array([90, 45, 0, -45, -90]) true_lons = np.array([np.nan, -180.0, -180.0, -180.0, -180.0]) @@ -119,8 +131,9 @@ def test_identity_rotation_sphere(self): ) def test_longitudinal_rotation(self): - # Rotation to pole (90.0, 90.0). - # The new pole is located at 270.0 after a 90.0 rotation. + """Test rotation to a new pole at latitude=90.0, longitude=90.0. + The new pole is located at latitude=90.0, longitude 270.0 in the rotated + pole coordinate system.""" points = np.array([[0.0, 45.0, 90.0, 135.0], [0.0, 0.0, 0.0, 0.0]]).T true_lats = np.array([0.0, 0.0, np.nan, 0.0]) @@ -138,8 +151,7 @@ def test_longitudinal_rotation(self): ) def test_rotation_to_equator_sphere(self): - # Rotation to pole at (0.0, 0.0). - # Latitudes become longitudes and vice versa. + """Test rotation to a new pole at latitude=0.0, longitude=0.0.""" true_lats = np.array([-45.0, 0.0, 45.0, 90.0, 45.0, 0.0, -45.0]) true_lons = np.array([90.0, 90.0, 90.0, np.nan, -90.0, -90.0, -90.0]) @@ -156,6 +168,8 @@ def test_rotation_to_equator_sphere(self): class Test__validate_args(ants.tests.TestCase): def test_error_raised(self): + """Test that an error is raised if only --source-cube is passed.""" + args = argparse.Namespace( json_file="json/path", output="output/path", @@ -189,6 +203,10 @@ def setUp(self): def test_transform_coords_not_called( self, mock_lsm, mock_source, mock_transform_coord, mock_validate ): + """Test that when _transform_coordinates is not called when an unrotated pole + is given and only target_lsm is given. The assert statements check that only + appropriate functions are called in this case. Also checks that the unrotated + points are correctly returned.""" mock_lsm.return_value = self.sphere_identity_target_lsm mock_validate.return_value = False @@ -204,6 +222,8 @@ def test_transform_coords_not_called( @mock.patch("ants.cli.ancil_create_shapefile.CubeBuilder") @mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") def test_transform_called(self, mock_lsm, mock_source, mock_transform_coord): + """Test that _transform_coordinates is called when a valid target_lsm is given + and no source_cube is given.""" mock_lsm.return_value = self.sphere_equator_target_lsm @@ -220,6 +240,9 @@ def test_transform_called(self, mock_lsm, mock_source, mock_transform_coord): def test_source_cube_loaded( self, mock_lsm, mock_source, mock_transform_coord, mock_load_cube ): + """Test that _transform_coordinates is called when both target_lsm and + source_cube are given. Check that the source comes from loading from + the cube path.""" mock_lsm.return_value = self.sphere_equator_target_lsm @@ -239,6 +262,8 @@ def setUp(self): @mock.patch("builtins.open", new_callable=mock.mock_open) @mock.patch("ants.cli.ancil_create_shapefile.json.load") def test_json_load(self, mock_json, mock_open, mock_transform): + """Test that the default arguments load a json file, + but do not try to call _transform_if_required.""" mock_json.return_value = self.json_values @@ -251,6 +276,8 @@ def test_json_load(self, mock_json, mock_open, mock_transform): @mock.patch("builtins.open", new_callable=mock.mock_open) @mock.patch("ants.cli.ancil_create_shapefile.json.load") def test_transform_called(self, mock_json, mock_open, mock_transform): + """Test that when target_lsm is given but not source_cube, + _transform_if_required is called once with the correct arguments.""" mock_json.return_value = self.json_values @@ -267,6 +294,8 @@ def test_transform_called(self, mock_json, mock_open, mock_transform): @mock.patch("builtins.open", new_callable=mock.mock_open) @mock.patch("ants.cli.ancil_create_shapefile.json.load") def test_transform_called_with_source(self, mock_json, mock_open, mock_transform): + """Test that when both target_lsm and source_cube are given, + _transform_if_required is called once with the correct arguments.""" mock_json.return_value = self.json_values From 4ad741242ba39442b2cebf79cebc0dd1899468bb Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Tue, 21 Jul 2026 09:26:06 +0100 Subject: [PATCH 13/29] Changed some mocking in tests to improve readability. --- .../tests/cli/test_ancil_create_shapefile.py | 92 ++++++++++--------- 1 file changed, 47 insertions(+), 45 deletions(-) diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index 0163ea6..4260711 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -181,6 +181,9 @@ def test_error_raised(self): _validate_args(args.target_lsm, args.source_cube) +@mock.patch("ants.cli.ancil_create_shapefile._transform_coordinates") +@mock.patch("ants.cli.ancil_create_shapefile.CubeBuilder") +@mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") class Test__transform_if_required(ants.tests.TestCase): def setUp(self): self.points = np.array([[1, 2], [3, 4], [5, 6], [7, 8]]) @@ -197,16 +200,17 @@ def setUp(self): )._cube @mock.patch("ants.cli.ancil_create_shapefile._validate_orientation") - @mock.patch("ants.cli.ancil_create_shapefile._transform_coordinates") - @mock.patch("ants.cli.ancil_create_shapefile.CubeBuilder") - @mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") def test_transform_coords_not_called( - self, mock_lsm, mock_source, mock_transform_coord, mock_validate + self, + mock_validate, + mock_lsm, + mock_source, + mock_transform_coord, ): - """Test that when _transform_coordinates is not called when an unrotated pole - is given and only target_lsm is given. The assert statements check that only - appropriate functions are called in this case. Also checks that the unrotated - points are correctly returned.""" + """Test that _transform_coordinates is not called when an unrotated pole + is given. The assert statements check that only appropriate functions + are called in this case. Also checks that the unrotated points are + correctly returned.""" mock_lsm.return_value = self.sphere_identity_target_lsm mock_validate.return_value = False @@ -218,12 +222,9 @@ def test_transform_coords_not_called( mock_transform_coord.assert_not_called() self.assertArrayEqual(rotated_points, self.points) - @mock.patch("ants.cli.ancil_create_shapefile._transform_coordinates") - @mock.patch("ants.cli.ancil_create_shapefile.CubeBuilder") - @mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") def test_transform_called(self, mock_lsm, mock_source, mock_transform_coord): """Test that _transform_coordinates is called when a valid target_lsm is given - and no source_cube is given.""" + but no source_cube is given.""" mock_lsm.return_value = self.sphere_equator_target_lsm @@ -234,11 +235,12 @@ def test_transform_called(self, mock_lsm, mock_source, mock_transform_coord): mock_transform_coord.assert_called_once() @mock.patch("ants.cli.ancil_create_shapefile.load_cube") - @mock.patch("ants.cli.ancil_create_shapefile._transform_coordinates") - @mock.patch("ants.cli.ancil_create_shapefile.CubeBuilder") - @mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") def test_source_cube_loaded( - self, mock_lsm, mock_source, mock_transform_coord, mock_load_cube + self, + mock_load_cube, + mock_lsm, + mock_source, + mock_transform_coord, ): """Test that _transform_coordinates is called when both target_lsm and source_cube are given. Check that the source comes from loading from @@ -254,56 +256,56 @@ def test_source_cube_loaded( mock_transform_coord.assert_called_once() +@mock.patch("builtins.open", new_callable=mock.mock_open) class Test__load_polygon_from_json(ants.tests.TestCase): def setUp(self): self.json_values = [[1, 2], [3, 4], [5, 6], [7, 8]] - @mock.patch("ants.cli.ancil_create_shapefile._transform_if_required") - @mock.patch("builtins.open", new_callable=mock.mock_open) + def test_json_load(self, *args): + """Test that the default arguments load a json file.""" + + with mock.patch("ants.cli.ancil_create_shapefile.json.load") as mock_json: + _ = _load_polygon_from_json("json/path", None, None) + + mock_json.assert_called_once() + @mock.patch("ants.cli.ancil_create_shapefile.json.load") - def test_json_load(self, mock_json, mock_open, mock_transform): - """Test that the default arguments load a json file, - but do not try to call _transform_if_required.""" + def test_transform_not_called(self, mock_json, *args): + """Test that transformation is not performed if target_lsm + is not passed.""" mock_json.return_value = self.json_values - _ = _load_polygon_from_json("json/path", None, None) + with mock.patch( + "ants.cli.ancil_create_shapefile._transform_if_required" + ) as mock_transform: + _ = _load_polygon_from_json("json/path", None, None) - mock_json.assert_called_once() mock_transform.assert_not_called() - @mock.patch("ants.cli.ancil_create_shapefile._transform_if_required") - @mock.patch("builtins.open", new_callable=mock.mock_open) @mock.patch("ants.cli.ancil_create_shapefile.json.load") - def test_transform_called(self, mock_json, mock_open, mock_transform): - """Test that when target_lsm is given but not source_cube, - _transform_if_required is called once with the correct arguments.""" + def test_transform_called(self, mock_json, *args): + """Test that transformation is called if target_lsm_path is passed.""" mock_json.return_value = self.json_values - _ = _load_polygon_from_json("json/path", "lsm/path", None) + with mock.patch( + "ants.cli.ancil_create_shapefile._transform_if_required" + ) as mock_transform: + _ = _load_polygon_from_json("json/path", "lsm/path", None) - target_path, source_path, points = mock_transform.call_args.args - mock_json.assert_called_once() mock_transform.assert_called_once() - self.assertArrayEqual(points, self.json_values) - self.assertEqual("lsm/path", target_path) - self.assertEqual(None, source_path) - @mock.patch("ants.cli.ancil_create_shapefile._transform_if_required") - @mock.patch("builtins.open", new_callable=mock.mock_open) @mock.patch("ants.cli.ancil_create_shapefile.json.load") - def test_transform_called_with_source(self, mock_json, mock_open, mock_transform): - """Test that when both target_lsm and source_cube are given, - _transform_if_required is called once with the correct arguments.""" + def test_transform_called_with_source(self, mock_json, *args): + """Test that transformation is called if target_lsm_path and source_cube_path + are passed.""" mock_json.return_value = self.json_values - _ = _load_polygon_from_json("json/path", "lsm/path", "source/path") + with mock.patch( + "ants.cli.ancil_create_shapefile._transform_if_required" + ) as mock_transform: + _ = _load_polygon_from_json("json/path", "lsm/path", "source/path") - target_path, source_path, points = mock_transform.call_args.args - mock_json.assert_called_once() mock_transform.assert_called_once() - self.assertArrayEqual(points, self.json_values) - self.assertEqual("lsm/path", target_path) - self.assertEqual("source/path", source_path) From 5ea3e014f87fd0870d878d5256b3260b1d84e014 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Tue, 21 Jul 2026 12:32:47 +0100 Subject: [PATCH 14/29] Corrected spelling mistakes. --- lib/ants/cli/ancil_create_shapefile.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index c903bee..f0e63f8 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -161,7 +161,7 @@ def _transform_if_required(target_lsm_path, source_cube_path, points): Parameters ---------- target_lsm_path : str - File path for a land sea mask that provides the new rotated pole + File path to a land sea mask that provides the new rotated pole coordinates to which the longitude, latitude pairs will be transformed. source_cube_path : str File path to an iris cube specifying the coordinate system of the @@ -206,7 +206,7 @@ def _load_polygon_from_json(json_file, target_lsm_path, source_cube_path): json_file : str Path to json file target_lsm_path : str - File path for a land sea mask that provides the new rotated pole + File path to a land sea mask that provides the new rotated pole coordinates to which the longitude, latitude pairs will be transformed. source_cube_path : str File path to an iris cube specifying the coordinate system of the @@ -246,10 +246,10 @@ def main(json_file, output, target_lsm_path, source_cube_path): output : str Location to store generated shape file target_lsm_path : str - File path for a land sea mask that provides the new rotated pole + File path to a land sea mask that provides the new rotated pole coordinates to which the longitude, latitude pairs will be transformed. source_cube_path : str - File path to an iris cube specifying the cosordinate system of the + File path to an iris cube specifying the coordinate system of the input json file. """ From c848fb5b9b2baad56ba75313af0c38f2a3882833 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Thu, 23 Jul 2026 12:54:49 +0100 Subject: [PATCH 15/29] Moved loading out of function to make testing easier. --- lib/ants/cli/ancil_create_shapefile.py | 96 +++++++++----- .../tests/cli/test_ancil_create_shapefile.py | 125 +++++++++--------- 2 files changed, 124 insertions(+), 97 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index f0e63f8..9c4302d 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -61,8 +61,8 @@ def _check_coord_system_type(target_lsm): def _validate_orientation(target_lsm): """ - Check if the rotated pole coordinate system specified by the target lsm - is coincident with the north pole. + Validate that the rotated pole coordinate system is not coincident with the + geographic north pole. Parameters ---------- @@ -71,14 +71,13 @@ def _validate_orientation(target_lsm): Returns ------- - valid_crs : bool - True if target_lsm is a rotated pole, - False otherwise. + bool + True if target_lsm is a valid rotated pole, False otherwise. Warns - ------ + ----- UserWarning - If target_lsm is already at latitude=90.0 and longitude=0.0. + If the rotated pole is coincident with the geographic north pole. """ target_crs = target_lsm.coord_system() @@ -110,26 +109,25 @@ def _validate_args(target_lsm_path, source_cube_path): def _transform_coordinates(target_lsm, source_cube, points): """ - Transforms the longitude, latitude points in the source coordinate - system to the rotated pole coordinate system defined by target_lsm. + Transform the longitude-latitude points to the rotated pole. - The source coordinate system is assumed to be unrotated geodetic - defined on a sphere, unless otherwise specified. + The source coordinate system is assumed to be unrotated geodetic defined + on a sphere, unless otherwise specified. Parameters ---------- target_lsm : :class:`iris.cube.Cube` - The lsm cube specifiying the rotated pole coordinate system. + The lsm cube specifying the rotated pole coordinate system. source_cube : :class:`iris.cube.Cube` - An iris cube specifying the coordinate system of the - input json file. - points : :class:`np.ndarry` - An (m,2) shaped numpy array of m longitude and m latitude points. + An iris cube specifying the coordinate system of the input json + file. + points : :class:`numpy.ndarray` + An ``(m, 2)`` numpy array of m longitude-latitude pairs. Returns ------- - rotated_points : :class:`np.ndarry` - An (m,2) shaped numpy array of transformed longitude and latitude points. + : :class:`numpy.ndarray` + An ``(m, 2)`` numpy array of transformed longitude-latitude pairs. """ source_crs = source_cube.coord_system().as_cartopy_crs() @@ -151,32 +149,25 @@ def _transform_coordinates(target_lsm, source_cube, points): return rotated_points -def _transform_if_required(target_lsm_path, source_cube_path, points): +def _load_cubes(target_lsm_path, source_cube_path): """ - Load the target lsm and create a source cube if not provided. The points - are then transformed to the specified rotated pole coordinate system - if the provided coordinate system is a rotated pole and is not - coincident with the north pole. + Load the target lsm and create a source cube if not provided. Parameters ---------- target_lsm_path : str - File path to a land sea mask that provides the new rotated pole - coordinates to which the longitude, latitude pairs will be transformed. + File path to a land sea mask that provides the new rotated pole. source_cube_path : str File path to an iris cube specifying the coordinate system of the input json file. - points : :class:`np.ndarry` - An (m,2) shaped numpy array of longitude and latitude points. Returns ------- - rotated_points : :class:`np.ndarry` - An (m,2) shaped numpy array of transformed longitude and latitude points. + : tuple(:class:`iris.cube.Cube`, :class:`iris.cube.Cube`) + A tuple containing the target lsm and the source cube respectively. """ target_lsm = load_landsea_mask(target_lsm_path) - _check_coord_system_type(target_lsm) if source_cube_path is None: crs = iris.coord_systems.GeogCS(6371229.0) @@ -184,6 +175,38 @@ def _transform_if_required(target_lsm_path, source_cube_path, points): else: source_cube = load_cube(source_cube_path) + return target_lsm, source_cube + + +def _transform_if_required(target_lsm_path, source_cube_path, points): + """ + Perform the transformation to a rotated pole if the coordinate system + is valid. + + The points are transformed to the specified rotated pole coordinate system + if the provided coordinate system is a rotated pole and is not coincident + with the geographic north pole. + + Parameters + ---------- + target_lsm_path : str + File path to a land sea mask that provides the new rotated pole. + source_cube_path : str + File path to an iris cube specifying the coordinate system of the + input json file. + points : :class:`numpy.ndarray` + An ``(m, 2)`` numpy array of m longitude-latitude pairs. + + Returns + ------- + : :class:`numpy.ndarray` + An ``(m, 2)`` numpy array of transformed longitude-latitude pairs. + """ + + target_lsm, source_cube = _load_cubes(target_lsm_path, source_cube_path) + + _check_coord_system_type(target_lsm) + if _validate_orientation(target_lsm): rotated_points = _transform_coordinates(target_lsm, source_cube, points) else: @@ -194,7 +217,7 @@ def _transform_if_required(target_lsm_path, source_cube_path, points): def _load_polygon_from_json(json_file, target_lsm_path, source_cube_path): """ - Load a json file containing a list of pairs of longitude, latitude points + Load a json file containing a list of pairs of longitude-latitude points to create a polygon from. If a valid rotated pole is specifed in the target lsm, first transforms @@ -206,8 +229,7 @@ def _load_polygon_from_json(json_file, target_lsm_path, source_cube_path): json_file : str Path to json file target_lsm_path : str - File path to a land sea mask that provides the new rotated pole - coordinates to which the longitude, latitude pairs will be transformed. + File path to a land sea mask that provides the new rotated pole. source_cube_path : str File path to an iris cube specifying the coordinate system of the input json file. @@ -216,6 +238,7 @@ def _load_polygon_from_json(json_file, target_lsm_path, source_cube_path): ------- : :class:`~shapely.geometry.Polygon` """ + with open(json_file, "r") as polygon_json: points = json.load(polygon_json) points = np.array(points) @@ -230,12 +253,14 @@ def _load_polygon_from_json(json_file, target_lsm_path, source_cube_path): def main(json_file, output, target_lsm_path, source_cube_path): """ + Create a shape file from pairs of longitude, latitude points. + Loads in a provided json file that defines pairs of longitude, latitude points to create a polygon from. That polygon is then used to create a shape file that is saved to the specified output location. If target_lsm_path is provided, the points are first transformed from an - unrotated geodetic coordinate system to a rotated pole coordinate system + source geodetic coordinate system to a rotated pole coordinate system specified by the lsm. It is assumed that the points in the json file are on an unrotated geodetic grid, unless otherwise specified. @@ -246,8 +271,7 @@ def main(json_file, output, target_lsm_path, source_cube_path): output : str Location to store generated shape file target_lsm_path : str - File path to a land sea mask that provides the new rotated pole - coordinates to which the longitude, latitude pairs will be transformed. + File path to a land sea mask that provides the new rotated pole. source_cube_path : str File path to an iris cube specifying the coordinate system of the input json file. diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index 4260711..4d4483f 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -12,6 +12,7 @@ import numpy as np from ants.cli.ancil_create_shapefile import ( _check_coord_system_type, + _load_cubes, _load_polygon_from_json, _transform_coordinates, _transform_if_required, @@ -93,7 +94,7 @@ def setUp(self): self.sphere_identity_target_lsm = CubeBuilder( self.sphere_identity_crs, (2, 2) )._cube - self.sphere_source = geodetic((2, 2), crs=self.sphere_crs) + self.sphere_source = CubeBuilder(self.sphere_crs, (2, 2))._cube self.sphere_equator_target_lsm = geodetic( (2, 2), north_pole_lat=0.0, north_pole_lon=0.0, crs=self.sphere_crs ) @@ -111,7 +112,7 @@ def setUp(self): ).T def test_identity_rotation_sphere(self): - """Test that a rotation to a pole at latitude=90.0, longitude=0.0 + """Test that rotation to a pole at latitude=90.0, longitude=0.0 returns the same points, except a 180.0 degree rotation in longitude. By convention the new pole is defined 180.0 rotated from the provided longitude.""" @@ -132,7 +133,7 @@ def test_identity_rotation_sphere(self): def test_longitudinal_rotation(self): """Test rotation to a new pole at latitude=90.0, longitude=90.0. - The new pole is located at latitude=90.0, longitude 270.0 in the rotated + The new pole is located at latitude=90.0, longitude=270.0 in the rotated pole coordinate system.""" points = np.array([[0.0, 45.0, 90.0, 135.0], [0.0, 0.0, 0.0, 0.0]]).T @@ -181,14 +182,44 @@ def test_error_raised(self): _validate_args(args.target_lsm, args.source_cube) -@mock.patch("ants.cli.ancil_create_shapefile._transform_coordinates") -@mock.patch("ants.cli.ancil_create_shapefile.CubeBuilder") -@mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") +class Test__load_cubes(ants.tests.TestCase): + def test_load_lsm_called(self): + """Test that loading the landsea mask is successfully called.""" + + with mock.patch( + "ants.cli.ancil_create_shapefile.load_landsea_mask" + ) as mock_lsm: + _ = _load_cubes("target/path", None) + + mock_lsm.assert_called_once_with("target/path") + + @mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") + def test_cubebuilder_called(self, *args): + """Test source cube created if no path given.""" + + with mock.patch( + "ants.cli.ancil_create_shapefile.CubeBuilder" + ) as mock_build_cube: + _ = _load_cubes("target/path", None) + + mock_build_cube.assert_called_once() + + @mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") + def test_load_cube_called(self, *args): + """Test load cube called if source path given.""" + + with mock.patch("ants.cli.ancil_create_shapefile.load_cube") as mock_load_cube: + _ = _load_cubes("target/path", "source/cube") + + mock_load_cube.assert_called_once() + + +@mock.patch("ants.cli.ancil_create_shapefile._load_cubes") class Test__transform_if_required(ants.tests.TestCase): def setUp(self): self.points = np.array([[1, 2], [3, 4], [5, 6], [7, 8]]) self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) - self.sphere_source = geodetic((2, 2), crs=self.sphere_crs) + self.sphere_source = CubeBuilder(self.sphere_crs, (2, 2))._cube self.sphere_equator_target_lsm = geodetic( (2, 2), north_pole_lat=0.0, north_pole_lon=0.0, crs=self.sphere_crs ) @@ -199,61 +230,33 @@ def setUp(self): self.sphere_identity_crs, (2, 2) )._cube - @mock.patch("ants.cli.ancil_create_shapefile._validate_orientation") - def test_transform_coords_not_called( - self, - mock_validate, - mock_lsm, - mock_source, - mock_transform_coord, - ): - """Test that _transform_coordinates is not called when an unrotated pole - is given. The assert statements check that only appropriate functions - are called in this case. Also checks that the unrotated points are - correctly returned.""" - - mock_lsm.return_value = self.sphere_identity_target_lsm - mock_validate.return_value = False - - rotated_points = _transform_if_required("lsm/path", None, self.points) - - mock_lsm.assert_called_once() - mock_source.assert_called_once() - mock_transform_coord.assert_not_called() - self.assertArrayEqual(rotated_points, self.points) - - def test_transform_called(self, mock_lsm, mock_source, mock_transform_coord): - """Test that _transform_coordinates is called when a valid target_lsm is given - but no source_cube is given.""" + def test_non_rotated_returned(self, mock_load): + """Test that the input points are returned if a non-rotated pole is + provided as the target_lsm.""" - mock_lsm.return_value = self.sphere_equator_target_lsm + mock_load.return_value = self.sphere_identity_target_lsm, self.sphere_source - _ = _transform_if_required("lsm/path", None, self.points) + with self.assertWarns(UserWarning): + rotated_points = _transform_if_required("target/lsm", None, self.points) - mock_lsm.assert_called_once() - mock_source.assert_called_once() - mock_transform_coord.assert_called_once() + self.assertArrayEqual(rotated_points, self.points) - @mock.patch("ants.cli.ancil_create_shapefile.load_cube") - def test_source_cube_loaded( - self, - mock_load_cube, - mock_lsm, - mock_source, - mock_transform_coord, - ): - """Test that _transform_coordinates is called when both target_lsm and - source_cube are given. Check that the source comes from loading from - the cube path.""" + def test_transform_called(self, mock_load): + """Test that transform coordinates is called with the correct arguments + when a valid target_lsm is given.""" - mock_lsm.return_value = self.sphere_equator_target_lsm + mock_load.return_value = self.sphere_equator_target_lsm, self.sphere_source - _ = _transform_if_required("lsm/path", "source/path", self.points) + with mock.patch( + "ants.cli.ancil_create_shapefile._transform_coordinates" + ) as mock_transform: + _ = _transform_if_required("target/lsm", None, self.points) - mock_lsm.assert_called_once() - mock_source.assert_not_called() - mock_load_cube.assert_called_once() - mock_transform_coord.assert_called_once() + target_lsm, source_cube, points = mock_transform.call_args.args + mock_transform.assert_called_once() + self.assertEqual(target_lsm, self.sphere_equator_target_lsm) + self.assertEqual(source_cube, self.sphere_source) + self.assertArrayEqual(points, self.points) @mock.patch("builtins.open", new_callable=mock.mock_open) @@ -278,10 +281,10 @@ def test_transform_not_called(self, mock_json, *args): with mock.patch( "ants.cli.ancil_create_shapefile._transform_if_required" - ) as mock_transform: + ) as mock_transform_if: _ = _load_polygon_from_json("json/path", None, None) - mock_transform.assert_not_called() + mock_transform_if.assert_not_called() @mock.patch("ants.cli.ancil_create_shapefile.json.load") def test_transform_called(self, mock_json, *args): @@ -291,10 +294,10 @@ def test_transform_called(self, mock_json, *args): with mock.patch( "ants.cli.ancil_create_shapefile._transform_if_required" - ) as mock_transform: + ) as mock_transform_if: _ = _load_polygon_from_json("json/path", "lsm/path", None) - mock_transform.assert_called_once() + mock_transform_if.assert_called_once() @mock.patch("ants.cli.ancil_create_shapefile.json.load") def test_transform_called_with_source(self, mock_json, *args): @@ -305,7 +308,7 @@ def test_transform_called_with_source(self, mock_json, *args): with mock.patch( "ants.cli.ancil_create_shapefile._transform_if_required" - ) as mock_transform: + ) as mock_transform_if: _ = _load_polygon_from_json("json/path", "lsm/path", "source/path") - mock_transform.assert_called_once() + mock_transform_if.assert_called_once() From 5f07ac18f56b6e89daf955d15f80421d2d989ae2 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Mon, 10 Aug 2026 10:26:38 +0100 Subject: [PATCH 16/29] Moved transformation out of polygon loading. --- lib/ants/cli/ancil_create_shapefile.py | 48 +++++--------- .../tests/cli/test_ancil_create_shapefile.py | 66 ++++--------------- 2 files changed, 30 insertions(+), 84 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index 9c4302d..10c4f63 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -178,7 +178,7 @@ def _load_cubes(target_lsm_path, source_cube_path): return target_lsm, source_cube -def _transform_if_required(target_lsm_path, source_cube_path, points): +def _transform_if_required(target_lsm, source_cube, points): """ Perform the transformation to a rotated pole if the coordinate system is valid. @@ -189,11 +189,11 @@ def _transform_if_required(target_lsm_path, source_cube_path, points): Parameters ---------- - target_lsm_path : str - File path to a land sea mask that provides the new rotated pole. - source_cube_path : str - File path to an iris cube specifying the coordinate system of the - input json file. + target_lsm : :class:`iris.cube.Cube` + The lsm cube specifying the rotated pole coordinate system. + source_cube : :class:`iris.cube.Cube` + An iris cube specifying the coordinate system of the input json + file. points : :class:`numpy.ndarray` An ``(m, 2)`` numpy array of m longitude-latitude pairs. @@ -203,8 +203,6 @@ def _transform_if_required(target_lsm_path, source_cube_path, points): An ``(m, 2)`` numpy array of transformed longitude-latitude pairs. """ - target_lsm, source_cube = _load_cubes(target_lsm_path, source_cube_path) - _check_coord_system_type(target_lsm) if _validate_orientation(target_lsm): @@ -215,40 +213,25 @@ def _transform_if_required(target_lsm_path, source_cube_path, points): return rotated_points -def _load_polygon_from_json(json_file, target_lsm_path, source_cube_path): +def _load_points_from_json(json_file): """ - Load a json file containing a list of pairs of longitude-latitude points - to create a polygon from. - - If a valid rotated pole is specifed in the target lsm, first transforms - the points from the source coordinate system to the rotated pole before - creating a polygon. + Load a json file containing a list of pairs of longitude-latitude points. Parameters ---------- json_file : str - Path to json file - target_lsm_path : str - File path to a land sea mask that provides the new rotated pole. - source_cube_path : str - File path to an iris cube specifying the coordinate system of the - input json file. + Path to json file. Returns ------- - : :class:`~shapely.geometry.Polygon` + : :class:`numpy.ndarray` """ with open(json_file, "r") as polygon_json: points = json.load(polygon_json) points = np.array(points) - if target_lsm_path is not None: - points = _transform_if_required(target_lsm_path, source_cube_path, points) - - polygon = Polygon(points) - - return polygon + return points def main(json_file, output, target_lsm_path, source_cube_path): @@ -259,7 +242,7 @@ def main(json_file, output, target_lsm_path, source_cube_path): points to create a polygon from. That polygon is then used to create a shape file that is saved to the specified output location. - If target_lsm_path is provided, the points are first transformed from an + If target_lsm_path is provided, the points are first transformed from a source geodetic coordinate system to a rotated pole coordinate system specified by the lsm. It is assumed that the points in the json file are on an unrotated geodetic grid, unless otherwise specified. @@ -278,8 +261,13 @@ def main(json_file, output, target_lsm_path, source_cube_path): """ # Load a json and make a polygon - polygon = _load_polygon_from_json(json_file, target_lsm_path, source_cube_path) + points = _load_points_from_json(json_file) + if target_lsm_path is not None: + target_lsm, source_cube = _load_cubes(target_lsm_path, source_cube_path) + points = _transform_if_required(target_lsm, source_cube, points) + + polygon = Polygon(points) # Now convert it to a shapefile with OGR driver = ogr.GetDriverByName("Esri Shapefile") datasource = driver.CreateDataSource(output) diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index 4d4483f..1115792 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -13,7 +13,7 @@ from ants.cli.ancil_create_shapefile import ( _check_coord_system_type, _load_cubes, - _load_polygon_from_json, + _load_points_from_json, _transform_coordinates, _transform_if_required, _validate_args, @@ -214,7 +214,6 @@ def test_load_cube_called(self, *args): mock_load_cube.assert_called_once() -@mock.patch("ants.cli.ancil_create_shapefile._load_cubes") class Test__transform_if_required(ants.tests.TestCase): def setUp(self): self.points = np.array([[1, 2], [3, 4], [5, 6], [7, 8]]) @@ -230,27 +229,27 @@ def setUp(self): self.sphere_identity_crs, (2, 2) )._cube - def test_non_rotated_returned(self, mock_load): + def test_non_rotated_returned(self): """Test that the input points are returned if a non-rotated pole is provided as the target_lsm.""" - mock_load.return_value = self.sphere_identity_target_lsm, self.sphere_source - with self.assertWarns(UserWarning): - rotated_points = _transform_if_required("target/lsm", None, self.points) + rotated_points = _transform_if_required( + self.sphere_identity_target_lsm, self.sphere_source, self.points + ) self.assertArrayEqual(rotated_points, self.points) - def test_transform_called(self, mock_load): + def test_transform_called(self): """Test that transform coordinates is called with the correct arguments when a valid target_lsm is given.""" - mock_load.return_value = self.sphere_equator_target_lsm, self.sphere_source - with mock.patch( "ants.cli.ancil_create_shapefile._transform_coordinates" ) as mock_transform: - _ = _transform_if_required("target/lsm", None, self.points) + _ = _transform_if_required( + self.sphere_equator_target_lsm, self.sphere_source, self.points + ) target_lsm, source_cube, points = mock_transform.call_args.args mock_transform.assert_called_once() @@ -260,55 +259,14 @@ def test_transform_called(self, mock_load): @mock.patch("builtins.open", new_callable=mock.mock_open) -class Test__load_polygon_from_json(ants.tests.TestCase): +class Test__load_points_from_json(ants.tests.TestCase): def setUp(self): self.json_values = [[1, 2], [3, 4], [5, 6], [7, 8]] def test_json_load(self, *args): - """Test that the default arguments load a json file.""" + """Test that a json file is loaded.""" with mock.patch("ants.cli.ancil_create_shapefile.json.load") as mock_json: - _ = _load_polygon_from_json("json/path", None, None) + _ = _load_points_from_json("json/path") mock_json.assert_called_once() - - @mock.patch("ants.cli.ancil_create_shapefile.json.load") - def test_transform_not_called(self, mock_json, *args): - """Test that transformation is not performed if target_lsm - is not passed.""" - - mock_json.return_value = self.json_values - - with mock.patch( - "ants.cli.ancil_create_shapefile._transform_if_required" - ) as mock_transform_if: - _ = _load_polygon_from_json("json/path", None, None) - - mock_transform_if.assert_not_called() - - @mock.patch("ants.cli.ancil_create_shapefile.json.load") - def test_transform_called(self, mock_json, *args): - """Test that transformation is called if target_lsm_path is passed.""" - - mock_json.return_value = self.json_values - - with mock.patch( - "ants.cli.ancil_create_shapefile._transform_if_required" - ) as mock_transform_if: - _ = _load_polygon_from_json("json/path", "lsm/path", None) - - mock_transform_if.assert_called_once() - - @mock.patch("ants.cli.ancil_create_shapefile.json.load") - def test_transform_called_with_source(self, mock_json, *args): - """Test that transformation is called if target_lsm_path and source_cube_path - are passed.""" - - mock_json.return_value = self.json_values - - with mock.patch( - "ants.cli.ancil_create_shapefile._transform_if_required" - ) as mock_transform_if: - _ = _load_polygon_from_json("json/path", "lsm/path", "source/path") - - mock_transform_if.assert_called_once() From 52dd4283f7d0767e3619f8d437d4fb6f6141833e Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Mon, 10 Aug 2026 16:15:54 +0100 Subject: [PATCH 17/29] Add code to address antimeridian crossing. --- lib/ants/cli/ancil_create_shapefile.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index 10c4f63..86a249c 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -111,8 +111,9 @@ def _transform_coordinates(target_lsm, source_cube, points): """ Transform the longitude-latitude points to the rotated pole. - The source coordinate system is assumed to be unrotated geodetic defined - on a sphere, unless otherwise specified. + If the distance between transformed longitude points exceeds 180.0, + it is assumed these points cross the antimeridian. The points are + changed so they are instead defined in the interval [0.0, 360.0]. Parameters ---------- @@ -138,8 +139,15 @@ def _transform_coordinates(target_lsm, source_cube, points): source_crs, points[:, 0], points[:, 1] )[:, :2] + closed_longitudes = np.vstack([points, points[0, :]]) + long_diff = np.abs(closed_longitudes[:-1, 0] - closed_longitudes[1:, 0]) + + if np.any(long_diff > 180.0): + neg_indices = np.where(points[:, 0] < 0) + points[neg_indices, 0] += 360.0 + _LOGGER.info( - "Input json file transformed to new pole rotated coordinate system at" + "Input json file transformed to new pole rotated coordinate system at " "pole longitude=%s, pole latitude=%s, central rotated longitude=%s.", target_coord.grid_north_pole_longitude, target_coord.grid_north_pole_latitude, From 8f9799cebf5012c35a8cca11156acdc4f335db90 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Tue, 11 Aug 2026 09:59:29 +0100 Subject: [PATCH 18/29] Fixed using wrong points for antimeridian crossing. --- lib/ants/cli/ancil_create_shapefile.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index 86a249c..c1713a7 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -139,12 +139,12 @@ def _transform_coordinates(target_lsm, source_cube, points): source_crs, points[:, 0], points[:, 1] )[:, :2] - closed_longitudes = np.vstack([points, points[0, :]]) + closed_longitudes = np.vstack([rotated_points, rotated_points[0, :]]) long_diff = np.abs(closed_longitudes[:-1, 0] - closed_longitudes[1:, 0]) if np.any(long_diff > 180.0): - neg_indices = np.where(points[:, 0] < 0) - points[neg_indices, 0] += 360.0 + neg_indices = np.where(rotated_points[:, 0] < 0) + rotated_points[neg_indices, 0] += 360.0 _LOGGER.info( "Input json file transformed to new pole rotated coordinate system at " From 2cc209b489017dc88953a9be539c1cd57bf89e95 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Tue, 11 Aug 2026 10:07:09 +0100 Subject: [PATCH 19/29] Changed variable names to correct long misname. --- lib/ants/cli/ancil_create_shapefile.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index c1713a7..d26738b 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -139,10 +139,10 @@ def _transform_coordinates(target_lsm, source_cube, points): source_crs, points[:, 0], points[:, 1] )[:, :2] - closed_longitudes = np.vstack([rotated_points, rotated_points[0, :]]) - long_diff = np.abs(closed_longitudes[:-1, 0] - closed_longitudes[1:, 0]) + closed_lon = np.vstack([rotated_points, rotated_points[0, :]]) + lon_diff = np.abs(closed_lon[:-1, 0] - closed_lon[1:, 0]) - if np.any(long_diff > 180.0): + if np.any(lon_diff > 180.0): neg_indices = np.where(rotated_points[:, 0] < 0) rotated_points[neg_indices, 0] += 360.0 From 2534d8b109cc5c8545ab25f59f44d12a561f8bbb Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Tue, 11 Aug 2026 10:17:07 +0100 Subject: [PATCH 20/29] Removed double negative. --- lib/ants/cli/ancil_create_shapefile.py | 11 +++++------ lib/ants/tests/cli/test_ancil_create_shapefile.py | 2 +- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index d26738b..8a41642 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -86,18 +86,17 @@ def _validate_orientation(target_lsm): grid_lat = target_crs.grid_north_pole_latitude pole_lon_rotation = target_crs.north_pole_grid_longitude - valid_crs = not ( - ants.utils.ndarray.allclose( - [grid_lon, grid_lat, pole_lon_rotation], [0.0, 90.0, 0.0] - ) + invalid_crs = ants.utils.ndarray.allclose( + [grid_lon, grid_lat, pole_lon_rotation], [0.0, 90.0, 0.0] ) - if not valid_crs: + + if invalid_crs: warnings.warn( "target_lsm has a geodetic coordinate system with pole located" f"at grid_longitude={grid_lon}, grid_latitude={grid_lat}." "No transformation will be carried out." ) - return valid_crs + return not invalid_crs def _validate_args(target_lsm_path, source_cube_path): diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index 1115792..59d14ba 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -138,7 +138,7 @@ def test_longitudinal_rotation(self): points = np.array([[0.0, 45.0, 90.0, 135.0], [0.0, 0.0, 0.0, 0.0]]).T true_lats = np.array([0.0, 0.0, np.nan, 0.0]) - true_lons = np.array([90.0, 135.0, -180.0, -135.0]) + true_lons = np.array([90.0, 135.0, 180.0, 225.0]) expected_rotation = np.array([true_lons, true_lats]).T check_mask = ~np.isnan(expected_rotation) From 362d65384d7eea47657349caae1a5806c1d7bb19 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Thu, 13 Aug 2026 16:18:56 +0100 Subject: [PATCH 21/29] Added extra checks and tests for polygon validity and orientation. --- lib/ants/cli/ancil_create_shapefile.py | 39 +++++++++++++++++-- .../tests/cli/test_ancil_create_shapefile.py | 36 ++++++++++++++++- 2 files changed, 70 insertions(+), 5 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index 8a41642..caa8422 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -28,6 +28,7 @@ from ants.utils.cube import CubeBuilder from osgeo import ogr from shapely.geometry import Polygon +from shapely.validation import explain_validity _LOGGER = logging.getLogger(__name__) @@ -59,6 +60,34 @@ def _check_coord_system_type(target_lsm): ) +def _check_polygon_validity(polygon, ccw_expected): + """ + Check if polygon is valid and has the expected orientation. + + Parameters + ---------- + polygon : :class:`~shapely.geometry.Polygon` + The polygon made from the input json file. + ccw_expected : bool + Expected polygon orientation before rotation. + + Raises + ------ + ValueError : + If the polygon is invalid (for example intersecting edges). + ValueError : + If the polygon has a different orientation after rotation. + """ + + if not polygon.is_valid: + raise ValueError(f"Polygon is invalid: {explain_validity(polygon)}") + if polygon.exterior.is_ccw != ccw_expected: + raise ValueError( + f"Polygon orientation has changed. Expected is_ccw={ccw_expected}, " + f"current polygon has is_ccw={polygon.exterior.is_ccw}." + ) + + def _validate_orientation(target_lsm): """ Validate that the rotated pole coordinate system is not coincident with the @@ -93,8 +122,8 @@ def _validate_orientation(target_lsm): if invalid_crs: warnings.warn( "target_lsm has a geodetic coordinate system with pole located" - f"at grid_longitude={grid_lon}, grid_latitude={grid_lat}." - "No transformation will be carried out." + f" at grid_longitude={grid_lon}, grid_latitude={grid_lat}." + " No transformation will be carried out." ) return not invalid_crs @@ -269,12 +298,16 @@ def main(json_file, output, target_lsm_path, source_cube_path): # Load a json and make a polygon points = _load_points_from_json(json_file) + ccw_expected = Polygon(points).exterior.is_ccw if target_lsm_path is not None: target_lsm, source_cube = _load_cubes(target_lsm_path, source_cube_path) points = _transform_if_required(target_lsm, source_cube, points) polygon = Polygon(points) + + _check_polygon_validity(polygon, ccw_expected) + # Now convert it to a shapefile with OGR driver = ogr.GetDriverByName("Esri Shapefile") datasource = driver.CreateDataSource(output) @@ -316,7 +349,7 @@ def _get_parser(): "--source-cube", type=ants.config.filepath_readable, required=False, - help="Path to an iris cube which specifies the co-ordinate" + help="Path to an iris cube which specifies the coordinate" " system of the json file.", ) return parser diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index 59d14ba..0d66ccb 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -12,6 +12,7 @@ import numpy as np from ants.cli.ancil_create_shapefile import ( _check_coord_system_type, + _check_polygon_validity, _load_cubes, _load_points_from_json, _transform_coordinates, @@ -21,6 +22,7 @@ ) from ants.tests.stock import geodetic from ants.utils.cube import CubeBuilder +from shapely.geometry import Polygon class Test__check_coord_system_type(ants.tests.TestCase): @@ -41,6 +43,36 @@ def test_unrotated_target_lsm(self): _check_coord_system_type(target_lsm) +class Test__check_polygon_validity(ants.tests.TestCase): + def test_invalid_polygon(self): + """Test that an invalid polygon with intersecting boundaries + raises an error.""" + + self.invalid_points = np.array([[1, 1], [-1, -1], [1, -1], [-1, 1]]) + self.invalid_poly = Polygon(self.invalid_points) + self.is_ccw = False + + error_msg = "Polygon is invalid: " + + with self.assertRaisesRegex(ValueError, error_msg): + _check_polygon_validity(self.invalid_poly, self.is_ccw) + + def test_changed_orientation(self): + """Test that a change in polygon orientation raises an error.""" + + self.ccw_points = np.array([[1, 1], [-1, 1], [-1, -1], [1, -1]]) + self.cw_poly = Polygon(self.ccw_points) + self.is_ccw = False + + error_msg = ( + "Polygon orientation has changed. Expected is_ccw=False, " + "current polygon has is_ccw=True." + ) + + with self.assertRaisesRegex(ValueError, error_msg): + _check_polygon_validity(self.cw_poly, self.is_ccw) + + class Test__validate_orientation(ants.tests.TestCase): def setUp(self): self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) @@ -61,8 +93,8 @@ def setUp(self): ) self.warning_msg = ( "target_lsm has a geodetic coordinate system with pole located" - "at grid_longitude=0.0, grid_latitude=90.0." - "No transformation will be carried out." + " at grid_longitude=0.0, grid_latitude=90.0." + " No transformation will be carried out." ) def test_warning_raised(self): From 368b878c2965d4ad4bbfca59ce7d5d13388888cc Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Tue, 18 Aug 2026 10:54:05 +0100 Subject: [PATCH 22/29] Renamed polygon variable name in test to match orientation. --- lib/ants/tests/cli/test_ancil_create_shapefile.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index 0d66ccb..4b9a8f1 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -61,7 +61,7 @@ def test_changed_orientation(self): """Test that a change in polygon orientation raises an error.""" self.ccw_points = np.array([[1, 1], [-1, 1], [-1, -1], [1, -1]]) - self.cw_poly = Polygon(self.ccw_points) + self.ccw_poly = Polygon(self.ccw_points) self.is_ccw = False error_msg = ( @@ -70,7 +70,7 @@ def test_changed_orientation(self): ) with self.assertRaisesRegex(ValueError, error_msg): - _check_polygon_validity(self.cw_poly, self.is_ccw) + _check_polygon_validity(self.ccw_poly, self.is_ccw) class Test__validate_orientation(ants.tests.TestCase): From e1d61ea1272782449100dd9e76627a8c8a5d4a1e Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Thu, 20 Aug 2026 13:38:53 +0100 Subject: [PATCH 23/29] Addressed review on readability and docstrings. --- lib/ants/cli/ancil_create_shapefile.py | 79 ++++++------ .../tests/cli/test_ancil_create_shapefile.py | 113 +++++++++--------- 2 files changed, 98 insertions(+), 94 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index caa8422..3f9b8dc 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -69,14 +69,14 @@ def _check_polygon_validity(polygon, ccw_expected): polygon : :class:`~shapely.geometry.Polygon` The polygon made from the input json file. ccw_expected : bool - Expected polygon orientation before rotation. + Expected polygon orientation before transformation. Raises ------ ValueError : If the polygon is invalid (for example intersecting edges). ValueError : - If the polygon has a different orientation after rotation. + If the polygon has a different orientation after transformation. """ if not polygon.is_valid: @@ -128,26 +128,29 @@ def _validate_orientation(target_lsm): return not invalid_crs -def _validate_args(target_lsm_path, source_cube_path): - if source_cube_path is not None and target_lsm_path is None: - raise ValueError( - "If --source-cube is passed then --target-lsm must also be given." - ) +def _validate_args(target_lsm_path, source_path): + if source_path is not None and target_lsm_path is None: + raise ValueError("If --source is passed then --target-lsm must also be given.") -def _transform_coordinates(target_lsm, source_cube, points): +def _transform_coordinates(target_lsm, source, points): """ Transform the longitude-latitude points to the rotated pole. - If the distance between transformed longitude points exceeds 180.0, - it is assumed these points cross the antimeridian. The points are - changed so they are instead defined in the interval [0.0, 360.0]. + If the transformation results in the points spanning the antimeridian, this + can cause incorrect polygons to be created because shapely does not wrap + between 180.0 and -180.0 degrees. + + If the distance between sequential transformed longitude points exceeds + 180.0 degrees, it is assumed these points cross the antimeridian. The points + are instead defined in the interval [0.0, 360.0] by adding 360.0 degrees + to the negative longitudes. Parameters ---------- target_lsm : :class:`iris.cube.Cube` The lsm cube specifying the rotated pole coordinate system. - source_cube : :class:`iris.cube.Cube` + source : :class:`iris.cube.Cube` An iris cube specifying the coordinate system of the input json file. points : :class:`numpy.ndarray` @@ -159,7 +162,7 @@ def _transform_coordinates(target_lsm, source_cube, points): An ``(m, 2)`` numpy array of transformed longitude-latitude pairs. """ - source_crs = source_cube.coord_system().as_cartopy_crs() + source_crs = source.coord_system().as_cartopy_crs() target_coord = target_lsm.coord_system() target_crs = target_coord.as_cartopy_crs() @@ -167,10 +170,17 @@ def _transform_coordinates(target_lsm, source_cube, points): source_crs, points[:, 0], points[:, 1] )[:, :2] + # Create a wrapped array of longitudes and find the point-wise difference. closed_lon = np.vstack([rotated_points, rotated_points[0, :]]) lon_diff = np.abs(closed_lon[:-1, 0] - closed_lon[1:, 0]) + # Add 360.0 degrees to negative longitudes if the points cross the antimeridian. if np.any(lon_diff > 180.0): + warnings.warn( + "The transformed points are assumed to cross the antimeridian. " + "The longitudinal points will instead be defined in the interval" + " [0, 360.0] degrees." + ) neg_indices = np.where(rotated_points[:, 0] < 0) rotated_points[neg_indices, 0] += 360.0 @@ -185,7 +195,7 @@ def _transform_coordinates(target_lsm, source_cube, points): return rotated_points -def _load_cubes(target_lsm_path, source_cube_path): +def _load_cubes(target_lsm_path, source_path): """ Load the target lsm and create a source cube if not provided. @@ -193,8 +203,8 @@ def _load_cubes(target_lsm_path, source_cube_path): ---------- target_lsm_path : str File path to a land sea mask that provides the new rotated pole. - source_cube_path : str - File path to an iris cube specifying the coordinate system of the + source_path : str + File path to a source file specifying the coordinate system of the input json file. Returns @@ -205,16 +215,16 @@ def _load_cubes(target_lsm_path, source_cube_path): target_lsm = load_landsea_mask(target_lsm_path) - if source_cube_path is None: + if source_path is None: crs = iris.coord_systems.GeogCS(6371229.0) - source_cube = CubeBuilder(crs, (2, 2))._cube + source = CubeBuilder(crs, (2, 2))._cube else: - source_cube = load_cube(source_cube_path) + source = load_cube(source_path) - return target_lsm, source_cube + return target_lsm, source -def _transform_if_required(target_lsm, source_cube, points): +def _transform_if_required(target_lsm, source, points): """ Perform the transformation to a rotated pole if the coordinate system is valid. @@ -227,7 +237,7 @@ def _transform_if_required(target_lsm, source_cube, points): ---------- target_lsm : :class:`iris.cube.Cube` The lsm cube specifying the rotated pole coordinate system. - source_cube : :class:`iris.cube.Cube` + source : :class:`iris.cube.Cube` An iris cube specifying the coordinate system of the input json file. points : :class:`numpy.ndarray` @@ -242,7 +252,7 @@ def _transform_if_required(target_lsm, source_cube, points): _check_coord_system_type(target_lsm) if _validate_orientation(target_lsm): - rotated_points = _transform_coordinates(target_lsm, source_cube, points) + rotated_points = _transform_coordinates(target_lsm, source, points) else: rotated_points = np.copy(points) @@ -270,7 +280,7 @@ def _load_points_from_json(json_file): return points -def main(json_file, output, target_lsm_path, source_cube_path): +def main(json_file, output, target_lsm_path, source_path): """ Create a shape file from pairs of longitude, latitude points. @@ -280,8 +290,8 @@ def main(json_file, output, target_lsm_path, source_cube_path): If target_lsm_path is provided, the points are first transformed from a source geodetic coordinate system to a rotated pole coordinate system - specified by the lsm. It is assumed that the points in the json - file are on an unrotated geodetic grid, unless otherwise specified. + specified by the lsm. It is assumed that the points in the json file + are on an unrotated spherical geodetic grid, unless otherwise specified. Parameters ---------- @@ -291,8 +301,8 @@ def main(json_file, output, target_lsm_path, source_cube_path): Location to store generated shape file target_lsm_path : str File path to a land sea mask that provides the new rotated pole. - source_cube_path : str - File path to an iris cube specifying the coordinate system of the + source_path : str + File path to a source file specifying the coordinate system of the input json file. """ @@ -301,8 +311,8 @@ def main(json_file, output, target_lsm_path, source_cube_path): ccw_expected = Polygon(points).exterior.is_ccw if target_lsm_path is not None: - target_lsm, source_cube = _load_cubes(target_lsm_path, source_cube_path) - points = _transform_if_required(target_lsm, source_cube, points) + target_lsm, source = _load_cubes(target_lsm_path, source_path) + points = _transform_if_required(target_lsm, source, points) polygon = Polygon(points) @@ -346,11 +356,10 @@ def _get_parser(): " coordinate system.", ) parser.add_argument( - "--source-cube", + "--source", type=ants.config.filepath_readable, required=False, - help="Path to an iris cube which specifies the coordinate" - " system of the json file.", + help="Path to a source containing the coordinate system of the json file.", ) return parser @@ -359,8 +368,8 @@ def cli_interface(): parser = _get_parser() args = parser.parse_args() - _validate_args(args.target_lsm, args.source_cube) - main(args.json_file, args.output, args.target_lsm, args.source_cube) + _validate_args(args.target_lsm, args.source) + main(args.json_file, args.output, args.target_lsm, args.source) if __name__ == "__main__": diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index 4b9a8f1..2de4d07 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -14,7 +14,6 @@ _check_coord_system_type, _check_polygon_validity, _load_cubes, - _load_points_from_json, _transform_coordinates, _transform_if_required, _validate_args, @@ -79,28 +78,21 @@ def setUp(self): self.sphere_identity_crs = iris.coord_systems.RotatedGeogCS( 90.0, 0.0, ellipsoid=self.sphere_crs ) - self.sphere_rotated_crs = iris.coord_systems.RotatedGeogCS( - 90.0, 0.0, 180.0, ellipsoid=self.sphere_crs - ) + + def test_warning_raised(self): + """Test that a warning is raised with a coordinate system at latitude=90.0, + longitude=0.0 is passed. Check the function returns False.""" + self.sphere_identity_target_lsm = CubeBuilder( self.sphere_identity_crs, (2, 2) )._cube - self.sphere_rotated_target_lsm = CubeBuilder( - self.sphere_rotated_crs, (2, 2) - )._cube - self.sphere_rotate_lon = geodetic( - (2, 2), north_pole_lat=90.0, north_pole_lon=90.0, crs=self.sphere_crs - ) + self.warning_msg = ( "target_lsm has a geodetic coordinate system with pole located" " at grid_longitude=0.0, grid_latitude=90.0." " No transformation will be carried out." ) - def test_warning_raised(self): - """Test that a warning is raised with a coordinate system at latitude=90.0, - longitude=0.0 is passed. Check the function returns False.""" - with self.assertWarnsRegex(UserWarning, self.warning_msg): valid_crs = _validate_orientation(self.sphere_identity_target_lsm) self.assertFalse(valid_crs) @@ -108,12 +100,23 @@ def test_warning_raised(self): def test_true_returned(self): """Test that True is returned when target_lsm has a valid rotated pole.""" + self.sphere_rotate_lon = geodetic( + (2, 2), north_pole_lat=90.0, north_pole_lon=90.0, crs=self.sphere_crs + ) + self.assertTrue(_validate_orientation(self.sphere_rotate_lon)) def test_no_warning_raised_rotated(self): """Test that passing a non-zero central rotated longitude is still accounted for when checking if the coordinate system is rotated.""" + self.sphere_rotated_crs = iris.coord_systems.RotatedGeogCS( + 90.0, 0.0, 180.0, ellipsoid=self.sphere_crs + ) + self.sphere_rotated_target_lsm = CubeBuilder( + self.sphere_rotated_crs, (2, 2) + )._cube + self.assertTrue(_validate_orientation(self.sphere_rotated_target_lsm)) @@ -123,25 +126,7 @@ def setUp(self): self.sphere_identity_crs = iris.coord_systems.RotatedGeogCS( 90.0, 0.0, ellipsoid=self.sphere_crs ) - self.sphere_identity_target_lsm = CubeBuilder( - self.sphere_identity_crs, (2, 2) - )._cube self.sphere_source = CubeBuilder(self.sphere_crs, (2, 2))._cube - self.sphere_equator_target_lsm = geodetic( - (2, 2), north_pole_lat=0.0, north_pole_lon=0.0, crs=self.sphere_crs - ) - self.sphere_rotate_lon = geodetic( - (2, 2), north_pole_lat=90.0, north_pole_lon=90.0, crs=self.sphere_crs - ) - self.vary_latitudes = np.array( - [[0.0, 0.0, 0.0, 0.0, 0.0], [90.0, 45.0, 0.0, -45.0, -90.0]] - ).T - self.vary_longitudes = np.array( - [ - [-135.0, -90.0, -45.0, 0.0, 45.0, 90.0, 135.0], - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], - ] - ).T def test_identity_rotation_sphere(self): """Test that rotation to a pole at latitude=90.0, longitude=0.0 @@ -149,6 +134,13 @@ def test_identity_rotation_sphere(self): convention the new pole is defined 180.0 rotated from the provided longitude.""" + self.sphere_identity_target_lsm = CubeBuilder( + self.sphere_identity_crs, (2, 2) + )._cube + + self.vary_latitudes = np.array( + [[0.0, 0.0, 0.0, 0.0, 0.0], [90.0, 45.0, 0.0, -45.0, -90.0]] + ).T true_lats = np.array([90, 45, 0, -45, -90]) true_lons = np.array([np.nan, -180.0, -180.0, -180.0, -180.0]) expected_rotation = np.array([true_lons, true_lats]).T @@ -163,11 +155,17 @@ def test_identity_rotation_sphere(self): expected_rotation[check_mask], rotated_coords[check_mask] ) + # TODO this test doesn't really make sense. These are really just a + # line of points and not a polygon. def test_longitudinal_rotation(self): """Test rotation to a new pole at latitude=90.0, longitude=90.0. The new pole is located at latitude=90.0, longitude=270.0 in the rotated pole coordinate system.""" + self.sphere_rotate_lon = geodetic( + (2, 2), north_pole_lat=90.0, north_pole_lon=90.0, crs=self.sphere_crs + ) + points = np.array([[0.0, 45.0, 90.0, 135.0], [0.0, 0.0, 0.0, 0.0]]).T true_lats = np.array([0.0, 0.0, np.nan, 0.0]) true_lons = np.array([90.0, 135.0, 180.0, 225.0]) @@ -185,11 +183,22 @@ def test_longitudinal_rotation(self): def test_rotation_to_equator_sphere(self): """Test rotation to a new pole at latitude=0.0, longitude=0.0.""" + self.vary_longitudes = np.array( + [ + [-135.0, -90.0, -45.0, 0.0, 45.0, 90.0, 135.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + ] + ).T + + self.sphere_equator_target_lsm = geodetic( + (2, 2), north_pole_lat=0.0, north_pole_lon=0.0, crs=self.sphere_crs + ) true_lats = np.array([-45.0, 0.0, 45.0, 90.0, 45.0, 0.0, -45.0]) true_lons = np.array([90.0, 90.0, 90.0, np.nan, -90.0, -90.0, -90.0]) expected_rotation = np.array([true_lons, true_lats]).T check_mask = ~np.isnan(expected_rotation) + rotated_coords = _transform_coordinates( self.sphere_equator_target_lsm, self.sphere_source, self.vary_longitudes ) @@ -201,17 +210,17 @@ def test_rotation_to_equator_sphere(self): class Test__validate_args(ants.tests.TestCase): def test_error_raised(self): - """Test that an error is raised if only --source-cube is passed.""" + """Test that an error is raised if only --source is passed.""" args = argparse.Namespace( json_file="json/path", output="output/path", target_lsm=None, - source_cube="source/path", + source="source/path", ) - error_msg = "If --source-cube is passed then --target-lsm must also be given." + error_msg = "If --source is passed then --target-lsm must also be given." with self.assertRaisesRegex(ValueError, error_msg): - _validate_args(args.target_lsm, args.source_cube) + _validate_args(args.target_lsm, args.source) class Test__load_cubes(ants.tests.TestCase): @@ -251,9 +260,10 @@ def setUp(self): self.points = np.array([[1, 2], [3, 4], [5, 6], [7, 8]]) self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) self.sphere_source = CubeBuilder(self.sphere_crs, (2, 2))._cube - self.sphere_equator_target_lsm = geodetic( - (2, 2), north_pole_lat=0.0, north_pole_lon=0.0, crs=self.sphere_crs - ) + + def test_non_rotated_returned(self): + """Test that the input points are returned if a non-rotated pole is + provided as the target_lsm.""" self.sphere_identity_crs = iris.coord_systems.RotatedGeogCS( 90.0, 0.0, ellipsoid=self.sphere_crs ) @@ -261,10 +271,6 @@ def setUp(self): self.sphere_identity_crs, (2, 2) )._cube - def test_non_rotated_returned(self): - """Test that the input points are returned if a non-rotated pole is - provided as the target_lsm.""" - with self.assertWarns(UserWarning): rotated_points = _transform_if_required( self.sphere_identity_target_lsm, self.sphere_source, self.points @@ -276,6 +282,9 @@ def test_transform_called(self): """Test that transform coordinates is called with the correct arguments when a valid target_lsm is given.""" + self.sphere_equator_target_lsm = geodetic( + (2, 2), north_pole_lat=0.0, north_pole_lon=0.0, crs=self.sphere_crs + ) with mock.patch( "ants.cli.ancil_create_shapefile._transform_coordinates" ) as mock_transform: @@ -283,22 +292,8 @@ def test_transform_called(self): self.sphere_equator_target_lsm, self.sphere_source, self.points ) - target_lsm, source_cube, points = mock_transform.call_args.args + target_lsm, source, points = mock_transform.call_args.args mock_transform.assert_called_once() self.assertEqual(target_lsm, self.sphere_equator_target_lsm) - self.assertEqual(source_cube, self.sphere_source) + self.assertEqual(source, self.sphere_source) self.assertArrayEqual(points, self.points) - - -@mock.patch("builtins.open", new_callable=mock.mock_open) -class Test__load_points_from_json(ants.tests.TestCase): - def setUp(self): - self.json_values = [[1, 2], [3, 4], [5, 6], [7, 8]] - - def test_json_load(self, *args): - """Test that a json file is loaded.""" - - with mock.patch("ants.cli.ancil_create_shapefile.json.load") as mock_json: - _ = _load_points_from_json("json/path") - - mock_json.assert_called_once() From 8e9ae5b262dea377115e0655f1d93a0cbc27ce86 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Thu, 20 Aug 2026 16:36:24 +0100 Subject: [PATCH 24/29] Added integration test for unrotated ITE. --- lib/ants/cli/ancil_create_shapefile.py | 2 +- .../tests/cli/test_ancil_create_shapefile.py | 67 +++++++++++++++++++ 2 files changed, 68 insertions(+), 1 deletion(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index 3f9b8dc..db366b3 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -11,7 +11,7 @@ points defining a single polygon in a specified polygon file. Rotated pole domains can be specified using the land sea mask argument, -where the longitude, latitude pairs are rotated to the new pole. +where the longitude, latitude pairs are transformed to the new pole. Unless a source cube is provided with the source coordinate reference system, it assumed that the points defined in the json file are on a standard spherical unrotated geodetic coordinate reference system. diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index 2de4d07..a62c8e9 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -297,3 +297,70 @@ def test_transform_called(self): self.assertEqual(target_lsm, self.sphere_equator_target_lsm) self.assertEqual(source, self.sphere_source) self.assertArrayEqual(points, self.points) + + +class Test_ite_transform(ants.tests.TestCase): + + def test_uk_pole_rotation(self): + """ + Test the unrotated polygon enclosing regions of valid data from ITE + expressed in true longitude and latitude is rotated correctly. The target + pole is lon=177.5, lat=37.5. + + The points have been obtained by unrotating the coordinates in + $UMDIR/ancil/data/shapefiles/ite_ukv_polygon/runme.py + """ + self.ite_points = np.array( + [ + [1.63160953, 51.09703216], + [-0.29556461, 50.36332573], + [-5.37076475, 49.91322904], + [-6.01981579, 50.15822986], + [-5.16157337, 53.5014088], + [-3.80608343, 53.99550534], + [-4.04624456, 54.56480733], + [-5.10119435, 54.47567927], + [-5.95760936, 55.29939785], + [-6.998891, 55.85406809], + [-7.98467448, 56.73347598], + [-7.65168129, 58.36692001], + [-3.71924355, 61.09496898], + [0.11346329, 61.0768747], + [2.0703864, 52.72809453], + [1.7594481, 51.23560201], + ] + ) + + self.expected_rotated_points = np.array( + [ + [362.594, -1.32876], + [361.407, -2.1152], + [358.15, -2.55], + [357.744, -2.28679], + [358.417, 1.03058], + [359.232, 1.50245], + [359.103, 2.07441], + [358.488, 2.00291], + [358.03, 2.84656], + [357.472, 3.43282], + [356.986, 4.34795], + [357.286, 5.96374], + [359.404, 8.6], + [361.278, 8.6], + [362.766, 0.315629], + [362.666, -1.18577], + ] + ) + self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) + self.sphere_source = CubeBuilder(self.sphere_crs, (2, 2))._cube + + target_lsm = geodetic( + (2, 2), north_pole_lat=37.5, north_pole_lon=177.5, crs=self.sphere_crs + ) + + points = _transform_if_required(target_lsm, self.sphere_source, self.ite_points) + + # Adding 360.0 to compare to the expected points + points[:, 0] += 360.0 + + self.assertArrayAlmostEqual(points, self.expected_rotated_points) From b76caf7a6a069e012f03e6479b6c9cea69087902 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Fri, 21 Aug 2026 13:28:17 +0100 Subject: [PATCH 25/29] Added new ite test and made other tests more consistent. --- .../tests/cli/test_ancil_create_shapefile.py | 198 ++++++++---------- 1 file changed, 85 insertions(+), 113 deletions(-) diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index a62c8e9..f2e73c2 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -47,21 +47,21 @@ def test_invalid_polygon(self): """Test that an invalid polygon with intersecting boundaries raises an error.""" - self.invalid_points = np.array([[1, 1], [-1, -1], [1, -1], [-1, 1]]) - self.invalid_poly = Polygon(self.invalid_points) - self.is_ccw = False + invalid_points = np.array([[1, 1], [-1, -1], [1, -1], [-1, 1]]) + invalid_poly = Polygon(invalid_points) + is_ccw = False error_msg = "Polygon is invalid: " with self.assertRaisesRegex(ValueError, error_msg): - _check_polygon_validity(self.invalid_poly, self.is_ccw) + _check_polygon_validity(invalid_poly, is_ccw) def test_changed_orientation(self): """Test that a change in polygon orientation raises an error.""" - self.ccw_points = np.array([[1, 1], [-1, 1], [-1, -1], [1, -1]]) - self.ccw_poly = Polygon(self.ccw_points) - self.is_ccw = False + ccw_points = np.array([[1, 1], [-1, 1], [-1, -1], [1, -1]]) + ccw_poly = Polygon(ccw_points) + is_ccw = False error_msg = ( "Polygon orientation has changed. Expected is_ccw=False, " @@ -69,7 +69,7 @@ def test_changed_orientation(self): ) with self.assertRaisesRegex(ValueError, error_msg): - _check_polygon_validity(self.ccw_poly, self.is_ccw) + _check_polygon_validity(ccw_poly, is_ccw) class Test__validate_orientation(ants.tests.TestCase): @@ -83,129 +83,108 @@ def test_warning_raised(self): """Test that a warning is raised with a coordinate system at latitude=90.0, longitude=0.0 is passed. Check the function returns False.""" - self.sphere_identity_target_lsm = CubeBuilder( - self.sphere_identity_crs, (2, 2) - )._cube + target_lsm = CubeBuilder(self.sphere_identity_crs, (2, 2))._cube - self.warning_msg = ( + warning_msg = ( "target_lsm has a geodetic coordinate system with pole located" " at grid_longitude=0.0, grid_latitude=90.0." " No transformation will be carried out." ) - with self.assertWarnsRegex(UserWarning, self.warning_msg): - valid_crs = _validate_orientation(self.sphere_identity_target_lsm) + with self.assertWarnsRegex(UserWarning, warning_msg): + valid_crs = _validate_orientation(target_lsm) self.assertFalse(valid_crs) def test_true_returned(self): """Test that True is returned when target_lsm has a valid rotated pole.""" - self.sphere_rotate_lon = geodetic( + target_lsm = geodetic( (2, 2), north_pole_lat=90.0, north_pole_lon=90.0, crs=self.sphere_crs ) - self.assertTrue(_validate_orientation(self.sphere_rotate_lon)) + self.assertTrue(_validate_orientation(target_lsm)) def test_no_warning_raised_rotated(self): """Test that passing a non-zero central rotated longitude is still accounted for when checking if the coordinate system is rotated.""" - self.sphere_rotated_crs = iris.coord_systems.RotatedGeogCS( + crs = iris.coord_systems.RotatedGeogCS( 90.0, 0.0, 180.0, ellipsoid=self.sphere_crs ) - self.sphere_rotated_target_lsm = CubeBuilder( - self.sphere_rotated_crs, (2, 2) - )._cube + target_lsm = CubeBuilder(crs, (2, 2))._cube - self.assertTrue(_validate_orientation(self.sphere_rotated_target_lsm)) + self.assertTrue(_validate_orientation(target_lsm)) class Test__transform_coordinates(ants.tests.TestCase): def setUp(self): self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) - self.sphere_identity_crs = iris.coord_systems.RotatedGeogCS( - 90.0, 0.0, ellipsoid=self.sphere_crs - ) self.sphere_source = CubeBuilder(self.sphere_crs, (2, 2))._cube def test_identity_rotation_sphere(self): """Test that rotation to a pole at latitude=90.0, longitude=0.0 - returns the same points, except a 180.0 degree rotation in longitude. By - convention the new pole is defined 180.0 rotated from the provided - longitude.""" - - self.sphere_identity_target_lsm = CubeBuilder( - self.sphere_identity_crs, (2, 2) - )._cube - - self.vary_latitudes = np.array( - [[0.0, 0.0, 0.0, 0.0, 0.0], [90.0, 45.0, 0.0, -45.0, -90.0]] - ).T - true_lats = np.array([90, 45, 0, -45, -90]) - true_lons = np.array([np.nan, -180.0, -180.0, -180.0, -180.0]) - expected_rotation = np.array([true_lons, true_lats]).T - - # We cannot check the longitude at the pole (any longitude is valid). - check_mask = ~np.isnan(expected_rotation) - - rotated_coords = _transform_coordinates( - self.sphere_identity_target_lsm, self.sphere_source, self.vary_latitudes - ) - self.assertArrayAlmostEqual( - expected_rotation[check_mask], rotated_coords[check_mask] - ) + returns the same points. - # TODO this test doesn't really make sense. These are really just a - # line of points and not a polygon. - def test_longitudinal_rotation(self): - """Test rotation to a new pole at latitude=90.0, longitude=90.0. - The new pole is located at latitude=90.0, longitude=270.0 in the rotated - pole coordinate system.""" + By convention, rotated pole co-ordinate systems will set the prime + meridian rotated 180.0 from the specified longitude. To place the + prime meridian at 0.0, we apply a further rotation of 180.0, following + rotation to the new pole. + """ - self.sphere_rotate_lon = geodetic( - (2, 2), north_pole_lat=90.0, north_pole_lon=90.0, crs=self.sphere_crs + crs = iris.coord_systems.RotatedGeogCS( + 90.0, 0.0, 180.0, ellipsoid=self.sphere_crs ) + target_lsm = CubeBuilder(crs, (2, 2))._cube + + points = np.array([[10, 10], [10, -10], [-10, -10], [-10, 10]]) - points = np.array([[0.0, 45.0, 90.0, 135.0], [0.0, 0.0, 0.0, 0.0]]).T - true_lats = np.array([0.0, 0.0, np.nan, 0.0]) - true_lons = np.array([90.0, 135.0, 180.0, 225.0]) - expected_rotation = np.array([true_lons, true_lats]).T + expected_points = np.copy(points) - check_mask = ~np.isnan(expected_rotation) + rotated_coords = _transform_coordinates(target_lsm, self.sphere_source, points) + self.assertArrayAlmostEqual(expected_points, rotated_coords) + + def test_longitudinal_rotation(self): + """Test rotation to a new pole at latitude=90.0, longitude=90.0.""" - rotated_coords = _transform_coordinates( - self.sphere_rotate_lon, self.sphere_source, points + crs = iris.coord_systems.RotatedGeogCS( + 90.0, 90.0, 180.0, ellipsoid=self.sphere_crs ) + target_lsm = CubeBuilder(crs, (2, 2))._cube - self.assertArrayAlmostEqual( - expected_rotation[check_mask], rotated_coords[check_mask] + points = np.array([[10, 10], [10, -10], [-10, -10], [-10, 10]]) + expected_points = np.array( + [[-80.0, 10.0], [-80.0, -10.0], [-100.0, -10.0], [-100.0, 10]] ) - def test_rotation_to_equator_sphere(self): - """Test rotation to a new pole at latitude=0.0, longitude=0.0.""" - self.vary_longitudes = np.array( - [ - [-135.0, -90.0, -45.0, 0.0, 45.0, 90.0, 135.0], - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], - ] - ).T + rotated_coords = _transform_coordinates(target_lsm, self.sphere_source, points) - self.sphere_equator_target_lsm = geodetic( - (2, 2), north_pole_lat=0.0, north_pole_lon=0.0, crs=self.sphere_crs - ) + self.assertArrayAlmostEqual(expected_points, rotated_coords) - true_lats = np.array([-45.0, 0.0, 45.0, 90.0, 45.0, 0.0, -45.0]) - true_lons = np.array([90.0, 90.0, 90.0, np.nan, -90.0, -90.0, -90.0]) - expected_rotation = np.array([true_lons, true_lats]).T - check_mask = ~np.isnan(expected_rotation) + def test_antimeridian_rotation(self): + """Test rotation to a new pole at latitude=90.0, longitude=180.0.""" - rotated_coords = _transform_coordinates( - self.sphere_equator_target_lsm, self.sphere_source, self.vary_longitudes + crs = iris.coord_systems.RotatedGeogCS( + 90.0, 180.0, 180.0, ellipsoid=self.sphere_crs ) - self.assertArrayAlmostEqual( - expected_rotation[check_mask], rotated_coords[check_mask] + target_lsm = CubeBuilder(crs, (2, 2))._cube + + points = np.array([[10, 10], [10, -10], [-10, -10], [-10, 10]]) + expected_points = np.array( + [[190.0, 10.0], [190.0, -10.0], [170.0, -10.0], [170.0, 10]] ) + warning_msg = re.escape( + "The transformed points are assumed to cross the antimeridian. " + "The longitudinal points will instead be defined in the interval" + " [0, 360.0] degrees." + ) + with self.assertWarnsRegex(UserWarning, warning_msg): + rotated_coords = _transform_coordinates( + target_lsm, self.sphere_source, points + ) + + self.assertArrayAlmostEqual(expected_points, rotated_coords) + class Test__validate_args(ants.tests.TestCase): @@ -264,16 +243,13 @@ def setUp(self): def test_non_rotated_returned(self): """Test that the input points are returned if a non-rotated pole is provided as the target_lsm.""" - self.sphere_identity_crs = iris.coord_systems.RotatedGeogCS( - 90.0, 0.0, ellipsoid=self.sphere_crs - ) - self.sphere_identity_target_lsm = CubeBuilder( - self.sphere_identity_crs, (2, 2) - )._cube + + crs = iris.coord_systems.RotatedGeogCS(90.0, 0.0, ellipsoid=self.sphere_crs) + target_lsm = CubeBuilder(crs, (2, 2))._cube with self.assertWarns(UserWarning): rotated_points = _transform_if_required( - self.sphere_identity_target_lsm, self.sphere_source, self.points + target_lsm, self.sphere_source, self.points ) self.assertArrayEqual(rotated_points, self.points) @@ -282,21 +258,19 @@ def test_transform_called(self): """Test that transform coordinates is called with the correct arguments when a valid target_lsm is given.""" - self.sphere_equator_target_lsm = geodetic( - (2, 2), north_pole_lat=0.0, north_pole_lon=0.0, crs=self.sphere_crs - ) + crs = iris.coord_systems.RotatedGeogCS(0.0, 0.0, ellipsoid=self.sphere_crs) + target_lsm = CubeBuilder(crs, (2, 2))._cube + with mock.patch( "ants.cli.ancil_create_shapefile._transform_coordinates" ) as mock_transform: - _ = _transform_if_required( - self.sphere_equator_target_lsm, self.sphere_source, self.points - ) + _ = _transform_if_required(target_lsm, self.sphere_source, self.points) - target_lsm, source, points = mock_transform.call_args.args + received_lsm, received_source, received_points = mock_transform.call_args.args mock_transform.assert_called_once() - self.assertEqual(target_lsm, self.sphere_equator_target_lsm) - self.assertEqual(source, self.sphere_source) - self.assertArrayEqual(points, self.points) + self.assertEqual(target_lsm, received_lsm) + self.assertEqual(received_source, self.sphere_source) + self.assertArrayEqual(received_points, self.points) class Test_ite_transform(ants.tests.TestCase): @@ -310,7 +284,7 @@ def test_uk_pole_rotation(self): The points have been obtained by unrotating the coordinates in $UMDIR/ancil/data/shapefiles/ite_ukv_polygon/runme.py """ - self.ite_points = np.array( + points = np.array( [ [1.63160953, 51.09703216], [-0.29556461, 50.36332573], @@ -331,7 +305,7 @@ def test_uk_pole_rotation(self): ] ) - self.expected_rotated_points = np.array( + expected_points = np.array( [ [362.594, -1.32876], [361.407, -2.1152], @@ -351,16 +325,14 @@ def test_uk_pole_rotation(self): [362.666, -1.18577], ] ) - self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) - self.sphere_source = CubeBuilder(self.sphere_crs, (2, 2))._cube - - target_lsm = geodetic( - (2, 2), north_pole_lat=37.5, north_pole_lon=177.5, crs=self.sphere_crs - ) + sphere_crs = iris.coord_systems.GeogCS(6371229.0) + source = CubeBuilder(sphere_crs, (2, 2))._cube + crs = iris.coord_systems.RotatedGeogCS(37.5, 177.5, ellipsoid=sphere_crs) + target_lsm = CubeBuilder(crs, (2, 2))._cube - points = _transform_if_required(target_lsm, self.sphere_source, self.ite_points) + rotated_points = _transform_if_required(target_lsm, source, points) - # Adding 360.0 to compare to the expected points - points[:, 0] += 360.0 + # Adding 360.0 to longitude to compare to the expected points + rotated_points[:, 0] += 360.0 - self.assertArrayAlmostEqual(points, self.expected_rotated_points) + self.assertArrayAlmostEqual(rotated_points, expected_points) From 12f4707fd3d3c641110749c51a27e27f3fc14eaa Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Tue, 25 Aug 2026 10:06:12 +0100 Subject: [PATCH 26/29] Added warning if the polygon lands outside of the local domain. --- lib/ants/cli/ancil_create_shapefile.py | 24 ++ ...test_moore_neighbourhood_search.py.isorted | 281 ++++++++++++++++++ .../tests/cli/test_ancil_create_shapefile.py | 78 +++-- 3 files changed, 365 insertions(+), 18 deletions(-) create mode 100644 lib/ants/tests/analysis/merge/test_moore_neighbourhood_search.py.isorted diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index db366b3..89fbef6 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -160,12 +160,20 @@ def _transform_coordinates(target_lsm, source, points): ------- : :class:`numpy.ndarray` An ``(m, 2)`` numpy array of transformed longitude-latitude pairs. + Warns + ----- + UserWarning + If the transformed points are assumed to cross the antimeridian. + UserWarning + If the transformed polygon lies at least partially outside of the + domain specified in the target lsm. """ source_crs = source.coord_system().as_cartopy_crs() target_coord = target_lsm.coord_system() target_crs = target_coord.as_cartopy_crs() + # We only return longitude and latitude and discard the z coordinate. rotated_points = target_crs.transform_points( source_crs, points[:, 0], points[:, 1] )[:, :2] @@ -184,6 +192,22 @@ def _transform_coordinates(target_lsm, source, points): neg_indices = np.where(rotated_points[:, 0] < 0) rotated_points[neg_indices, 0] += 360.0 + # Find the bounds of the target lsm. + bounds_lon = target_lsm.coord(axis="X").bounds + bounds_lat = target_lsm.coord(axis="Y").bounds + min_lon, max_lon = bounds_lon.min(), bounds_lon.max() + min_lat, max_lat = bounds_lat.min(), bounds_lat.max() + lons, lats = rotated_points[:, 0], rotated_points[:, 1] + + # Check if the polygon lies in the domain specified by the target lsm. + if ( + lons.min() < min_lon + or lons.max() > max_lon + or lats.min() < min_lat + or lats.max() > max_lat + ): + warnings.warn("The transformed points lie outside the target lsm domain.") + _LOGGER.info( "Input json file transformed to new pole rotated coordinate system at " "pole longitude=%s, pole latitude=%s, central rotated longitude=%s.", diff --git a/lib/ants/tests/analysis/merge/test_moore_neighbourhood_search.py.isorted b/lib/ants/tests/analysis/merge/test_moore_neighbourhood_search.py.isorted new file mode 100644 index 0000000..1a5e120 --- /dev/null +++ b/lib/ants/tests/analysis/merge/test_moore_neighbourhood_search.py.isorted @@ -0,0 +1,281 @@ +# (C) Crown Copyright, Met Office. All rights reserved. +# +# This file is part of ANTS and is released under the BSD 3-Clause license. +# See LICENSE.txt in the root of the repository for full licensing details. +import ants.tests +import iris +import numpy as np +from ants.analysis._merge import moore_neighbourhood_search + + +class _MaskApplication(object): + def get_cube(self, data): + # Return a global cube with -1 values the assumed masked values + cube = ants.tests.stock.geodetic(data.shape, with_bounds=True) + cube.coord(axis="x").circular = False + data = np.ma.masked_values(data, -1, copy=False) + cube.data = data + return cube + + def _assert_expected(self, cube, target, land_binary_mask): + self.assertArrayEqual(cube.data.data, target) + self.assertArrayEqual(cube.data.mask, land_binary_mask.data == 0) + + +class Test_inappropriate_input(_MaskApplication, ants.tests.TestCase): + def test_not_2d_source(self): + cube = iris.cube.Cube(np.arange(8).reshape(2, 2, 2)) + land_binary_mask = iris.cube.Cube(np.arange(8).reshape(2, 2, 2)) + msg = ( + "Currently, only mask application to 2D grids are supported " + "and with no broadcasting." + ) + with self.assertRaisesRegex(RuntimeError, msg): + moore_neighbourhood_search(cube, land_binary_mask) + + def test_different_grids(self): + data = np.arange(4).reshape(2, 2) + cube = self.get_cube(data) + land_binary_mask = self.get_cube(data) + land_binary_mask.coord(axis="x").points = cube.coord(axis="x").points - 10 + msg = "Both source and land_binary_mask must be defined on identical " "grids" + with self.assertRaisesRegex(RuntimeError, msg): + moore_neighbourhood_search(cube, land_binary_mask) + + def test_alternate_orientation(self): + # Currently we do not support automatically transposing the data where + # coordinates are mapped to alternate dimensions. Make sure we catch + # this case. + data = np.arange(4).reshape(2, 2) + cube = self.get_cube(data) + land_binary_mask = self.get_cube(data) + lx = land_binary_mask.coord(axis="x") + ly = land_binary_mask.coord(axis="y") + land_binary_mask.remove_coord(lx) + land_binary_mask.remove_coord(ly) + land_binary_mask.add_dim_coord(lx, 0) + land_binary_mask.add_dim_coord(ly, 1) + + msg = ( + "Currently, the source and the land_binary_mask must be " + "defined in the same orientation" + ) + with self.assertRaisesRegex(RuntimeError, msg): + moore_neighbourhood_search(cube, land_binary_mask) + + +class Test_missing_coastlines(_MaskApplication, ants.tests.TestCase): + def test_no_change(self): + data = np.ma.array([[-1, 2, 3], [4, -1, -1], [7, 8, 9]], dtype=float) + cube = self.get_cube(data) + land_binary_mask = cube.copy((data > 0).astype(int)) + + moore_neighbourhood_search(cube, land_binary_mask, 10) + + target = data.copy() + self._assert_expected(cube, target, land_binary_mask) + + def test_missing_top_left_edge(self): + data = np.ma.array([[-1, 2, 3], [2, 5, 6], [7, 8, 9]], dtype=float) + cube = self.get_cube(data) + land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) + + moore_neighbourhood_search(cube, land_binary_mask, 10) + + target = data.data.copy() + target[0, 0] = 3 + self._assert_expected(cube, target, land_binary_mask) + + def test_missing_top_centre_edge(self): + data = np.ma.array([[4, -1, 3], [2, 5, 6], [7, 8, 9]], dtype=float) + cube = self.get_cube(data) + land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) + + moore_neighbourhood_search(cube, land_binary_mask, 10) + + target = data.data.copy() + target[0, 1] = 4 + self._assert_expected(cube, target, land_binary_mask) + + def test_missing_top_right_edge(self): + data = np.ma.array([[4, 2, -1], [2, 5, 2], [7, 8, 9]], dtype=float) + cube = self.get_cube(data) + land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) + + moore_neighbourhood_search(cube, land_binary_mask, 10) + + target = data.data.copy() + target[0, 2] = 3 + self._assert_expected(cube, target, land_binary_mask) + + def test_missing_centre_left_edge(self): + data = np.ma.array([[4, 2, 2], [-1, 5, 2], [3, 6, 9]], dtype=float) + cube = self.get_cube(data) + land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) + + moore_neighbourhood_search(cube, land_binary_mask, 10) + + target = data.data.copy() + target[1, 0] = 4 + self._assert_expected(cube, target, land_binary_mask) + + def test_missing_centre_centre(self): + data = np.ma.array([[2, 4, 2], [4, -1, 4], [2, 4, 2]], dtype=float) + cube = self.get_cube(data) + land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) + + moore_neighbourhood_search(cube, land_binary_mask, 10) + + target = data.data.copy() + target[1, 1] = 3 + self._assert_expected(cube, target, land_binary_mask) + + def test_missing_centre_right_edge(self): + data = np.ma.array([[2, 4, 2], [4, 3, -1], [2, 4, 2]], dtype=float) + cube = self.get_cube(data) + land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) + + moore_neighbourhood_search(cube, land_binary_mask, 10) + + target = data.data.copy() + target[1, 2] = 3 + self._assert_expected(cube, target, land_binary_mask) + + def test_missing_bottom_left_edge(self): + data = np.ma.array([[2, 4, 2], [4, 3, 6], [-1, 5, 2]], dtype=float) + cube = self.get_cube(data) + land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) + + moore_neighbourhood_search(cube, land_binary_mask, 10) + + target = data.data.copy() + target[2, 0] = 4 + self._assert_expected(cube, target, land_binary_mask) + + def test_missing_bottom_centre_edge(self): + data = np.ma.array([[2, 4, 2], [4, 3, 6], [5, -1, 2]], dtype=float) + cube = self.get_cube(data) + land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) + + moore_neighbourhood_search(cube, land_binary_mask, 10) + + target = data.data.copy() + target[2, 1] = 4 + self._assert_expected(cube, target, land_binary_mask) + + def test_missing_bottom_right_edge(self): + data = np.ma.array([[2, 4, 2], [4, 3, 6], [5, 3, -1]], dtype=float) + cube = self.get_cube(data) + land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) + + moore_neighbourhood_search(cube, land_binary_mask, 10) + + target = data.data.copy() + target[2, 2] = 4 + self._assert_expected(cube, target, land_binary_mask) + + def test_missing_including_zero_values(self): + # Ensure that the value 0 is included in the count for valid + # neighbours. + data = np.ma.array([[2, 4, 2], [4, 0, 7], [5, 3, -1]], dtype=float) + cube = self.get_cube(data) + land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) + + moore_neighbourhood_search(cube, land_binary_mask, 10) + + target = data.data.copy() + target[2, 2] = 10 / 3.0 + self._assert_expected(cube, target, land_binary_mask) + + def test_global_circular(self): + # Ensure that we correctly handle wraparound coordinates + data = np.ma.array([[2, 4, 2], [4, 3, 5], [5, 3, -1]], dtype=float) + cube = self.get_cube(data) + cube.coord(axis="x").circular = True + land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) + + moore_neighbourhood_search(cube, land_binary_mask, 10) + + target = data.data.copy() + target[2, 2] = 4 + self._assert_expected(cube, target, land_binary_mask) + + def test_global_circular_transposed(self): + # Ensure that we correctly handle wraparound coordinates where x and y + # are mapped to alternative dimensions. + data = np.ma.array([[2, 4, 2], [4, 3, 5], [5, 6, -1]], dtype=float) + cube = self.get_cube(data) + coord_x = cube.coord(axis="x") + coord_x.circular = True + coord_y = cube.coord(axis="y") + cube.remove_coord(coord_x) + cube.remove_coord(coord_y) + cube.add_dim_coord(coord_x, 0) + cube.add_dim_coord(coord_y, 1) + land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) + + moore_neighbourhood_search(cube, land_binary_mask, 10) + + target = data.data.copy() + target[2, 2] = 4 + self._assert_expected(cube, target, land_binary_mask) + + +class Test_missing_islands(_MaskApplication, ants.tests.TestCase): + def test_isolated_data(self): + data = np.ma.array([[-1, -1, -1], [-1, -1, -1], [-1, -1, -1]], dtype=float) + cube = self.get_cube(data) + land_binary_mask = cube.copy(np.zeros(data.shape, dtype=int)) + land_binary_mask.data[1, 1] = 1 + + moore_neighbourhood_search(cube, land_binary_mask, value=10) + + target = data.data.copy() + target[1, 1] = 10 + self._assert_expected(cube, target, land_binary_mask) + + def test_default_fill_value_no_valid_data(self): + # Ensure that islands stay as missing if we cannot infer a reasonable + # fill value (i.e. no valid data present to calculate a mean). + data = np.ma.array([[-1, -1, -1], [-1, -1, -1], [-1, -1, -1]], dtype=float) + cube = self.get_cube(data) + land_binary_mask = cube.copy(np.zeros(data.shape, dtype=int)) + land_binary_mask.data[1, 1] = 1 + + moore_neighbourhood_search(cube, land_binary_mask) + + target = data.data.copy() + target_land_binary_mask = cube.copy(np.zeros(data.shape, dtype=int)) + self._assert_expected(cube, target, target_land_binary_mask) + + def test_default_fill_value_valid_data(self): + # Ensure that the default fill value with fill the island with the data + # mean if there is valid data present. + data = np.ma.array([[-1, -1, 5], [-1, -1, 4], [1, 2, 3]], dtype=float) + cube = self.get_cube(data) + land_binary_mask = cube.copy() + land_binary_mask.data = (land_binary_mask.data > 0).astype(int) + land_binary_mask.data[0, 0] = 1 + + moore_neighbourhood_search(cube, land_binary_mask) + + target = data.data.copy() + target[0, 0] = 3 + self._assert_expected(cube, target, land_binary_mask) + + +class Test_mask_prioritisation(_MaskApplication, ants.tests.TestCase): + def test_masked_source_not_masked_target(self): + # Ensure that we use all source data available to derive missing points + # but modify the mask to remove these points that are masked in the + # provided binary land mask. + data = np.ma.array([[1, 2, 3], [4, 5, 5], [7, 8, -1]], dtype=float) + cube = self.get_cube(data) + land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) + land_binary_mask.data[1, 1] = 0 + + moore_neighbourhood_search(cube, land_binary_mask) + + target = data.copy() + target[2, 2] = 6.5 + self._assert_expected(cube, target, land_binary_mask) diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index f2e73c2..c6f96c2 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -120,12 +120,13 @@ class Test__transform_coordinates(ants.tests.TestCase): def setUp(self): self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) self.sphere_source = CubeBuilder(self.sphere_crs, (2, 2))._cube + self.points = np.array([[10, 10], [10, -10], [-10, -10], [-10, 10]]) def test_identity_rotation_sphere(self): """Test that rotation to a pole at latitude=90.0, longitude=0.0 returns the same points. - By convention, rotated pole co-ordinate systems will set the prime + By convention, rotated pole coordinate systems will set the prime meridian rotated 180.0 from the specified longitude. To place the prime meridian at 0.0, we apply a further rotation of 180.0, following rotation to the new pole. @@ -135,12 +136,12 @@ def test_identity_rotation_sphere(self): 90.0, 0.0, 180.0, ellipsoid=self.sphere_crs ) target_lsm = CubeBuilder(crs, (2, 2))._cube + expected_points = np.copy(self.points) - points = np.array([[10, 10], [10, -10], [-10, -10], [-10, 10]]) - - expected_points = np.copy(points) + rotated_coords = _transform_coordinates( + target_lsm, self.sphere_source, self.points + ) - rotated_coords = _transform_coordinates(target_lsm, self.sphere_source, points) self.assertArrayAlmostEqual(expected_points, rotated_coords) def test_longitudinal_rotation(self): @@ -150,25 +151,25 @@ def test_longitudinal_rotation(self): 90.0, 90.0, 180.0, ellipsoid=self.sphere_crs ) target_lsm = CubeBuilder(crs, (2, 2))._cube - - points = np.array([[10, 10], [10, -10], [-10, -10], [-10, 10]]) expected_points = np.array( [[-80.0, 10.0], [-80.0, -10.0], [-100.0, -10.0], [-100.0, 10]] ) - rotated_coords = _transform_coordinates(target_lsm, self.sphere_source, points) + rotated_coords = _transform_coordinates( + target_lsm, self.sphere_source, self.points + ) self.assertArrayAlmostEqual(expected_points, rotated_coords) def test_antimeridian_rotation(self): - """Test rotation to a new pole at latitude=90.0, longitude=180.0.""" + """Test rotation to a new pole at latitude=90.0, longitude=180.0. + In this case, the points are assumed to cross the anti meridian + and so are instead mapped to the interval [0, 360.0].""" crs = iris.coord_systems.RotatedGeogCS( 90.0, 180.0, 180.0, ellipsoid=self.sphere_crs ) target_lsm = CubeBuilder(crs, (2, 2))._cube - - points = np.array([[10, 10], [10, -10], [-10, -10], [-10, 10]]) expected_points = np.array( [[190.0, 10.0], [190.0, -10.0], [170.0, -10.0], [170.0, 10]] ) @@ -180,7 +181,43 @@ def test_antimeridian_rotation(self): ) with self.assertWarnsRegex(UserWarning, warning_msg): rotated_coords = _transform_coordinates( - target_lsm, self.sphere_source, points + target_lsm, self.sphere_source, self.points + ) + + self.assertArrayAlmostEqual(expected_points, rotated_coords) + + def test_local_target_no_domain_warning(self): + """Test that no warning is emitted if the transformed polygon + lies inside the target domain.""" + + crs = iris.coord_systems.RotatedGeogCS( + 90.0, 0.0, 170.0, ellipsoid=self.sphere_crs + ) + target_lsm = CubeBuilder(crs, (2, 2), xlim=(-25.0, 5.0), ylim=(-12, 12))._cube + expected_points = np.array([[0, 10], [0, -10], [-20, -10.0], [-20, 10]]) + + rotated_coords = _transform_coordinates( + target_lsm, self.sphere_source, self.points + ) + + self.assertArrayAlmostEqual(expected_points, rotated_coords) + + def test_local_target_domain_warning(self): + """Test that a warning is emitted if the transformed polygon + lies at least partially outside of the valid target domain.""" + + crs = iris.coord_systems.RotatedGeogCS( + 90.0, 0.0, 170.0, ellipsoid=self.sphere_crs + ) + target_lsm = CubeBuilder(crs, (2, 2), xlim=(-10.0, 5.0), ylim=(-12, 12))._cube + expected_points = np.array([[0, 10], [0, -10], [-20, -10.0], [-20, 10]]) + + warning_msg = re.escape( + "The transformed points lie outside the target lsm domain." + ) + with self.assertWarnsRegex(UserWarning, warning_msg): + rotated_coords = _transform_coordinates( + target_lsm, self.sphere_source, self.points ) self.assertArrayAlmostEqual(expected_points, rotated_coords) @@ -198,6 +235,7 @@ def test_error_raised(self): source="source/path", ) error_msg = "If --source is passed then --target-lsm must also be given." + with self.assertRaisesRegex(ValueError, error_msg): _validate_args(args.target_lsm, args.source) @@ -267,6 +305,7 @@ def test_transform_called(self): _ = _transform_if_required(target_lsm, self.sphere_source, self.points) received_lsm, received_source, received_points = mock_transform.call_args.args + mock_transform.assert_called_once() self.assertEqual(target_lsm, received_lsm) self.assertEqual(received_source, self.sphere_source) @@ -277,13 +316,17 @@ class Test_ite_transform(ants.tests.TestCase): def test_uk_pole_rotation(self): """ - Test the unrotated polygon enclosing regions of valid data from ITE - expressed in true longitude and latitude is rotated correctly. The target - pole is lon=177.5, lat=37.5. + The following test uses points which create a validity polygon + over the UK. The input points describe the points in an unrotated + geodetic coordinate system as an array of (lon, lat) pairs. The + points are then rotated to the target pole at lon=177.5, lat=37.5 The points have been obtained by unrotating the coordinates in - $UMDIR/ancil/data/shapefiles/ite_ukv_polygon/runme.py + $UMDIR/ancil/data/shapefiles/ite_ukv_polygon/runme.py, which can + also be found as part of the ANTS rose-stem test suite in + /data/users/ants/sources/ANTS/developer/core/ancil_create_shapefile/ """ + points = np.array( [ [1.63160953, 51.09703216], @@ -304,7 +347,6 @@ def test_uk_pole_rotation(self): [1.7594481, 51.23560201], ] ) - expected_points = np.array( [ [362.594, -1.32876], @@ -325,6 +367,7 @@ def test_uk_pole_rotation(self): [362.666, -1.18577], ] ) + sphere_crs = iris.coord_systems.GeogCS(6371229.0) source = CubeBuilder(sphere_crs, (2, 2))._cube crs = iris.coord_systems.RotatedGeogCS(37.5, 177.5, ellipsoid=sphere_crs) @@ -334,5 +377,4 @@ def test_uk_pole_rotation(self): # Adding 360.0 to longitude to compare to the expected points rotated_points[:, 0] += 360.0 - self.assertArrayAlmostEqual(rotated_points, expected_points) From 613372e80b849d116b04bf07c0e035f3597d45a8 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Thu, 27 Aug 2026 09:32:23 +0100 Subject: [PATCH 27/29] Added polygon domain warning and added json metadata file. --- lib/ants/cli/ancil_create_shapefile.py | 32 +- ...test_moore_neighbourhood_search.py.isorted | 281 ------------------ .../tests/cli/test_ancil_create_shapefile.py | 13 +- 3 files changed, 35 insertions(+), 291 deletions(-) delete mode 100644 lib/ants/tests/analysis/merge/test_moore_neighbourhood_search.py.isorted diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index 89fbef6..9585f26 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -19,12 +19,13 @@ import argparse import json import logging +import os import warnings import ants import iris.coord_systems import numpy as np -from ants.io.load import load_cube, load_landsea_mask +from ants.io.load import load_cube from ants.utils.cube import CubeBuilder from osgeo import ogr from shapely.geometry import Polygon @@ -237,7 +238,7 @@ def _load_cubes(target_lsm_path, source_path): A tuple containing the target lsm and the source cube respectively. """ - target_lsm = load_landsea_mask(target_lsm_path) + target_lsm = load_cube(target_lsm_path) if source_path is None: crs = iris.coord_systems.GeogCS(6371229.0) @@ -283,6 +284,32 @@ def _transform_if_required(target_lsm, source, points): return rotated_points +def _save_json(output, target_lsm_path, source_path, target_lsm, source): + """ + Save a json file containing metadata about the target lsm and source. + """ + + parent_path = os.path.dirname(output) + target_coord = target_lsm.coord_system() + target_proj_params = target_coord.as_cartopy_crs().proj4_params + source_proj_params = source.coord_system().as_cartopy_crs().proj4_params + + prj_metadata = { + "target_lsm": { + "target_lsm_path": target_lsm_path, + "proj4_params": target_proj_params, + "grid_north_pole_longitude": target_coord.grid_north_pole_longitude, + "grid_north_pole_latitude": target_coord.grid_north_pole_latitude, + "north_pole_grid_longitude": target_coord.north_pole_grid_longitude, + }, + "source": {"source_path": source_path, "proj4_params": source_proj_params}, + } + json_name = "prj_metadata" + + with open(os.path.join(parent_path, json_name + ".json"), "w") as json_file: + json.dump(prj_metadata, json_file, indent=4) + + def _load_points_from_json(json_file): """ Load a json file containing a list of pairs of longitude-latitude points. @@ -337,6 +364,7 @@ def main(json_file, output, target_lsm_path, source_path): if target_lsm_path is not None: target_lsm, source = _load_cubes(target_lsm_path, source_path) points = _transform_if_required(target_lsm, source, points) + _save_json(output, target_lsm_path, source_path, target_lsm, source) polygon = Polygon(points) diff --git a/lib/ants/tests/analysis/merge/test_moore_neighbourhood_search.py.isorted b/lib/ants/tests/analysis/merge/test_moore_neighbourhood_search.py.isorted deleted file mode 100644 index 1a5e120..0000000 --- a/lib/ants/tests/analysis/merge/test_moore_neighbourhood_search.py.isorted +++ /dev/null @@ -1,281 +0,0 @@ -# (C) Crown Copyright, Met Office. All rights reserved. -# -# This file is part of ANTS and is released under the BSD 3-Clause license. -# See LICENSE.txt in the root of the repository for full licensing details. -import ants.tests -import iris -import numpy as np -from ants.analysis._merge import moore_neighbourhood_search - - -class _MaskApplication(object): - def get_cube(self, data): - # Return a global cube with -1 values the assumed masked values - cube = ants.tests.stock.geodetic(data.shape, with_bounds=True) - cube.coord(axis="x").circular = False - data = np.ma.masked_values(data, -1, copy=False) - cube.data = data - return cube - - def _assert_expected(self, cube, target, land_binary_mask): - self.assertArrayEqual(cube.data.data, target) - self.assertArrayEqual(cube.data.mask, land_binary_mask.data == 0) - - -class Test_inappropriate_input(_MaskApplication, ants.tests.TestCase): - def test_not_2d_source(self): - cube = iris.cube.Cube(np.arange(8).reshape(2, 2, 2)) - land_binary_mask = iris.cube.Cube(np.arange(8).reshape(2, 2, 2)) - msg = ( - "Currently, only mask application to 2D grids are supported " - "and with no broadcasting." - ) - with self.assertRaisesRegex(RuntimeError, msg): - moore_neighbourhood_search(cube, land_binary_mask) - - def test_different_grids(self): - data = np.arange(4).reshape(2, 2) - cube = self.get_cube(data) - land_binary_mask = self.get_cube(data) - land_binary_mask.coord(axis="x").points = cube.coord(axis="x").points - 10 - msg = "Both source and land_binary_mask must be defined on identical " "grids" - with self.assertRaisesRegex(RuntimeError, msg): - moore_neighbourhood_search(cube, land_binary_mask) - - def test_alternate_orientation(self): - # Currently we do not support automatically transposing the data where - # coordinates are mapped to alternate dimensions. Make sure we catch - # this case. - data = np.arange(4).reshape(2, 2) - cube = self.get_cube(data) - land_binary_mask = self.get_cube(data) - lx = land_binary_mask.coord(axis="x") - ly = land_binary_mask.coord(axis="y") - land_binary_mask.remove_coord(lx) - land_binary_mask.remove_coord(ly) - land_binary_mask.add_dim_coord(lx, 0) - land_binary_mask.add_dim_coord(ly, 1) - - msg = ( - "Currently, the source and the land_binary_mask must be " - "defined in the same orientation" - ) - with self.assertRaisesRegex(RuntimeError, msg): - moore_neighbourhood_search(cube, land_binary_mask) - - -class Test_missing_coastlines(_MaskApplication, ants.tests.TestCase): - def test_no_change(self): - data = np.ma.array([[-1, 2, 3], [4, -1, -1], [7, 8, 9]], dtype=float) - cube = self.get_cube(data) - land_binary_mask = cube.copy((data > 0).astype(int)) - - moore_neighbourhood_search(cube, land_binary_mask, 10) - - target = data.copy() - self._assert_expected(cube, target, land_binary_mask) - - def test_missing_top_left_edge(self): - data = np.ma.array([[-1, 2, 3], [2, 5, 6], [7, 8, 9]], dtype=float) - cube = self.get_cube(data) - land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) - - moore_neighbourhood_search(cube, land_binary_mask, 10) - - target = data.data.copy() - target[0, 0] = 3 - self._assert_expected(cube, target, land_binary_mask) - - def test_missing_top_centre_edge(self): - data = np.ma.array([[4, -1, 3], [2, 5, 6], [7, 8, 9]], dtype=float) - cube = self.get_cube(data) - land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) - - moore_neighbourhood_search(cube, land_binary_mask, 10) - - target = data.data.copy() - target[0, 1] = 4 - self._assert_expected(cube, target, land_binary_mask) - - def test_missing_top_right_edge(self): - data = np.ma.array([[4, 2, -1], [2, 5, 2], [7, 8, 9]], dtype=float) - cube = self.get_cube(data) - land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) - - moore_neighbourhood_search(cube, land_binary_mask, 10) - - target = data.data.copy() - target[0, 2] = 3 - self._assert_expected(cube, target, land_binary_mask) - - def test_missing_centre_left_edge(self): - data = np.ma.array([[4, 2, 2], [-1, 5, 2], [3, 6, 9]], dtype=float) - cube = self.get_cube(data) - land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) - - moore_neighbourhood_search(cube, land_binary_mask, 10) - - target = data.data.copy() - target[1, 0] = 4 - self._assert_expected(cube, target, land_binary_mask) - - def test_missing_centre_centre(self): - data = np.ma.array([[2, 4, 2], [4, -1, 4], [2, 4, 2]], dtype=float) - cube = self.get_cube(data) - land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) - - moore_neighbourhood_search(cube, land_binary_mask, 10) - - target = data.data.copy() - target[1, 1] = 3 - self._assert_expected(cube, target, land_binary_mask) - - def test_missing_centre_right_edge(self): - data = np.ma.array([[2, 4, 2], [4, 3, -1], [2, 4, 2]], dtype=float) - cube = self.get_cube(data) - land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) - - moore_neighbourhood_search(cube, land_binary_mask, 10) - - target = data.data.copy() - target[1, 2] = 3 - self._assert_expected(cube, target, land_binary_mask) - - def test_missing_bottom_left_edge(self): - data = np.ma.array([[2, 4, 2], [4, 3, 6], [-1, 5, 2]], dtype=float) - cube = self.get_cube(data) - land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) - - moore_neighbourhood_search(cube, land_binary_mask, 10) - - target = data.data.copy() - target[2, 0] = 4 - self._assert_expected(cube, target, land_binary_mask) - - def test_missing_bottom_centre_edge(self): - data = np.ma.array([[2, 4, 2], [4, 3, 6], [5, -1, 2]], dtype=float) - cube = self.get_cube(data) - land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) - - moore_neighbourhood_search(cube, land_binary_mask, 10) - - target = data.data.copy() - target[2, 1] = 4 - self._assert_expected(cube, target, land_binary_mask) - - def test_missing_bottom_right_edge(self): - data = np.ma.array([[2, 4, 2], [4, 3, 6], [5, 3, -1]], dtype=float) - cube = self.get_cube(data) - land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) - - moore_neighbourhood_search(cube, land_binary_mask, 10) - - target = data.data.copy() - target[2, 2] = 4 - self._assert_expected(cube, target, land_binary_mask) - - def test_missing_including_zero_values(self): - # Ensure that the value 0 is included in the count for valid - # neighbours. - data = np.ma.array([[2, 4, 2], [4, 0, 7], [5, 3, -1]], dtype=float) - cube = self.get_cube(data) - land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) - - moore_neighbourhood_search(cube, land_binary_mask, 10) - - target = data.data.copy() - target[2, 2] = 10 / 3.0 - self._assert_expected(cube, target, land_binary_mask) - - def test_global_circular(self): - # Ensure that we correctly handle wraparound coordinates - data = np.ma.array([[2, 4, 2], [4, 3, 5], [5, 3, -1]], dtype=float) - cube = self.get_cube(data) - cube.coord(axis="x").circular = True - land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) - - moore_neighbourhood_search(cube, land_binary_mask, 10) - - target = data.data.copy() - target[2, 2] = 4 - self._assert_expected(cube, target, land_binary_mask) - - def test_global_circular_transposed(self): - # Ensure that we correctly handle wraparound coordinates where x and y - # are mapped to alternative dimensions. - data = np.ma.array([[2, 4, 2], [4, 3, 5], [5, 6, -1]], dtype=float) - cube = self.get_cube(data) - coord_x = cube.coord(axis="x") - coord_x.circular = True - coord_y = cube.coord(axis="y") - cube.remove_coord(coord_x) - cube.remove_coord(coord_y) - cube.add_dim_coord(coord_x, 0) - cube.add_dim_coord(coord_y, 1) - land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) - - moore_neighbourhood_search(cube, land_binary_mask, 10) - - target = data.data.copy() - target[2, 2] = 4 - self._assert_expected(cube, target, land_binary_mask) - - -class Test_missing_islands(_MaskApplication, ants.tests.TestCase): - def test_isolated_data(self): - data = np.ma.array([[-1, -1, -1], [-1, -1, -1], [-1, -1, -1]], dtype=float) - cube = self.get_cube(data) - land_binary_mask = cube.copy(np.zeros(data.shape, dtype=int)) - land_binary_mask.data[1, 1] = 1 - - moore_neighbourhood_search(cube, land_binary_mask, value=10) - - target = data.data.copy() - target[1, 1] = 10 - self._assert_expected(cube, target, land_binary_mask) - - def test_default_fill_value_no_valid_data(self): - # Ensure that islands stay as missing if we cannot infer a reasonable - # fill value (i.e. no valid data present to calculate a mean). - data = np.ma.array([[-1, -1, -1], [-1, -1, -1], [-1, -1, -1]], dtype=float) - cube = self.get_cube(data) - land_binary_mask = cube.copy(np.zeros(data.shape, dtype=int)) - land_binary_mask.data[1, 1] = 1 - - moore_neighbourhood_search(cube, land_binary_mask) - - target = data.data.copy() - target_land_binary_mask = cube.copy(np.zeros(data.shape, dtype=int)) - self._assert_expected(cube, target, target_land_binary_mask) - - def test_default_fill_value_valid_data(self): - # Ensure that the default fill value with fill the island with the data - # mean if there is valid data present. - data = np.ma.array([[-1, -1, 5], [-1, -1, 4], [1, 2, 3]], dtype=float) - cube = self.get_cube(data) - land_binary_mask = cube.copy() - land_binary_mask.data = (land_binary_mask.data > 0).astype(int) - land_binary_mask.data[0, 0] = 1 - - moore_neighbourhood_search(cube, land_binary_mask) - - target = data.data.copy() - target[0, 0] = 3 - self._assert_expected(cube, target, land_binary_mask) - - -class Test_mask_prioritisation(_MaskApplication, ants.tests.TestCase): - def test_masked_source_not_masked_target(self): - # Ensure that we use all source data available to derive missing points - # but modify the mask to remove these points that are masked in the - # provided binary land mask. - data = np.ma.array([[1, 2, 3], [4, 5, 5], [7, 8, -1]], dtype=float) - cube = self.get_cube(data) - land_binary_mask = cube.copy(np.ones(data.shape, dtype=int)) - land_binary_mask.data[1, 1] = 0 - - moore_neighbourhood_search(cube, land_binary_mask) - - target = data.copy() - target[2, 2] = 6.5 - self._assert_expected(cube, target, land_binary_mask) diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index c6f96c2..3d216b0 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -244,14 +244,12 @@ class Test__load_cubes(ants.tests.TestCase): def test_load_lsm_called(self): """Test that loading the landsea mask is successfully called.""" - with mock.patch( - "ants.cli.ancil_create_shapefile.load_landsea_mask" - ) as mock_lsm: + with mock.patch("ants.cli.ancil_create_shapefile.load_cube") as mock_load: _ = _load_cubes("target/path", None) - mock_lsm.assert_called_once_with("target/path") + mock_load.assert_called_once_with("target/path") - @mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") + @mock.patch("ants.cli.ancil_create_shapefile.load_cube") def test_cubebuilder_called(self, *args): """Test source cube created if no path given.""" @@ -262,14 +260,13 @@ def test_cubebuilder_called(self, *args): mock_build_cube.assert_called_once() - @mock.patch("ants.cli.ancil_create_shapefile.load_landsea_mask") def test_load_cube_called(self, *args): """Test load cube called if source path given.""" - with mock.patch("ants.cli.ancil_create_shapefile.load_cube") as mock_load_cube: + with mock.patch("ants.cli.ancil_create_shapefile.load_cube") as mock_load: _ = _load_cubes("target/path", "source/cube") - mock_load_cube.assert_called_once() + self.assertEqual(2, mock_load.call_count) class Test__transform_if_required(ants.tests.TestCase): From c0cc24fcc0be3890add0382b8589be192c9667f3 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Wed, 2 Sep 2026 10:58:31 +0100 Subject: [PATCH 28/29] Minor format changes. --- lib/ants/cli/ancil_create_shapefile.py | 28 +++++++++++------- .../tests/cli/test_ancil_create_shapefile.py | 29 ++++++++++++------- 2 files changed, 37 insertions(+), 20 deletions(-) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index 9585f26..bb8ba96 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -290,21 +290,29 @@ def _save_json(output, target_lsm_path, source_path, target_lsm, source): """ parent_path = os.path.dirname(output) - target_coord = target_lsm.coord_system() - target_proj_params = target_coord.as_cartopy_crs().proj4_params - source_proj_params = source.coord_system().as_cartopy_crs().proj4_params + filename = os.path.splitext(os.path.basename(output))[0] + + target_crs = target_lsm.coord_system() + target_proj_params = target_crs.as_cartopy_crs().proj4_params + source_crs = source.coord_system() + source_proj_params = source_crs.as_cartopy_crs().proj4_params prj_metadata = { "target_lsm": { "target_lsm_path": target_lsm_path, + "coord_system": type(target_crs).__name__, + "grid_north_pole_longitude": target_crs.grid_north_pole_longitude, + "grid_north_pole_latitude": target_crs.grid_north_pole_latitude, + "north_pole_grid_longitude": target_crs.north_pole_grid_longitude, "proj4_params": target_proj_params, - "grid_north_pole_longitude": target_coord.grid_north_pole_longitude, - "grid_north_pole_latitude": target_coord.grid_north_pole_latitude, - "north_pole_grid_longitude": target_coord.north_pole_grid_longitude, }, - "source": {"source_path": source_path, "proj4_params": source_proj_params}, + "source": { + "source_path": source_path, + "coord_system": type(source_crs).__name__, + "proj4_params": source_proj_params, + }, } - json_name = "prj_metadata" + json_name = filename + "_prj" with open(os.path.join(parent_path, json_name + ".json"), "w") as json_file: json.dump(prj_metadata, json_file, indent=4) @@ -361,13 +369,13 @@ def main(json_file, output, target_lsm_path, source_path): points = _load_points_from_json(json_file) ccw_expected = Polygon(points).exterior.is_ccw + # Transform points to a rotated pole if required if target_lsm_path is not None: target_lsm, source = _load_cubes(target_lsm_path, source_path) points = _transform_if_required(target_lsm, source, points) _save_json(output, target_lsm_path, source_path, target_lsm, source) polygon = Polygon(points) - _check_polygon_validity(polygon, ccw_expected) # Now convert it to a shapefile with OGR @@ -411,7 +419,7 @@ def _get_parser(): "--source", type=ants.config.filepath_readable, required=False, - help="Path to a source containing the coordinate system of the json file.", + help="Path to a source specifying the coordinate system of the json file.", ) return parser diff --git a/lib/ants/tests/cli/test_ancil_create_shapefile.py b/lib/ants/tests/cli/test_ancil_create_shapefile.py index 3d216b0..e339307 100644 --- a/lib/ants/tests/cli/test_ancil_create_shapefile.py +++ b/lib/ants/tests/cli/test_ancil_create_shapefile.py @@ -191,7 +191,7 @@ def test_local_target_no_domain_warning(self): lies inside the target domain.""" crs = iris.coord_systems.RotatedGeogCS( - 90.0, 0.0, 170.0, ellipsoid=self.sphere_crs + 90.0, 10.0, 180.0, ellipsoid=self.sphere_crs ) target_lsm = CubeBuilder(crs, (2, 2), xlim=(-25.0, 5.0), ylim=(-12, 12))._cube expected_points = np.array([[0, 10], [0, -10], [-20, -10.0], [-20, 10]]) @@ -207,7 +207,7 @@ def test_local_target_domain_warning(self): lies at least partially outside of the valid target domain.""" crs = iris.coord_systems.RotatedGeogCS( - 90.0, 0.0, 170.0, ellipsoid=self.sphere_crs + 90.0, 10.0, 180.0, ellipsoid=self.sphere_crs ) target_lsm = CubeBuilder(crs, (2, 2), xlim=(-10.0, 5.0), ylim=(-12, 12))._cube expected_points = np.array([[0, 10], [0, -10], [-20, -10.0], [-20, 10]]) @@ -271,7 +271,7 @@ def test_load_cube_called(self, *args): class Test__transform_if_required(ants.tests.TestCase): def setUp(self): - self.points = np.array([[1, 2], [3, 4], [5, 6], [7, 8]]) + self.points = np.array([[10, 10], [10, -10], [-10, -10], [-10, 10]]) self.sphere_crs = iris.coord_systems.GeogCS(6371229.0) self.sphere_source = CubeBuilder(self.sphere_crs, (2, 2))._cube @@ -313,15 +313,24 @@ class Test_ite_transform(ants.tests.TestCase): def test_uk_pole_rotation(self): """ - The following test uses points which create a validity polygon - over the UK. The input points describe the points in an unrotated - geodetic coordinate system as an array of (lon, lat) pairs. The - points are then rotated to the target pole at lon=177.5, lat=37.5 + The following test uses points that define a validity polygon over the UK. + + The input points are specified in an unrotated geodetic coordinate system + as (lon, lat) pairs. They are then rotated to a pole at lon=177.5, + lat=37.5. Note that PROJ applies an additional 180.0 degree rotation + to the specified pole location. - The points have been obtained by unrotating the coordinates in + The points have been obtained by unrotating the coordinates under $UMDIR/ancil/data/shapefiles/ite_ukv_polygon/runme.py, which can - also be found as part of the ANTS rose-stem test suite in + also be found as part of the ANTS rose-stem test suite under /data/users/ants/sources/ANTS/developer/core/ancil_create_shapefile/ + + Notes + ----- + The transformed points differ from the expected points by 360 degrees in + longitude. Cartopy treats these longitudes as equivalent, but operations + that perform Cartesian comparisons may not consider the resulting polygons + equivalent. """ points = np.array( @@ -372,6 +381,6 @@ def test_uk_pole_rotation(self): rotated_points = _transform_if_required(target_lsm, source, points) - # Adding 360.0 to longitude to compare to the expected points + # Adding 360.0 to longitude to compare to the expected points. rotated_points[:, 0] += 360.0 self.assertArrayAlmostEqual(rotated_points, expected_points) From 4b68e1847fecc9c47fa6281bd1a30a72469e0e22 Mon Sep 17 00:00:00 2001 From: Alasdair Roy <292044538+alasdair-roy@users.noreply.github.com> Date: Wed, 2 Sep 2026 11:25:33 +0100 Subject: [PATCH 29/29] Added comment on polygon orientation. --- lib/ants/cli/ancil_create_shapefile.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/lib/ants/cli/ancil_create_shapefile.py b/lib/ants/cli/ancil_create_shapefile.py index bb8ba96..3cd817a 100755 --- a/lib/ants/cli/ancil_create_shapefile.py +++ b/lib/ants/cli/ancil_create_shapefile.py @@ -65,6 +65,9 @@ def _check_polygon_validity(polygon, ccw_expected): """ Check if polygon is valid and has the expected orientation. + A rotation should not change the orientation of the points. If the orientation + does change, it is a sign that the transformation has not behaved as expected. + Parameters ---------- polygon : :class:`~shapely.geometry.Polygon`