diff --git a/cf/__init__.py b/cf/__init__.py index bbde8d4002..8c9a767518 100644 --- a/cf/__init__.py +++ b/cf/__init__.py @@ -254,6 +254,9 @@ from .domainaxis import DomainAxis from .fieldancillary import FieldAncillary from .field import Field + +from .gridmappings import * # noqa: F403 + from .data import Data from .data.array import ( CFANetCDFArray, diff --git a/cf/constants.py b/cf/constants.py index ea335f63ee..93b86dee59 100644 --- a/cf/constants.py +++ b/cf/constants.py @@ -72,26 +72,123 @@ _stash2standard_name = {} + # --------------------------------------------------------------------- -# Coordinate reference constants TODO: turn these into functions +# Coordinate reference and grid mapping related constants # --------------------------------------------------------------------- +# Grid Mapping valid attribute names. Taken from the list given in +# 'Table F.1. Grid Mapping Attributes' in Appendix F: Grid Mappings' +# of the Conventions document. +# +# Values indicate if attribute values are expected to be numeric (instead +# of string) for the given attribute key (the table defines this via N, S). +# For any value where this is True, canonical units must be provided +# in the 'cr_canonical_units' dict below. +cr_gm_valid_attr_names_are_numeric = { + "grid_mapping_name": False, + # *Those which describe the ellipsoid and prime meridian:* + "earth_radius": True, + "inverse_flattening": True, + "longitude_of_prime_meridian": True, + "prime_meridian_name": False, + "reference_ellipsoid_name": False, + "semi_major_axis": True, + "semi_minor_axis": True, + # *Specific/applicable to only given grid mapping(s):* + # ...projection origin related: + "longitude_of_projection_origin": True, + "latitude_of_projection_origin": True, + "scale_factor_at_projection_origin": True, + # ...false-Xings: + "false_easting": True, + "false_northing": True, + # ...angle axis related: + "sweep_angle_axis": False, + "fixed_angle_axis": False, + # ...central meridian related: + "longitude_of_central_meridian": True, + "scale_factor_at_central_meridian": True, + # ...pole coordinates related: + "grid_north_pole_latitude": True, + "grid_north_pole_longitude": True, + "north_pole_grid_longitude": True, + # ...other: + "standard_parallel": True, + "perspective_point_height": True, + "azimuth_of_central_line": True, + "straight_vertical_longitude_from_pole": True, + # *Other, not needed for a specific grid mapping but also listed + # in 'Table F.1. Grid Mapping Attributes':* + "crs_wkt": False, + "geographic_crs_name": False, + "geoid_name": False, + "geopotential_datum_name": False, + "horizontal_datum_name": False, + "projected_crs_name": False, + "towgs84": True, +} cr_canonical_units = { + # *Those which describe the ellipsoid and prime meridian:* "earth_radius": Units("m"), - "grid_north_pole_latitude": Units("degrees_north"), - "grid_north_pole_longitude": Units("degrees_east"), "inverse_flattening": Units("1"), - "latitude_of_projection_origin": Units("degrees_north"), - "longitude_of_central_meridian": Units("degrees_east"), "longitude_of_prime_meridian": Units("degrees_east"), - "longitude_of_projection_origin": Units("degrees_east"), - "north_pole_grid_longitude": Units("degrees"), - "perspective_point_height": Units("m"), - "scale_factor_at_central_meridian": Units("1"), - "scale_factor_at_projection_origin": Units("1"), "semi_major_axis": Units("m"), "semi_minor_axis": Units("m"), + # *Specific/applicable to only given grid mapping(s):* + # ...projection origin related: + "longitude_of_projection_origin": Units("degrees_east"), + "latitude_of_projection_origin": Units("degrees_north"), + "scale_factor_at_projection_origin": Units("1"), + # ...false-Xings: + "false_easting": Units("m"), + "false_northing": Units("m"), + # ...central meridian related: + "longitude_of_central_meridian": Units("degrees_east"), + "scale_factor_at_central_meridian": Units("1"), + # ...pole coordinates related: + "grid_north_pole_latitude": Units("degrees_north"), + "grid_north_pole_longitude": Units("degrees_east"), + "north_pole_grid_longitude": Units("degrees"), + # ...other: "standard_parallel": Units("degrees_north"), - "straight_vertical_longitude_from_pole": Units("degrees_north"), + "perspective_point_height": Units("m"), + "azimuth_of_central_line": Units("degrees"), + # TODOGM check the below is correct, was 'degrees_north' before but it + # is a longitude so surely that was wrong?: + "straight_vertical_longitude_from_pole": Units("degrees_east"), + # *Other, not needed for a specific grid mapping but also listed + # in 'Table F.1. Grid Mapping Attributes':* + # TODOGM towgs84 is very complicated, with multiple components + # of translations, roations and scaling each with different units + # requirements, so we should probably sort this out later, see: + # https://proj.org/en/9.3/operations/transformations/ + # helmert.html#helmert-transform + # where translations are in meters, rotations in arc-seconds, and + # scaling is in parts per million. + "towgs84": None, +} + +# Indicates what PROJ string ID component(s) refer(s) to the grid mapping +# attribute in question when '+' is added as a suffix, for example +# 'earth_radius' is '+R' and 'inverse_flattening' is '+rf' +cr_gm_attr_to_proj_string_comps = { + "earth_radius": "R", + "inverse_flattening": "rf", + "prime_meridian_name": "pm", + "reference_ellipsoid_name": "ellps", + "semi_major_axis": "a", + "semi_minor_axis": "b", + "longitude_of_projection_origin": "lon_0", + "latitude_of_projection_origin": "lat_0", + "scale_factor_at_projection_origin": "k_0", + "false_easting": "x_0", + "false_northing": "y_0", + "sweep_angle_axis": "sweep", + "standard_parallel": ("lat_1", "lat_2"), + "perspective_point_height": "h", + "azimuth_of_central_line": ("alpha", "gamma"), + "crs_wkt": "crs_wkt", + "towgs84": "towgs84", } cr_coordinates = { diff --git a/cf/domain.py b/cf/domain.py index fe43779025..cf68c7a19a 100644 --- a/cf/domain.py +++ b/cf/domain.py @@ -18,6 +18,7 @@ indices_shape, parse_indices, ) +from .gridmappings import GM, InvalidGridMapping _empty_set = set() @@ -585,6 +586,59 @@ def get_data(self, default=ValueError(), _units=None, _fill_value=True): default, message=f"{self.__class__.__name__} has no data" ) + def get_grid_mappings(self, as_class=False): + """Returns coordinate conversions with their grid mappings. + + .. versionadded:: GMVER + + :Parameters: + + as_class: `bool`, optional + If `True`, return the grid mapping as the equivalent + CF Grid Mapping class, for example + cf.RotatedLatitudeLongitude, rather than as a string + corresponding to the value of the 'grid_mapping_name' + attribute, for example 'rotated_latitude_longitude'. + By default the 'grid_mapping_name' value is returned. + + :Returns: + + `dict` + CoordinateConversion construct identifiers with + values of their 'grid_mapping_name' attribute, + or corresponding CF Grid Mapping class if + as_class is `True`, for all CoordinateConversions + of the domain that have a 'grid_mapping_name' + parameter defined. + + **Examples** + + >>> f.get_grid_mappings() + {'coordinatereference1': "rotated_latitude_longitude"} + >>> f.get_grid_mappings(as_class=True) + {'coordinatereference1': cf.gridmappings.RotatedLatitudeLongitude} + >>> g.get_grid_mappings() + {} + + """ + gms = {} + for cref_name, cref in self.coordinate_references().items(): + # If there is no coordinate conversion or parameters set on one, + # this will give None, so is safe + gm = cref.coordinate_conversion.get_parameter( + "grid_mapping_name", default=None + ) + if gm: + if as_class: + try: + gm_class = GM(cref) + gms[cref_name] = gm_class + except InvalidGridMapping: + pass # not a supported GM so don't add + else: + gms[cref_name] = gm + return gms + def identity(self, default="", strict=False, relaxed=False, nc_only=False): """Return the canonical identity. diff --git a/cf/field.py b/cf/field.py index 0ec52aa34e..6fe8b05824 100644 --- a/cf/field.py +++ b/cf/field.py @@ -15667,3 +15667,40 @@ def unlimited(self, *args): version="3.0.0", removed_at="4.0.0", ) # pragma: no cover + + def get_grid_mappings(self, as_class=False): + """Returns coordinate conversions with their grid mappings. + + .. versionadded:: GMVER + + :Parameters: + + as_class: `bool`, optional + If `True`, return the grid mapping as the equivalent + CF Grid Mapping class, for example + cf.RotatedLatitudeLongitude, rather than as a string + corresponding to the value of the 'grid_mapping_name' + attribute, for example 'rotated_latitude_longitude'. + By default the 'grid_mapping_name' value is returned. + + :Returns: + + `dict` + CoordinateConversion construct identifiers with + values of their 'grid_mapping_name' attribute, + or corresponding CF Grid Mapping class if + as_class is `True`, for all CoordinateConversions + of the domain that have a 'grid_mapping_name' + parameter defined. + + **Examples** + + >>> f.get_grid_mappings() + {'coordinatereference1': "rotated_latitude_longitude"} + >>> f.get_grid_mappings(as_class=True) + {'coordinatereference1': cf.gridmappings.RotatedLatitudeLongitude} + >>> g.get_grid_mappings() + {} + + """ + return self.domain.get_grid_mappings(as_class=as_class) diff --git a/cf/gridmappings/__init__.py b/cf/gridmappings/__init__.py new file mode 100644 index 0000000000..1878ab2657 --- /dev/null +++ b/cf/gridmappings/__init__.py @@ -0,0 +1,40 @@ +""" +Module for Grid Mappings supported by the CF Conventions. + +For the full list of supported Grid Mappings and details, see Appendix F +of the canonical document: + +https://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ +cf-conventions.html#appendix-grid-mappings + +This module should be kept up to date with the Appendix, by adding or +amending appropriate classes. + +Note that abstract classes to support the creation of concrete classes +for concrete Grid Mappinds are defined in the 'abstract' module, so +not included in the listing below. + +""" + + +from .abstract import * +from .gridmapping import ( + GM, + InvalidGridMapping, + AlbersEqualArea, + AzimuthalEquidistant, + Geostationary, + LambertAzimuthalEqualArea, + LambertConformalConic, + LambertCylindricalEqualArea, + Mercator, + ObliqueMercator, + Orthographic, + PolarStereographic, + RotatedLatitudeLongitude, + LatitudeLongitude, + Sinusoidal, + Stereographic, + TransverseMercator, + VerticalPerspective, +) diff --git a/cf/gridmappings/abstract/__init__.py b/cf/gridmappings/abstract/__init__.py new file mode 100644 index 0000000000..0712560fe5 --- /dev/null +++ b/cf/gridmappings/abstract/__init__.py @@ -0,0 +1,25 @@ +""" +Abstract classes representing categories of Grid Mapping. + +Categories correspond to the type of projection a Grid Mapping is +based upon, for example these are often based upon a geometric shape +that forms the developable surface that is used to flatten the map +of Earth (such as conic for a cone or cylindrical for a cylinder). + +The LatLonGridMapping case is special in that it (from the CF Conventions, +Appendix F) 'defines the canonical 2D geographical coordinate system +based upon latitude and longitude coordinates on a spherical Earth'. + +""" + +from .gridmappingbase import ( + GridMapping, + convert_proj_angular_data_to_cf, + convert_cf_angular_data_to_proj, + _validate_map_parameter, +) +from .azimuthal import AzimuthalGridMapping +from .conic import ConicGridMapping +from .cylindrical import CylindricalGridMapping +from .latlon import LatLonGridMapping +from .perspective import PerspectiveGridMapping diff --git a/cf/gridmappings/abstract/azimuthal.py b/cf/gridmappings/abstract/azimuthal.py new file mode 100644 index 0000000000..e6f646d4c3 --- /dev/null +++ b/cf/gridmappings/abstract/azimuthal.py @@ -0,0 +1,65 @@ +from .gridmappingbase import ( + GridMapping, + _validate_map_parameter, +) + + +class AzimuthalGridMapping(GridMapping): + """A Grid Mapping with Azimuthal classification. + + .. versionadded:: GMVER + + :Parameters: + + longitude_of_projection_origin: number or scalar `Data`, optional + The longitude of projection center (PROJ 'lon_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + latitude_of_projection_origin: number or scalar `Data`, optional + The latitude of projection center (PROJ 'lat_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_north', else the `Data` units are + taken and must be angular and compatible with latitude. + The default is 0.0 degrees_north. + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + """ + + def __init__( + self, + longitude_of_projection_origin=0.0, + latitude_of_projection_origin=0.0, + false_easting=0.0, + false_northing=0.0, + **kwargs, + ): + super().__init__(**kwargs) + + self.longitude_of_projection_origin = _validate_map_parameter( + "longitude_of_projection_origin", longitude_of_projection_origin + ) + self.latitude_of_projection_origin = _validate_map_parameter( + "latitude_of_projection_origin", latitude_of_projection_origin + ) + self.false_easting = _validate_map_parameter( + "false_easting", false_easting + ) + self.false_northing = _validate_map_parameter( + "false_northing", false_northing + ) diff --git a/cf/gridmappings/abstract/conic.py b/cf/gridmappings/abstract/conic.py new file mode 100644 index 0000000000..0c1989d947 --- /dev/null +++ b/cf/gridmappings/abstract/conic.py @@ -0,0 +1,85 @@ +from .gridmappingbase import ( + GridMapping, + _validate_map_parameter, +) + + +class ConicGridMapping(GridMapping): + """A Grid Mapping with Conic classification. + + .. versionadded:: GMVER + + :Parameters: + + standard_parallel: 2-`tuple` of number or scalar `Data` or `None` + The standard parallel value(s): the first (PROJ 'lat_1' + value) and/or the second (PROJ 'lat_2' value), given + as a 2-tuple of numbers or strings corresponding to + the first and then the second in order, where `None` + indicates that a value is not being set for either. + + If provided as a number or `Data` without units, the units + for each of the values are taken as 'degrees_north', else + the `Data` units are taken and must be angular and + compatible with latitude. + + The default is (0.0, 0.0), that is 0.0 degrees_north + for the first and second standard parallel values. + + longitude_of_central_meridian: number or scalar `Data`, optional + The longitude of (natural) origin i.e. central meridian. + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + latitude_of_projection_origin: number or scalar `Data`, optional + The latitude of projection center (PROJ 'lat_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_north', else the `Data` units are + taken and must be angular and compatible with latitude. + The default is 0.0 degrees_north. + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + """ + + def __init__( + self, + standard_parallel, + longitude_of_central_meridian=0.0, + latitude_of_projection_origin=0.0, + false_easting=0.0, + false_northing=0.0, + **kwargs, + ): + super().__init__(**kwargs) + + self.standard_parallel = ( + _validate_map_parameter("standard_parallel", standard_parallel[0]), + _validate_map_parameter("standard_parallel", standard_parallel[1]), + ) + self.longitude_of_central_meridian = _validate_map_parameter( + "longitude_of_central_meridian", longitude_of_central_meridian + ) + self.latitude_of_projection_origin = _validate_map_parameter( + "latitude_of_projection_origin", latitude_of_projection_origin + ) + self.false_easting = _validate_map_parameter( + "false_easting", false_easting + ) + self.false_northing = _validate_map_parameter( + "false_northing", false_northing + ) diff --git a/cf/gridmappings/abstract/cylindrical.py b/cf/gridmappings/abstract/cylindrical.py new file mode 100644 index 0000000000..47e44ee02a --- /dev/null +++ b/cf/gridmappings/abstract/cylindrical.py @@ -0,0 +1,38 @@ +from .gridmappingbase import ( + GridMapping, + _validate_map_parameter, +) + + +class CylindricalGridMapping(GridMapping): + """A Grid Mapping with Cylindrical classification. + + .. versionadded:: GMVER + + :Parameters: + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + """ + + def __init__(self, false_easting=0.0, false_northing=0.0, **kwargs): + super().__init__(**kwargs) + + self.false_easting = _validate_map_parameter( + "false_easting", false_easting + ) + self.false_northing = _validate_map_parameter( + "false_northing", false_northing + ) diff --git a/cf/gridmappings/abstract/gridmappingbase.py b/cf/gridmappings/abstract/gridmappingbase.py new file mode 100644 index 0000000000..b0179979fd --- /dev/null +++ b/cf/gridmappings/abstract/gridmappingbase.py @@ -0,0 +1,527 @@ +import itertools +import re + +from pyproj import CRS + +from ...constants import cr_canonical_units, cr_gm_valid_attr_names_are_numeric +from ...data import Data +from ...data.utils import is_numeric_dtype +from ...units import Units + + +PROJ_PREFIX = "+proj" + +""" +Define this first since it provides the default for several parameters, +e.g. WGS1984_CF_ATTR_DEFAULTS.semi_major_axis is 6378137.0, the radius +of the Earth in metres. Note we use the 'merc' projection to take these +from since that projection includes all the attributes given for +'latlon' instead and with identical values, but also includes further +map parameters with defaults applied across the projections. + +At the time of dedicating the code, the value of this is as follows, and +the values documented as defaults in the docstrings are taken from this: + +{'crs_wkt': '', + 'semi_major_axis': 6378137.0, + 'semi_minor_axis': 6356752.314245179, + 'inverse_flattening': 298.257223563, + 'reference_ellipsoid_name': 'WGS 84', + 'longitude_of_prime_meridian': 0.0, + 'prime_meridian_name': 'Greenwich', + 'geographic_crs_name': 'unknown', + 'horizontal_datum_name': 'World Geodetic System 1984', + 'projected_crs_name': 'unknown', + 'grid_mapping_name': 'mercator', + 'standard_parallel': 0.0, + 'longitude_of_projection_origin': 0.0, + 'false_easting': 0.0, + 'false_northing': 0.0, + 'scale_factor_at_projection_origin': 1.0} +""" +WGS1984_CF_ATTR_DEFAULTS = CRS.from_proj4("+proj=merc").to_cf() + + +def convert_proj_angular_data_to_cf(proj_data, context=None): + """Convert a PROJ angular data component into CF Data with CF Units. + + PROJ units for latitude and longitude are in + units of decimal degrees, where forming a string by adding + a suffix character indicates alternative units of + radians if the suffix is 'R' or 'r'. If a string, a suffix + of 'd', 'D' or '°' confirm units of decimal degrees. + + Note that `cf.convert_cf_angular_data_to_proj` and + this function are not strict inverse functions, + since the former will always convert to the *simplest* + way to specify the PROJ input, namely with no suffix for + degrees(_X) units and the 'R' suffix for radians, whereas + the input might have 'D' or 'r' etc. instead. + + .. versionadded:: GMVER + + :Parameters: + + proj_data: `str` + The PROJ angular data component, for example "90", "90.0", + "90.0D", "1.0R", or "1r". + + For details on valid PROJ data components in PROJ strings, + notably indicating units, see: + + https://proj.org/en/9.2/usage/projections.html#projection-units + + context: `str` or `None`, optional + The physical context of the conversion, where 'lat' indicates + a latitude value and 'lon' indicates a longitude, such that + indication of either context will return cf.Data with values + having units appropriate to that context, namely 'degrees_north' + or 'degrees_east' respectively. If None, 'degrees' or 'radians' + (depending on the input PROJ units) will be the units. + The default is None. + + :Returns: + + `Data` + A cf.Data object with CF-compliant units that corresponds + to the PROJ data and the context provided. + + """ + cf_compatible = True # unless find otherwise (True unless proven False) + if context == "lat": + cf_units = "degrees_north" + elif context == "lon": + cf_units = "degrees_east" + else: + # From the CF Conventions Document (3.1. Units): + # "The COARDS convention prohibits the unit degrees altogether, + # but this unit is not forbidden by the CF convention because + # it may in fact be appropriate for a variable containing, say, + # solar zenith angle. The unit degrees is also allowed on + # coordinate variables such as the latitude and longitude + # coordinates of a transformed grid. In this case the coordinate + # values are not true latitudes and longitudes which must always + # be identified using the more specific forms of degrees as + # described in Section 4.1. + cf_units = "degrees" + + # Only valid input is a valid float or integer (digit with zero or one + # decimal point only) optionally followed by a single suffix letter + # indicating decimal degrees or radians with PROJ. Be strict about an + # exact regex match, because anything not following the pattern (e.g. + # something with extra letters) will be ambiguous for PROJ units. + valid_form = re.compile(r"(-?\d+(\.\d*)?)([rRdD°]?)") + form = re.fullmatch(valid_form, proj_data) + if form: + comps = form.groups() + suffix = None + if len(comps) == 3: + value, float_comp, suffix = comps + else: + value, *float_comp = comps + + # Convert string value to relevant numeric form + if float_comp: + numeric_value = float(value) + else: + numeric_value = int(value) + + if suffix in ("r", "R"): # radians units + if context: + # Convert so we can apply degree_X form of the lat/lon context + numeric_value = Units.conform( + numeric_value, Units("radians"), Units("degrees") + ) + else: # Otherwise leave as radians to avoid rounding etc. + cf_units = "radians" + elif suffix and suffix not in ("d", "D", "°"): # 'decimal degrees' + cf_compatible = False + else: + cf_compatible = False + + if not cf_compatible: + raise ValueError( + f"PROJ data input not valid: {proj_data}. Ensure a valid " + "PROJ value and optionally units are supplied." + ) + + return Data(numeric_value, Units(cf_units)) + + +def convert_cf_angular_data_to_proj(data): + """Convert singleton angular CF Data into a PROJ data component. + + PROJ units for latitude and longitude are in + units of decimal degrees, where forming a string by adding + a suffix character indicates alternative units of + radians if the suffix is 'R' or 'r'. If a string, a suffix + of 'd', 'D' or '°' confirm units of decimal degrees. + + Note that this function and `convert_proj_angular_data_to_cf` + are not strict inverse functions, since this function + will always convert to the *simplest* + way to specify the PROJ input, namely with no suffix for + degrees(_X) units and the 'R' suffix for radians, whereas + the input might have 'D' or 'r' etc. instead. + + .. versionadded:: GMVER + + :Parameters: + + data: `Data` + A cf.Data object of size 1 containing an angular value + with CF-compliant angular units, for example + cf.Data(45, units="degrees_north"). + + :Returns: + + `str` + A PROJ angular data component that corresponds + to the Data provided. + + """ + if data.size != 1: + raise ValueError( + f"Input cf.Data must have size 1, got size: {data.size}" + ) + if not is_numeric_dtype(data): + raise TypeError( + f"Input cf.Data must have numeric data type, got: {data.dtype}" + ) + + units = data.Units + if not units: + raise ValueError( + "Must provide cf.Data with units for unambiguous conversion." + ) + units_str = units.units + + degrees_unit_prefix = ["degree", "degrees"] + # Taken from 4.1. Latitude Coordinate and 4.2. Longitude Coordinate: + # http://cfconventions.org/cf-conventions/cf-conventions.html... + # ...#latitude-coordinate and ...#longitude-coordinate + valid_cf_lat_lon_units = [ + s + e + for s, e in itertools.product( + degrees_unit_prefix, ("_north", "_N", "N", "_east", "_E", "E") + ) + ] + valid_degrees_units = degrees_unit_prefix + valid_cf_lat_lon_units + + if units_str in valid_degrees_units: + # No need for suffix 'D' for decimal degrees, as that is the default + # recognised when no suffix is given + proj_data = f"{data.data.array.item()}" + elif units_str == "radians": + proj_data = f"{data.data.array.item()}R" + else: + raise ValueError( + "Unrecognised angular units set on the cf.Data. Valid options " + f"are: {', '.join(valid_degrees_units)} and radians but got: " + f"{units_str}" + ) + + return proj_data + + +def _make_proj_string_comp(spec): + """Form a PROJ proj-string end from the given PROJ parameters. + + :Parameters: + + spec: `dict` + A dictionary providing the proj-string specifiers for + parameters, as keys, with their values as values. Values + must be convertible to strings. + + """ + proj_string = "" + for comp, value in spec.items(): + if not isinstance(value, str): + try: + value = str(value) + except TypeError: + raise TypeError( + "Can't create proj-string due to non-representable " + f"value {value} for key {comp}" + ) + proj_string += f" +{comp}={value}" + return proj_string + + +def _validate_map_parameter(mp_name, mp_value): + """Validate map parameters for correct type and canonical units. + + :Parameters: + + mp_name: `str` + The name of the map parameter to validate. It should be + a name valid in this way in CF, namely one listed under + the 'Table F.1. Grid Mapping Attributes' in Appendix + F: Grid Mappings', therefore listed as a key in + `cr_gm_valid_attr_names_are_numeric`, else a ValueError + will be raised early. + + mp_value: `str`, `Data`, numeric or `None` + The map parameter value being set for the given map + parameter name. The type, and if numeric or `Data`, + units, will be validated against the expected values + for the given map parameter name. + + :Returns: + + `Data`, `str`, or `None` + The map parameter value, assuming it passes validation, + conformed to the canonical units of the map parameter + name, if units are applicable. + + """ + # 0. Check input parameters are valid CF GM map parameters, not + # for any case something unrecognised that will silently do nothing. + if mp_name not in cr_gm_valid_attr_names_are_numeric: + raise ValueError( + "Unrecognised map parameter provided for the " + f"Grid Mapping: {mp_name}" + ) + + # 1. If None, can return early: + if mp_value is None: # distinguish from 0 or 0.0 etc. + return None + + # 2. Now ensure the type of the value is as expected. + expect_numeric = cr_gm_valid_attr_names_are_numeric[mp_name] + if expect_numeric: + if (isinstance(mp_value, Data) and not is_numeric_dtype(mp_value)) or ( + not isinstance(mp_value, (Data, int, float)) + ): + raise TypeError( + f"Map parameter {mp_name} has an incompatible " + "data type, expected numeric or Data but got " + f"{type(mp_value)}" + ) + elif not expect_numeric and not isinstance(mp_value, str): + raise TypeError( + f"Map parameter {mp_name} has an incompatible " + "data type, expected a string but got " + f"{type(mp_value)}" + ) + + # 3. Finally ensure the units are valid and conformed to the + # canonical units, where numeric. + if expect_numeric: + canon_units = cr_canonical_units[mp_name] + if isinstance(mp_value, Data): + # In this case, is Data which may have units which aren't equal to + # the canonical ones, so may need to conform the value. + units = mp_value.Units + # The units must be checked and might need to be conformed + if not units.equivalent(canon_units): + raise ValueError( + f"Map parameter {mp_name} value has units that " + "are incompatible with the expected units of " + f"{canon_units}: {units}" + ) + elif not units.equals(canon_units): + conforming_value = Units.conform( + mp_value.array.item(), units, canon_units + ) + else: + conforming_value = mp_value + else: + conforming_value = mp_value + + # Return numeric value as Data with conformed value and canonical units + return Data(conforming_value, units=canon_units) + else: + # Return string value + return mp_value + + +class GridMapping(): + """A container for a Grid Mapping recognised by the CF Conventions.""" + + # The value of the 'grid_mapping_name' attribute. + grid_mapping_name = None + # The PROJ projection identifier shorthand name. + proj_id = None + + def __init__( + self, + # i.e. WGS1984_CF_ATTR_DEFAULTS["reference_ellipsoid_name"], etc. + reference_ellipsoid_name="WGS 84", + # The next three parameters are non-zero floats so don't hard-code + # WGS84 defaults in case of future precision changes: + semi_major_axis=WGS1984_CF_ATTR_DEFAULTS["semi_major_axis"], + semi_minor_axis=WGS1984_CF_ATTR_DEFAULTS["semi_minor_axis"], + inverse_flattening=WGS1984_CF_ATTR_DEFAULTS["inverse_flattening"], + prime_meridian_name="Greenwich", + longitude_of_prime_meridian=0.0, + earth_radius=None, + **kwargs, + ): + """**Initialisation** + + :Parameters: + + reference_ellipsoid_name: `str` or `None`, optional + The name of a built-in ellipsoid definition. + The default is "WGS 84". + + .. note:: If used in conjunction with 'earth_radius', + the 'earth_radius' parameter takes precedence. + + inverse_flattening: number or scalar `Data`, optional + The reverse flattening of the ellipsoid (PROJ 'rf' + value), :math:`\frac{1}{f}`, where f corresponds to + the flattening value (PROJ 'f' value) for the + ellipsoid. Unitless, so `Data` must be unitless. + The default is 298.257223563. + + prime_meridian_name: `str`, optional + A predeclared name to define the prime meridian (PROJ + 'pm' value). The default is "Greenwich". Supported + names and corresponding longitudes are listed at: + + https://proj.org/en/9.2/usage/ + projections.html#prime-meridian + + .. note:: If used in conjunction with + 'longitude_of_prime_meridian', this + parameter takes precedence. + + longitude_of_prime_meridian: number or scalar `Data`, optional + The longitude relative to Greenwich of the + prime meridian. If provided as a number or `Data` without + units, the units are taken as 'degrees_east', else the + `Data` units are taken and must be angular and + compatible with longitude. The default is 0.0 + degrees_east. + + .. note:: If used in conjunction with + 'prime_meridian_name', the + 'prime_meridian_name' parameter takes + precedence. + + semi_major_axis: number, scalar `Data` or `None`, optional + The semi-major axis of the ellipsoid (PROJ 'a' value), + in units of meters unless units are otherwise specified + via `Data` units, in which case they must be conformable + to meters and will be converted. The default is 6378137.0. + + semi_minor_axis: number, scalar `Data` or `None`, optional + The semi-minor axis of the ellipsoid (PROJ 'b' value) + in units of meters unless units are otherwise specified + via `Data` units, in which case they must be conformable + to meters and will be converted. The default is + 6356752.314245179. + + earth_radius: number, scalar `Data` or `None`, optional + The radius of the ellipsoid, if a sphere (PROJ 'R' value), + in units of meters unless units are otherwise specified + via `Data` units, in which case they must be conformable + to meters and will be converted. + + If the ellipsoid is not a sphere, then + set as `None`, the default, to indicate that ellipsoid + parameters such as the reference_ellipsoid_name or + semi_major_axis and semi_minor_axis are being set, + since these will should be used to define a + non-spherical ellipsoid. + + .. note:: If used in conjunction with + 'reference_ellipsoid_name', this parameter + takes precedence. + + """ + # Validate map parameters that are valid for any GridMapping. + # These are attributes which describe the ellipsoid and prime meridian, + # which may be included, when applicable, with any grid mapping, as + # specified in Appendix F of the Conventions. + self.earth_radius = _validate_map_parameter( + "earth_radius", earth_radius + ) + self.inverse_flattening = _validate_map_parameter( + "inverse_flattening", inverse_flattening + ) + self.longitude_of_prime_meridian = _validate_map_parameter( + "longitude_of_prime_meridian", longitude_of_prime_meridian + ) + self.prime_meridian_name = _validate_map_parameter( + "prime_meridian_name", prime_meridian_name + ) + self.reference_ellipsoid_name = _validate_map_parameter( + "reference_ellipsoid_name", reference_ellipsoid_name + ) + self.semi_major_axis = _validate_map_parameter( + "semi_major_axis", semi_major_axis + ) + self.semi_minor_axis = _validate_map_parameter( + "semi_minor_axis", semi_minor_axis + ) + + @classmethod + def is_latlon_gm(cls): + """Whether the Grid Mapping is of LatitudeLongitude form. + + :Returns: + + `bool` + True only if the Grid Mapping is LatitudeLongitude. + + """ + return False + + def __repr__(self): + """x.__repr__() <==> repr(x)""" + # Report parent GridMapping class to indicate classification, + # but only if it has one (> 2 avoids own class and 'object') + # base. E.g. we get , + # , . + parent_gm = "" + if len(self.__class__.__mro__) > 2: + parent_gm = self.__class__.__mro__[1].__name__ + ": " + return f"" + + def __str__(self): + """x.__str__() <==> str(x)""" + return f"{self.__repr__()[:-1]} {self.get_proj_string()}>" + + def __eq__(self, other): + """The rich comparison operator ``==``.""" + return self.get_proj_crs() == other.get_proj_crs() + + def __hash__(self, other): + """The built-in function `hash`, x.__hash__() <==> hash(x).""" + return hash(self.get_proj_crs()) + + def get_proj_string(self, params=None): + """The value of the PROJ proj-string defining the projection.""" + # TODO enable parameter input and sync'ing + if not params: + params = "" + return f"{PROJ_PREFIX}={self.proj_id}{params}" + + def get_proj_crs(self): + """Get the PROJ Coordinate Reference System. + + :Returns: + + `pyproj.crs.CRS` + The PROJ Coordinate Reference System defined with + a `pyproj` `CRS` class that corresponds to the + Grid Mapping instance. + + """ + return CRS.from_proj4(self.get_proj_string()) + + def has_crs_wkt(self): + """True if the Grid Mapping has a valid crs_wkt attribute set. + + :Returns: + + `bool` + Whether the Grid Mapping instance has a crs_wkt + attribute set. + + """ + return self.crs_wkt is not None diff --git a/cf/gridmappings/abstract/latlon.py b/cf/gridmappings/abstract/latlon.py new file mode 100644 index 0000000000..44381df9c4 --- /dev/null +++ b/cf/gridmappings/abstract/latlon.py @@ -0,0 +1,18 @@ +from .gridmappingbase import ( + GridMapping, + _validate_map_parameter, +) + + +class LatLonGridMapping(GridMapping): + """A Grid Mapping with Latitude-Longitude nature. + + Such a Grid Mapping is based upon latitude and longitude coordinates + on a spherical Earth, defining the canonical 2D geographical coordinate + system so that the figure of the Earth can be described. + + .. versionadded:: GMVER + + """ + + pass diff --git a/cf/gridmappings/abstract/perspective.py b/cf/gridmappings/abstract/perspective.py new file mode 100644 index 0000000000..0a3bd38a16 --- /dev/null +++ b/cf/gridmappings/abstract/perspective.py @@ -0,0 +1,27 @@ +from .gridmappingbase import _validate_map_parameter +from .azimuthal import AzimuthalGridMapping + + +class PerspectiveGridMapping(AzimuthalGridMapping): + """A Grid Mapping with Azimuthal classification and perspective view. + + .. versionadded:: GMVER + + :Parameters: + + perspective_point_height: number or scalar `Data` + The height of the view point above the surface (PROJ + 'h') value, for example the height of a satellite above + the Earth, in units of meters 'm'. If provided + as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + + """ + + def __init__(self, perspective_point_height, **kwargs): + super().__init__(**kwargs) + + self.perspective_point_height = _validate_map_parameter( + "perspective_point_height", perspective_point_height + ) diff --git a/cf/gridmappings/albersequalarea.py b/cf/gridmappings/albersequalarea.py new file mode 100644 index 0000000000..de2e80a58b --- /dev/null +++ b/cf/gridmappings/albersequalarea.py @@ -0,0 +1,66 @@ +from .abstract import ConicGridMapping + + +class AlbersEqualArea(ConicGridMapping): + """The Albers Equal Area grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#_albers_equal_area + + or the corresponding PROJ projection page: + + https://proj.org/en/9.2/operations/projections/aea.html + + for more information. + + .. versionadded:: GMVER + + :Parameters: + + standard_parallel: 2-`tuple` of number or scalar `Data` or `None` + The standard parallel value(s): the first (PROJ 'lat_1' + value) and/or the second (PROJ 'lat_2' value), given + as a 2-tuple of numbers or strings corresponding to + the first and then the second in order, where `None` + indicates that a value is not being set for either. + + If provided as a number or `Data` without units, the units + for each of the values are taken as 'degrees_north', else + the `Data` units are taken and must be angular and + compatible with latitude. + + longitude_of_central_meridian: number or scalar `Data`, optional + The longitude of (natural) origin i.e. central meridian. + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + latitude_of_projection_origin: number or scalar `Data`, optional + The latitude of projection center (PROJ 'lat_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_north', else the `Data` units are + taken and must be angular and compatible with latitude. + The default is 0.0 degrees_north. + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + """ + + grid_mapping_name = "albers_conical_equal_area" + proj_id = "aea" diff --git a/cf/gridmappings/azimuthalequidistant.py b/cf/gridmappings/azimuthalequidistant.py new file mode 100644 index 0000000000..7b322a27cf --- /dev/null +++ b/cf/gridmappings/azimuthalequidistant.py @@ -0,0 +1,54 @@ +from .abstract import AzimuthalGridMapping + + +class AzimuthalEquidistant(AzimuthalGridMapping): + """The Azimuthal Equidistant grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#azimuthal-equidistant + + or the corresponding PROJ projection page: + + https://proj.org/en/9.2/operations/projections/aeqd.html + + for more information. + + .. versionadded:: GMVER + + :Parameters: + + longitude_of_projection_origin: number or scalar `Data`, optional + The longitude of projection center (PROJ 'lon_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + latitude_of_projection_origin: number or scalar `Data`, optional + The latitude of projection center (PROJ 'lat_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_north', else the `Data` units are + taken and must be angular and compatible with latitude. + The default is 0.0 degrees_north. + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + """ + + grid_mapping_name = "azimuthal_equidistant" + proj_id = "aeqd" diff --git a/cf/gridmappings/geostationary.py b/cf/gridmappings/geostationary.py new file mode 100644 index 0000000000..578b61a509 --- /dev/null +++ b/cf/gridmappings/geostationary.py @@ -0,0 +1,121 @@ +from .abstract import PerspectiveGridMapping +from .abstract.gridmappingbase import _validate_map_parameter + + +class Geostationary(PerspectiveGridMapping): + """The Geostationary Satellite View grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#_geostationary_projection + + or the corresponding PROJ projection page: + + https://proj.org/en/9.2/operations/projections/geos.html + + for more information. + + .. versionadded:: GMVER + + :Parameters: + + perspective_point_height: number or scalar `Data` + The height of the view point above the surface (PROJ + 'h') value, for example the height of a satellite above + the Earth, in units of meters 'm'. If provided + as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + + longitude_of_projection_origin: number or scalar `Data`, optional + The longitude of projection center (PROJ 'lon_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + latitude_of_projection_origin: number or scalar `Data`, optional + The latitude of projection center (PROJ 'lat_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_north', else the `Data` units are + taken and must be angular and compatible with latitude. + The default is 0.0 degrees_north. + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + sweep_angle_axis: `str`, optional + Sweep angle axis of the viewing instrument, which indicates + the axis on which the view sweeps. Valid options + are "x" and "y". The default is "y". + + For more information about the nature of this parameter, see: + + https://proj.org/en/9.2/operations/projections/ + geos.html#note-on-sweep-angle + + fixed_angle_axis: `str`, optional + The axis on which the view is fixed. It corresponds to the + inner-gimbal axis of the gimbal view model, whose axis of + rotation moves about the outer-gimbal axis. Valid options + are "x" and "y". The default is "x". + + .. note:: If the fixed_angle_axis is "x", sweep_angle_axis + is "y", and vice versa. + + """ + + grid_mapping_name = "geostationary" + proj_id = "geos" + + def __init__( + self, + perspective_point_height, + longitude_of_projection_origin=0.0, + latitude_of_projection_origin=0.0, + false_easting=0.0, + false_northing=0.0, + sweep_angle_axis="y", + fixed_angle_axis="x", + **kwargs, + ): + super().__init__( + perspective_point_height, + longitude_of_projection_origin=longitude_of_projection_origin, + latitude_of_projection_origin=latitude_of_projection_origin, + false_easting=false_easting, + false_northing=false_northing, + **kwargs, + ) + + # Values "x" and "y" are not case-sensitive, so convert to lower-case + self.sweep_angle_axis = _validate_map_parameter( + "sweep_angle_axis", sweep_angle_axis + ).lower() + self.fixed_angle_axis = _validate_map_parameter( + "fixed_angle_axis", fixed_angle_axis + ).lower() + + # sweep_angle_axis must be the opposite (of "x" and "y") to + # fixed_angle_axis. + if (self.sweep_angle_axis, self.fixed_angle_axis) not in [ + ("x", "y"), + ("y", "x"), + ]: + raise ValueError( + "The sweep_angle_axis must be the opposite value, from 'x' " + "and 'y', to the fixed_angle_axis." + ) diff --git a/cf/gridmappings/gridmapping.py b/cf/gridmappings/gridmapping.py new file mode 100644 index 0000000000..65e93d809f --- /dev/null +++ b/cf/gridmappings/gridmapping.py @@ -0,0 +1,130 @@ +# Concrete classes for all Grid Mappings supported by the CF Conventions. +# For the full listing, see: +# https://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ +# cf-conventions.html#appendix-grid-mappings +# from which these classes should be kept consistent and up-to-date. +from .albersequalarea import AlbersEqualArea +from .azimuthalequidistant import AzimuthalEquidistant +from .geostationary import Geostationary +from .lambertazimuthalequalarea import LambertAzimuthalEqualArea +from .lambertconformalconic import LambertConformalConic +from .lambertcylindricalequalarea import LambertCylindricalEqualArea +from .latitudelongitude import LatitudeLongitude +from .mercator import Mercator +from .obliquemercator import ObliqueMercator +from .orthographic import Orthographic +from .polarstereographic import PolarStereographic +from .rotatedlatitudelongitude import RotatedLatitudeLongitude +from .sinusoidal import Sinusoidal +from .stereographic import Stereographic +from .transversemercator import TransverseMercator +from .verticalperspective import VerticalPerspective + + +class InvalidGridMapping(Exception): + """Exception for a coordinate reference with unsupported Grid Mapping. + + .. versionadded:: GMVER + + :Parameters: + + TODOGM. + + """ + + def __init__(self, gm_name_attr, custom_message=None): + self.gm_name_attr = gm_name_attr + self.custom_message = custom_message + + def __str__(self): + if self.custom_message: + return self.custom_message + elif self.gm_name_attr: + return ( + f"Coordinate reference construct with grid_mapping_name " + f"{self.gm_name_attr} corresponds to a Grid Mapping that " + "is not supported by the CF Conventions." + ) + else: + return ( + f"A coordinate reference construct must have an attribute " + "'Coordinate conversion:grid_mapping_name' defined in order " + "to interpret it as a GM class. Missing 'grid_mapping_name'." + ) + + +class GM(): + """A validated Grid Mapping supported by the CF Conventions. + + .. versionadded:: GMVER + + :Parameters: + TODOGM. + + """ + + def __new__(cls, coordinate_reference, **override_kwargs): + """TODOGM.""" + if cls is GM: + # If there is no coordinate conversion or parameters set on one, + # None will result from the pair of lines below, so is safe. + conv = coordinate_reference.get_coordinate_conversion() + name = conv.get_parameter( + "grid_mapping_name", default=None + ) + + if not name: # Exit early before further parameter querying + raise InvalidGridMapping(name) + + kwargs = conv.parameters() + kwargs.pop("grid_mapping_name") # will be set in creation + # Left with those parameters to describe the GM, which must be + # validated as map parameters. Note x.update(y) will override x + # with any keys from y that are duplicates, avoiding issues + # from duplication of keyword inputs + kwargs.update(override_kwargs) + + # TODO: once cf Python minimum is v.3.10, use the new match/case + # syntax to consolidate this long if/elif. + if name == "albers_conical_equal_area": + return AlbersEqualArea(**kwargs) + elif name == "azimuthal_equidistant": + return AzimuthalEquidistant(**kwargs) + elif name == "geostationary": + return Geostationary(**kwargs) + elif name == "lambert_azimuthal_equal_area": + return LambertAzimuthalEqualArea(**kwargs) + elif name == "lambert_conformal_conic": + return LambertConformalConic(**kwargs) + elif name == "lambert_cylindrical_equal_area": + return LambertCylindricalEqualArea(**kwargs) + elif name == "latitude_longitude": + return LatitudeLongitude(**kwargs) + elif name == "mercator": + return Mercator(**kwargs) + elif name == "oblique_mercator": + return ObliqueMercator(**kwargs) + elif name == "orthographic": + return Orthographic(**kwargs) + elif name == "polar_stereographic": + return PolarStereographic(**kwargs) + elif name == "rotated_latitude_longitude": + return RotatedLatitudeLongitude(**kwargs) + elif name == "sinusoidal": + return Sinusoidal(**kwargs) + elif name == "stereographic": + return Stereographic(**kwargs) + elif name == "transverse_mercator": + return TransverseMercator(**kwargs) + elif name == "vertical_perspective": + return VerticalPerspective(**kwargs) + else: + raise InvalidGridMapping(name) + + def __init__(self, coordinate_reference): + """TODOGM.""" + pass # TODOGM + + def create_crs(): + """TODOGM.""" + pass # TODOGM diff --git a/cf/gridmappings/lambertazimuthalequalarea.py b/cf/gridmappings/lambertazimuthalequalarea.py new file mode 100644 index 0000000000..94b23aefb4 --- /dev/null +++ b/cf/gridmappings/lambertazimuthalequalarea.py @@ -0,0 +1,54 @@ +from .abstract import AzimuthalGridMapping + + +class LambertAzimuthalEqualArea(AzimuthalGridMapping): + """The Lambert Azimuthal Equal Area grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#lambert-azimuthal-equal-area + + or the corresponding PROJ projection page: + + https://proj.org/en/9.2/operations/projections/laea.html + + for more information. + + .. versionadded:: GMVER + + :Parameters: + + longitude_of_projection_origin: number or scalar `Data`, optional + The longitude of projection center (PROJ 'lon_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + latitude_of_projection_origin: number or scalar `Data`, optional + The latitude of projection center (PROJ 'lat_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_north', else the `Data` units are + taken and must be angular and compatible with latitude. + The default is 0.0 degrees_north. + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + """ + + grid_mapping_name = "lambert_azimuthal_equal_area" + proj_id = "laea" diff --git a/cf/gridmappings/lambertconformalconic.py b/cf/gridmappings/lambertconformalconic.py new file mode 100644 index 0000000000..6a0e9e4e61 --- /dev/null +++ b/cf/gridmappings/lambertconformalconic.py @@ -0,0 +1,69 @@ +from .abstract import ConicGridMapping + + +class LambertConformalConic(ConicGridMapping): + """The Lambert Conformal Conic grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#_lambert_conformal + + or the corresponding PROJ projection page: + + https://proj.org/en/9.2/operations/projections/lcc.html + + for more information. + + .. versionadded:: GMVER + + :Parameters: + + standard_parallel: 2-`tuple` of number or scalar `Data` or `None` + The standard parallel value(s): the first (PROJ 'lat_1' + value) and/or the second (PROJ 'lat_2' value), given + as a 2-tuple of numbers or strings corresponding to + the first and then the second in order, where `None` + indicates that a value is not being set for either. + + If provided as a number or `Data` without units, the units + for each of the values are taken as 'degrees_north', else + the `Data` units are taken and must be angular and + compatible with latitude. + + The default is (0.0, 0.0), that is 0.0 degrees_north + for the first and second standard parallel values. + + longitude_of_projection_origin: number or scalar `Data`, optional + The longitude of projection center (PROJ 'lon_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + latitude_of_projection_origin: number or scalar `Data`, optional + The latitude of projection center (PROJ 'lat_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_north', else the `Data` units are + taken and must be angular and compatible with latitude. + The default is 0.0 degrees_north. + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + """ + + grid_mapping_name = "lambert_conformal_conic" + proj_id = "lcc" diff --git a/cf/gridmappings/lambertcylindricalequalarea.py b/cf/gridmappings/lambertcylindricalequalarea.py new file mode 100644 index 0000000000..06c0b3c246 --- /dev/null +++ b/cf/gridmappings/lambertcylindricalequalarea.py @@ -0,0 +1,95 @@ +from .abstract import CylindricalGridMapping +from .abstract.gridmappingbase import _validate_map_parameter + + +class LambertCylindricalEqualArea(CylindricalGridMapping): + """The Equal Area Cylindrical grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#_lambert_cylindrical_equal_area + + or the corresponding PROJ projection page: + + https://proj.org/en/9.2/operations/projections/cea.html + + for more information. + + .. versionadded:: GMVER + + :Parameters: + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + standard_parallel: 2-`tuple` of number or scalar `Data` + or `None`, optional + The standard parallel value(s): the first (PROJ 'lat_1' + value) and/or the second (PROJ 'lat_2' value), given + as a 2-tuple of numbers or strings corresponding to + the first and then the second in order, where `None` + indicates that a value is not being set for either. + + If provided as a number or `Data` without units, the units + for each of the values are taken as 'degrees_north', else + the `Data` units are taken and must be angular and + compatible with latitude. + + The default is (0.0, 0.0), that is 0.0 degrees_north + for the first and second standard parallel values. + + longitude_of_central_meridian: number or scalar `Data`, optional + The longitude of (natural) origin i.e. central meridian. + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + scale_factor_at_projection_origin: number or scalar `Data`, optional + The scale factor used in the projection (PROJ 'k_0' value). + Unitless, so `Data` must be unitless. The default is 1.0. + + """ + + grid_mapping_name = "lambert_cylindrical_equal_area" + proj_id = "cea" + + def __init__( + self, + false_easting=0.0, + false_northing=0.0, + standard_parallel=(0.0, 0.0), + scale_factor_at_projection_origin=1.0, + longitude_of_central_meridian=0.0, + **kwargs, + ): + super().__init__( + false_easting=false_easting, + false_northing=false_northing, + **kwargs, + ) + + self.standard_parallel = ( + _validate_map_parameter("standard_parallel", standard_parallel[0]), + _validate_map_parameter("standard_parallel", standard_parallel[1]), + ) + self.longitude_of_central_meridian = _validate_map_parameter( + "longitude_of_central_meridian", longitude_of_central_meridian + ) + self.scale_factor_at_projection_origin = _validate_map_parameter( + "scale_factor_at_projection_origin", + scale_factor_at_projection_origin, + ) diff --git a/cf/gridmappings/latitudelongitude.py b/cf/gridmappings/latitudelongitude.py new file mode 100644 index 0000000000..4844f527c1 --- /dev/null +++ b/cf/gridmappings/latitudelongitude.py @@ -0,0 +1,32 @@ +from .abstract import LatLonGridMapping + + +class LatitudeLongitude(LatLonGridMapping): + """The Latitude-Longitude grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#_latitude_longitude + + for more information. + + .. versionadded:: GMVER + + """ + + grid_mapping_name = "latitude_longitude" + proj_id = "latlong" + + @classmethod + def is_latlon_gm(cls): + """Whether the Grid Mapping is of LatitudeLongitude form. + + :Returns: + + `bool` + True only if the Grid Mapping is LatitudeLongitude. + + """ + return True diff --git a/cf/gridmappings/mercator.py b/cf/gridmappings/mercator.py new file mode 100644 index 0000000000..117f5a411c --- /dev/null +++ b/cf/gridmappings/mercator.py @@ -0,0 +1,95 @@ +from .abstract import CylindricalGridMapping +from .abstract.gridmappingbase import _validate_map_parameter + + +class Mercator(CylindricalGridMapping): + """The Mercator grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#_mercator + + or the corresponding PROJ projection page: + + https://proj.org/en/9.2/operations/projections/merc.html + + for more information. + + .. versionadded:: GMVER + + :Parameters: + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + standard_parallel: 2-`tuple` of number or scalar `Data` + or `None`, optional + The standard parallel value(s): the first (PROJ 'lat_1' + value) and/or the second (PROJ 'lat_2' value), given + as a 2-tuple of numbers or strings corresponding to + the first and then the second in order, where `None` + indicates that a value is not being set for either. + + If provided as a number or `Data` without units, the units + for each of the values are taken as 'degrees_north', else + the `Data` units are taken and must be angular and + compatible with latitude. + + The default is (0.0, 0.0), that is 0.0 degrees_north + for the first and second standard parallel values. + + longitude_of_projection_origin: number or scalar `Data`, optional + The longitude of projection center (PROJ 'lon_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + scale_factor_at_projection_origin: number or scalar `Data`, optional + The scale factor used in the projection (PROJ 'k_0' value). + Unitless, so `Data` must be unitless. The default is 1.0. + + """ + + grid_mapping_name = "mercator" + proj_id = "merc" + + def __init__( + self, + false_easting=0.0, + false_northing=0.0, + standard_parallel=(0.0, 0.0), + longitude_of_projection_origin=0.0, + scale_factor_at_projection_origin=1.0, + **kwargs, + ): + super().__init__( + false_easting=false_easting, + false_northing=false_northing, + **kwargs, + ) + + self.standard_parallel = ( + _validate_map_parameter("standard_parallel", standard_parallel[0]), + _validate_map_parameter("standard_parallel", standard_parallel[1]), + ) + self.longitude_of_projection_origin = _validate_map_parameter( + "longitude_of_projection_origin", longitude_of_projection_origin + ) + self.scale_factor_at_projection_origin = _validate_map_parameter( + "scale_factor_at_projection_origin", + scale_factor_at_projection_origin, + ) diff --git a/cf/gridmappings/obliquemercator.py b/cf/gridmappings/obliquemercator.py new file mode 100644 index 0000000000..310906a0eb --- /dev/null +++ b/cf/gridmappings/obliquemercator.py @@ -0,0 +1,97 @@ +from .abstract import CylindricalGridMapping +from .abstract.gridmappingbase import _validate_map_parameter + + +class ObliqueMercator(CylindricalGridMapping): + """The Oblique Mercator grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#_oblique_mercator + + or the corresponding PROJ projection page: + + https://proj.org/en/9.2/operations/projections/omerc.html + + for more information. + + .. versionadded:: GMVER + + :Parameters: + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + azimuth_of_central_line: number or scalar `Data`, optional + The azimuth i.e. tilt angle of the centerline clockwise + from north at the center point of the line (PROJ 'alpha' + value). If provided as a number or `Data` without units, + the units are taken as 'degrees', else the `Data` + units are taken and must be angular and compatible. + The default is 0.0 degrees. + + longitude_of_projection_origin: number or scalar `Data`, optional + The longitude of projection center (PROJ 'lon_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + latitude_of_projection_origin: number or scalar `Data`, optional + The latitude of projection center (PROJ 'lat_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_north', else the `Data` units are + taken and must be angular and compatible with latitude. + The default is 0.0 degrees_north. + + scale_factor_at_projection_origin: number or scalar `Data`, optional + The scale factor used in the projection (PROJ 'k_0' value). + Unitless, so `Data` must be unitless. The default is 1.0. + + """ + + grid_mapping_name = "oblique_mercator" + proj_id = "omerc" + + def __init__( + self, + azimuth_of_central_line=0.0, + latitude_of_projection_origin=0.0, + longitude_of_projection_origin=0.0, + scale_factor_at_projection_origin=1.0, + false_easting=0.0, + false_northing=0.0, + **kwargs, + ): + super().__init__( + false_easting=false_easting, + false_northing=false_northing, + **kwargs, + ) + + self.azimuth_of_central_line = _validate_map_parameter( + "azimuth_of_central_line", azimuth_of_central_line + ) + self.latitude_of_projection_origin = _validate_map_parameter( + "latitude_of_projection_origin", latitude_of_projection_origin + ) + self.longitude_of_projection_origin = _validate_map_parameter( + "longitude_of_projection_origin", longitude_of_projection_origin + ) + self.scale_factor_at_projection_origin = _validate_map_parameter( + "scale_factor_at_projection_origin", + scale_factor_at_projection_origin, + ) diff --git a/cf/gridmappings/orthographic.py b/cf/gridmappings/orthographic.py new file mode 100644 index 0000000000..0ef0eea14a --- /dev/null +++ b/cf/gridmappings/orthographic.py @@ -0,0 +1,54 @@ +from .abstract import AzimuthalGridMapping + + +class Orthographic(AzimuthalGridMapping): + """The Orthographic grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#_orthographic + + or the corresponding PROJ projection page: + + https://proj.org/en/9.2/operations/projections/ortho.html + + for more information. + + .. versionadded:: GMVER + + :Parameters: + + longitude_of_projection_origin: number or scalar `Data`, optional + The longitude of projection center (PROJ 'lon_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + latitude_of_projection_origin: number or scalar `Data`, optional + The latitude of projection center (PROJ 'lat_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_north', else the `Data` units are + taken and must be angular and compatible with latitude. + The default is 0.0 degrees_north. + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + """ + + grid_mapping_name = "orthographic" + proj_id = "ortho" diff --git a/cf/gridmappings/polarstereographic.py b/cf/gridmappings/polarstereographic.py new file mode 100644 index 0000000000..29caad1096 --- /dev/null +++ b/cf/gridmappings/polarstereographic.py @@ -0,0 +1,127 @@ +from .abstract import AzimuthalGridMapping +from .abstract.gridmappingbase import _validate_map_parameter + + +class PolarStereographic(AzimuthalGridMapping): + """The Universal Polar Stereographic grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#polar-stereographic + + or the corresponding PROJ projection page: + + https://proj.org/en/9.2/operations/projections/ups.html + + for more information. + + .. versionadded:: GMVER + + :Parameters: + + straight_vertical_longitude_from_pole: number or scalar `Data`, optional + The longitude of (natural) origin i.e. central meridian, + oriented straight up from the North or South Pole. + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + longitude_of_projection_origin: number or scalar `Data`, optional + The longitude of projection center (PROJ 'lon_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + latitude_of_projection_origin: number or scalar `Data`, optional + The latitude of projection center (PROJ 'lat_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_north', else the `Data` units are + taken and must be angular and compatible with latitude. + The default is 0.0 degrees_north. + + scale_factor_at_projection_origin: number or scalar `Data`, optional + The scale factor used in the projection (PROJ 'k_0' value). + Unitless, so `Data` must be unitless. The default is 1.0. + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + standard_parallel: 2-`tuple` of number or scalar `Data` + or `None`, optional + The standard parallel value(s): the first (PROJ 'lat_1' + value) and/or the second (PROJ 'lat_2' value), given + as a 2-tuple of numbers or strings corresponding to + the first and then the second in order, where `None` + indicates that a value is not being set for either. + + If provided as a number or `Data` without units, the units + for each of the values are taken as 'degrees_north', else + the `Data` units are taken and must be angular and + compatible with latitude. + + The default is (0.0, 0.0), that is 0.0 degrees_north + for the first and second standard parallel values. + + """ + + grid_mapping_name = "polar_stereographic" + proj_id = "ups" + + def __init__( + self, + latitude_of_projection_origin=0.0, + longitude_of_projection_origin=0.0, + false_easting=0.0, + false_northing=0.0, + standard_parallel=(0.0, 0.0), + straight_vertical_longitude_from_pole=0.0, + scale_factor_at_projection_origin=1.0, + **kwargs, + ): + # TODO check defaults here, they do not appear for + # CRS.from_proj4("+proj=ups").to_cf() to cross reference! + super().__init__( + latitude_of_projection_origin=latitude_of_projection_origin, + longitude_of_projection_origin=longitude_of_projection_origin, + false_easting=false_easting, + false_northing=false_northing, + **kwargs, + ) + + # See: https://github.com/cf-convention/cf-conventions/issues/445 + if ( + longitude_of_projection_origin + and straight_vertical_longitude_from_pole + ): + raise ValueError( + "Only one of 'longitude_of_projection_origin' and " + "'straight_vertical_longitude_from_pole' can be set." + ) + + self.straight_vertical_longitude_from_pole = _validate_map_parameter( + "straight_vertical_longitude_from_pole", + straight_vertical_longitude_from_pole, + ) + self.standard_parallel = ( + _validate_map_parameter("standard_parallel", standard_parallel[0]), + _validate_map_parameter("standard_parallel", standard_parallel[1]), + ) + self.scale_factor_at_projection_origin = _validate_map_parameter( + "scale_factor_at_projection_origin", + scale_factor_at_projection_origin, + ) diff --git a/cf/gridmappings/rotatedlatitudelongitude.py b/cf/gridmappings/rotatedlatitudelongitude.py new file mode 100644 index 0000000000..eaba83e2eb --- /dev/null +++ b/cf/gridmappings/rotatedlatitudelongitude.py @@ -0,0 +1,62 @@ +from .abstract import LatLonGridMapping +from .abstract.gridmappingbase import _validate_map_parameter + + +class RotatedLatitudeLongitude(LatLonGridMapping): + """The Rotated Latitude-Longitude grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#_rotated_pole + + for more information. + + .. versionadded:: GMVER + + :Parameters: + + grid_north_pole_latitude: number or scalar `Data` + Latitude of the North pole of the unrotated source CRS, + expressed in the rotated geographic CRS. + If provided as a number or `Data` without units, the units + are taken as 'degrees_north', else the `Data` units are + taken and must be angular and compatible with latitude. + + grid_north_pole_longitude: number or scalar `Data` + Longitude of the North pole of the unrotated source CRS, + expressed in the rotated geographic CRS. + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + + north_pole_grid_longitude: number or scalar `Data`, optional + The longitude of projection center (PROJ 'lon_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + + """ + + grid_mapping_name = "rotated_latitude_longitude" + proj_id = "latlong" + + def __init__( + self, + grid_north_pole_latitude, + grid_north_pole_longitude, + north_pole_grid_longitude=0.0, + **kwargs, + ): + super().__init__(**kwargs) + + self.grid_north_pole_latitude = _validate_map_parameter( + "grid_north_pole_latitude", grid_north_pole_latitude + ) + self.grid_north_pole_longitude = _validate_map_parameter( + "grid_north_pole_longitude", grid_north_pole_longitude + ) + self.north_pole_grid_longitude = _validate_map_parameter( + "north_pole_grid_longitude", north_pole_grid_longitude + ) diff --git a/cf/gridmappings/sinusoidal.py b/cf/gridmappings/sinusoidal.py new file mode 100644 index 0000000000..9da1f93151 --- /dev/null +++ b/cf/gridmappings/sinusoidal.py @@ -0,0 +1,67 @@ +from .abstract import GridMapping +from .abstract.gridmappingbase import _validate_map_parameter + + +class Sinusoidal(GridMapping): + """The Sinusoidal (Sanson-Flamsteed) grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#_sinusoidal + + or the corresponding PROJ projection page: + + https://proj.org/en/9.2/operations/projections/sinu.html + + for more information. + + .. versionadded:: GMVER + + :Parameters: + + longitude_of_central_meridian: number or scalar `Data`, optional + The longitude of (natural) origin i.e. central meridian. + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + """ + + grid_mapping_name = "sinusoidal" + proj_id = "sinu" + + def __init__( + self, + longitude_of_projection_origin=0.0, + false_easting=0.0, + false_northing=0.0, + **kwargs, + ): + super().__init__(**kwargs) + + self.longitude_of_projection_origin = _validate_map_parameter( + "longitude_of_projection_origin", longitude_of_projection_origin + ) + self.false_easting = _validate_map_parameter( + "false_easting", false_easting + ) + self.false_northing = _validate_map_parameter( + "false_northing", false_northing + ) diff --git a/cf/gridmappings/stereographic.py b/cf/gridmappings/stereographic.py new file mode 100644 index 0000000000..69fed28481 --- /dev/null +++ b/cf/gridmappings/stereographic.py @@ -0,0 +1,81 @@ +from .abstract import AzimuthalGridMapping +from .abstract.gridmappingbase import _validate_map_parameter + + +class Stereographic(AzimuthalGridMapping): + """The Stereographic grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#_stereographic + + or the corresponding PROJ projection page: + + https://proj.org/en/9.2/operations/projections/stere.html + + for more information. + + .. versionadded:: GMVER + + :Parameters: + + longitude_of_projection_origin: number or scalar `Data`, optional + The longitude of projection center (PROJ 'lon_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + latitude_of_projection_origin: number or scalar `Data`, optional + The latitude of projection center (PROJ 'lat_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_north', else the `Data` units are + taken and must be angular and compatible with latitude. + The default is 0.0 degrees_north. + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + scale_factor_at_projection_origin: number or scalar `Data`, optional + The scale factor used in the projection (PROJ 'k_0' value). + Unitless, so `Data` must be unitless. The default is 1.0. + + """ + + grid_mapping_name = "stereographic" + proj_id = "stere" + + def __init__( + self, + false_easting=0.0, + false_northing=0.0, + longitude_of_projection_origin=0.0, + latitude_of_projection_origin=0.0, + scale_factor_at_projection_origin=1.0, + **kwargs, + ): + super().__init__( + false_easting=false_easting, + false_northing=false_northing, + longitude_of_projection_origin=longitude_of_projection_origin, + latitude_of_projection_origin=latitude_of_projection_origin, + **kwargs, + ) + + self.scale_factor_at_projection_origin = _validate_map_parameter( + "scale_factor_at_projection_origin", + scale_factor_at_projection_origin, + ) diff --git a/cf/gridmappings/transversemercator.py b/cf/gridmappings/transversemercator.py new file mode 100644 index 0000000000..22a952d18b --- /dev/null +++ b/cf/gridmappings/transversemercator.py @@ -0,0 +1,85 @@ +from .abstract import CylindricalGridMapping +from .abstract.gridmappingbase import _validate_map_parameter + + +class TransverseMercator(CylindricalGridMapping): + """The Transverse Mercator grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#_transverse_mercator + + or the corresponding PROJ projection page: + + https://proj.org/en/9.2/operations/projections/tmerc.html + + for more information. + + .. versionadded:: GMVER + + :Parameters: + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + scale_factor_at_projection_origin: number or scalar `Data`, optional + The scale factor used in the projection (PROJ 'k_0' value). + Unitless, so `Data` must be unitless. The default is 1.0. + + longitude_of_central_meridian: number or scalar `Data`, optional + The longitude of (natural) origin i.e. central meridian. + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + latitude_of_projection_origin: number or scalar `Data`, optional + The latitude of projection center (PROJ 'lat_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_north', else the `Data` units are + taken and must be angular and compatible with latitude. + The default is 0.0 degrees_north. + + """ + + grid_mapping_name = "transverse_mercator" + proj_id = "tmerc" + + def __init__( + self, + scale_factor_at_central_meridian=1.0, + longitude_of_central_meridian=0.0, + latitude_of_projection_origin=0.0, + false_easting=0.0, + false_northing=0.0, + **kwargs, + ): + super().__init__( + false_easting=false_easting, + false_northing=false_northing, + **kwargs, + ) + + self.scale_factor_at_central_meridian = _validate_map_parameter( + "scale_factor_at_central_meridian", + scale_factor_at_central_meridian, + ) + self.longitude_of_central_meridian = _validate_map_parameter( + "longitude_of_central_meridian", longitude_of_central_meridian + ) + self.latitude_of_projection_origin = _validate_map_parameter( + "latitude_of_projection_origin", latitude_of_projection_origin + ) diff --git a/cf/gridmappings/verticalperspective.py b/cf/gridmappings/verticalperspective.py new file mode 100644 index 0000000000..0050d63b1e --- /dev/null +++ b/cf/gridmappings/verticalperspective.py @@ -0,0 +1,62 @@ +from .abstract import PerspectiveGridMapping + + +class VerticalPerspective(PerspectiveGridMapping): + """The Vertical (or Near-sided) Perspective grid mapping. + + See the CF Conventions document 'Appendix F: Grid Mappings' sub-section on + this grid mapping: + + http://cfconventions.org/Data/cf-conventions/cf-conventions-1.10/ + cf-conventions.html#vertical-perspective + + or the corresponding PROJ projection page: + + https://proj.org/en/9.2/operations/projections/nsper.html + + for more information. + + .. versionadded:: GMVER + + :Parameters: + + perspective_point_height: number or scalar `Data` + The height of the view point above the surface (PROJ + 'h') value, for example the height of a satellite above + the Earth, in units of meters 'm'. If provided + as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + + longitude_of_projection_origin: number or scalar `Data`, optional + The longitude of projection center (PROJ 'lon_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_east', else the `Data` units are + taken and must be angular and compatible with longitude. + The default is 0.0 degrees_east. + + latitude_of_projection_origin: number or scalar `Data`, optional + The latitude of projection center (PROJ 'lat_0' value). + If provided as a number or `Data` without units, the units + are taken as 'degrees_north', else the `Data` units are + taken and must be angular and compatible with latitude. + The default is 0.0 degrees_north. + + false_easting: number or scalar `Data`, optional + The false easting (PROJ 'x_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + false_northing: number or scalar `Data`, optional + The false northing (PROJ 'y_0') value. + If provided as a number or `Data` without units, the units + are taken as metres 'm', else the `Data` units are + taken and must be compatible with distance. + The default is 0.0 metres. + + """ + + grid_mapping_name = "vertical_perspective" + proj_id = "nsper" diff --git a/cf/mixin/fielddomain.py b/cf/mixin/fielddomain.py index 8af4ba4f8c..9471e9f84b 100644 --- a/cf/mixin/fielddomain.py +++ b/cf/mixin/fielddomain.py @@ -2611,6 +2611,57 @@ def refs(self, *identities, **filter_kwargs): """Alias for `coordinate_references`.""" return self.coordinate_references(*identities, **filter_kwargs) + @_inplace_enabled(default=False) + def create_2d_lats_and_lons( + self, destination_crs="LatitudeLongitude", inplace=False): + """Create 2-dimensional latitude and longitude coordinates. + + The generated coordinates are added to the `{{class}}` as new + auxiliary coordinates if the operation is in-place, else, by + default, a new copy of the `{{class}}` is returned with the + auxiliary coordinates included. + + .. versionadded:: GMVER + + :Parameters: + + destination_crs: `str` + The coordinate reference system that defines the + grid on which to create the latitudes and longitudes. + It should be an identifier which corresponds to the + name of a CF Grid Mapping class. By default, the class + taken is cf.LatitudeLongitude with "LatitudeLongitude". + + .. note:: Creating coordinates on a grid other than + the default cf.LatitudeLongitude is not + yet supported. + + {{inplace: `bool`, optional}} + + :Returns: + + `{{class}}` or `None` + The {{class}} with the newly-created latitude and + longitude coordinates incorporated as new auxiliary + coordinate constructs, or `None` if the operation + was in-place. + + **Examples** + + TODOGM + + """ + if destination_crs != "LatitudeLongitude": + raise NotImplementedError( + "Creating latitude and longitude coordinates for " + "a destination grid that is not described by " + "cf.LatitudeLongitude, the default, is not yet supported." + ) + + f = _inplace_enabled_define_and_cleanup(self) + + # TODO functional code here to go from input to lat_data and lon_dat + def _create_ancillary_mask_component(mask_shape, ind, compress): """Create an ancillary mask component. diff --git a/cf/test/test_Domain.py b/cf/test/test_Domain.py index b548c553a0..0d097a3ff5 100644 --- a/cf/test/test_Domain.py +++ b/cf/test/test_Domain.py @@ -295,6 +295,21 @@ def test_Domain_transpose(self): def test_Domain_size(self): self.assertEqual(self.d.size, 90) + def test_Domain_get_grid_mappings(self): + self.assertEqual( + self.d.get_grid_mappings(), { + "coordinatereference1": "rotated_latitude_longitude" + } + ) + self.assertEqual( + self.d.get_grid_mappings(as_class=True), { + "coordinatereference1": cf.RotatedLatitudeLongitude( + grid_north_pole_latitude=38.0, + grid_north_pole_longitude=190.0, + ) + } + ) + def test_Domain_create_regular(self): domain = cf.Domain.create_regular((-180, 180, 1), (-90, 90, 1)) self.assertIsInstance(domain, cf.Domain) diff --git a/cf/test/test_Field.py b/cf/test/test_Field.py index 9292d76056..0943e79eb6 100644 --- a/cf/test/test_Field.py +++ b/cf/test/test_Field.py @@ -2630,6 +2630,26 @@ def test_Field_auxiliary_to_dimension_to_auxiliary(self): with self.assertRaises(ValueError): f.auxiliary_to_dimension("latitude") + def test_Field_get_grid_mappings(self): + self.assertEqual( + self.f.get_grid_mappings(), + {"coordinatereference1": "rotated_latitude_longitude"} + ) + self.assertEqual( + self.f.get_grid_mappings(as_class=True), { + "coordinatereference1": cf.RotatedLatitudeLongitude( + grid_north_pole_latitude=38.0, + grid_north_pole_longitude=190.0, + ) + } + ) + self.assertEqual( + self.f0.get_grid_mappings(), {} + ) + self.assertEqual( + self.f0.get_grid_mappings(as_class=True), {} + ) + if __name__ == "__main__": print("Run date:", datetime.datetime.now()) diff --git a/cf/test/test_gridmappings.py b/cf/test/test_gridmappings.py new file mode 100644 index 0000000000..81cddd67d8 --- /dev/null +++ b/cf/test/test_gridmappings.py @@ -0,0 +1,356 @@ +import datetime +import faulthandler +import unittest + +import numpy as np + +import cf + +faulthandler.enable() # to debug seg faults and timeouts + +pyproj_imported = False +try: + import pyproj # noqa: F401 + + pyproj_imported = True +except ImportError: + pass + + +_all_abstract_grid_mappings = ( + cf.GridMapping, + cf.AzimuthalGridMapping, + cf.ConicGridMapping, + cf.CylindricalGridMapping, + cf.LatLonGridMapping, + cf.PerspectiveGridMapping, +) +_all_concrete_grid_mappings = ( + cf.AlbersEqualArea, + cf.AzimuthalEquidistant, + cf.Geostationary, + cf.LambertAzimuthalEqualArea, + cf.LambertConformalConic, + cf.LambertCylindricalEqualArea, + cf.LatitudeLongitude, + cf.Mercator, + cf.ObliqueMercator, + cf.Orthographic, + cf.PolarStereographic, + cf.RotatedLatitudeLongitude, + cf.Sinusoidal, + cf.Stereographic, + cf.TransverseMercator, + cf.VerticalPerspective, +) + + +# These are those of the above which have required positional arguments +all_grid_mappings_required_args = { + "AlbersEqualArea": {"standard_parallel": (0.0, None)}, + "ConicGridMapping": {"standard_parallel": (0.0, None)}, + "Geostationary": {"perspective_point_height": 1000}, + "LambertConformalConic": {"standard_parallel": (1.0, 1.0)}, + "PerspectiveGridMapping": {"perspective_point_height": 1000}, + "RotatedLatitudeLongitude": { + "grid_north_pole_latitude": 0.0, + "grid_north_pole_longitude": 0.0, + }, + "VerticalPerspective": {"perspective_point_height": 1000}, +} + + +class GridMappingsTest(unittest.TestCase): + """Unit test for the GridMapping class and any derived classes.""" + + # Of the example fields, only 1, 6 and 7 have any coordinate references + # with a coordinate conversion, hence use these to test, plus 0 as an + # example case of a field without a coordinate reference at all. + f = cf.example_fields() + # TODOGM, could do with a new example field or two with a grid mapping + # other than those two below, for diversity and testing on etc. + f0 = f[0] # No coordinate reference + f1 = f[1] # 2, with grid mappings of [None, 'rotated_latitude_longitude'] + f6 = f[6] # 1, with grid mapping of ['latitude_longitude'] + f7 = f[7] # 1, with grid mapping of ['rotated_latitude_longitude'] + f_with_gm = (f1, f6, f7) + f_wth_gm_mapping = { + "f1": cf.RotatedLatitudeLongitude, + "f6": cf.LatitudeLongitude, + "f7": cf.RotatedLatitudeLongitude, + } + + @unittest.skipUnless(pyproj_imported, "Requires pyproj package.") + def test_grid_mapping__init__(self): + """Test GridMapping object initiation.""" + for cls in _all_concrete_grid_mappings: + if cls.__name__ not in all_grid_mappings_required_args: + g = cls() + g.grid_mapping_name + else: + example_minimal_args = all_grid_mappings_required_args[ + cls.__name__ + ] + g = cls(**example_minimal_args) + g.grid_mapping_name + + for cls in _all_abstract_grid_mappings: + if cls.__name__ not in all_grid_mappings_required_args: + g = cls() + self.assertEqual(g.grid_mapping_name, None) + else: + example_minimal_args = all_grid_mappings_required_args[ + cls.__name__ + ] + g = cls(**example_minimal_args) + self.assertEqual(g.grid_mapping_name, None) + + @unittest.skipUnless(pyproj_imported, "Requires pyproj package.") + def test_grid_mapping__repr__str__(self): + """Test all means of GridMapping inspection.""" + for cls in _all_concrete_grid_mappings: + if cls.__name__ not in all_grid_mappings_required_args: + g = cls() + else: + example_minimal_args = all_grid_mappings_required_args[ + cls.__name__ + ] + g = cls(**example_minimal_args) + repr(g) + str(g) + + g1 = cf.Mercator() + self.assertEqual(repr(g1), "") + self.assertEqual( + str(g1), "" + ) + + g2 = cf.Orthographic() + self.assertEqual(repr(g2), "") + self.assertEqual( + str(g2), "" + ) + + g3 = cf.Sinusoidal() + self.assertEqual(repr(g3), "") + self.assertEqual(str(g3), "") + + g4 = cf.Stereographic() + self.assertEqual(repr(g4), "") + self.assertEqual( + str(g4), "" + ) + + @unittest.skipUnless(pyproj_imported, "Requires pyproj package.") + def test_grid_mapping_is_latlon_gm(self): + """Test the 'is_latlon_gm' method on all GridMappings.""" + # In this one case we expect True... + g = cf.LatitudeLongitude + self.assertTrue(g.is_latlon_gm()) # check on class + self.assertTrue(g().is_latlon_gm()) # check on instance + + # ...and expect False for all other GridMappings + for cls in _all_concrete_grid_mappings: + if not issubclass(cls, cf.LatitudeLongitude): + self.assertFalse(cls.is_latlon_gm()) + + @unittest.skipUnless(pyproj_imported, "Requires pyproj package.") + def test_grid_mapping_GM_class(self): + """Test factory creation with the GM class.""" + # Note f1 has cr0 with no GM and cr1 with 'rotated_latitude_longitude' + cr0 = self.f1.coordinate_references('coordinatereference0').value() + cr1 = self.f1.coordinate_references('coordinatereference1').value() + # Note f6 has just one CR, grid mapping name of ['latitude_longitude'] + cr2 = self.f6.coordinate_references('coordinatereference0').value() + + r = cf.GM(cr1) + self.assertTrue(isinstance(r, cf.RotatedLatitudeLongitude)) + self.assertEqual(r.grid_mapping_name, "rotated_latitude_longitude") + self.assertEqual(r.grid_north_pole_latitude, 38.0) + self.assertEqual(r.grid_north_pole_longitude, 190.0) + + l = cf.GM(cr2) + self.assertTrue(isinstance(l, cf.LatitudeLongitude)) + self.assertEqual(l.grid_mapping_name, "latitude_longitude") + + with self.assertRaises(cf.InvalidGridMapping): + cf.GM(cr0) + + # Test creation with overriding of parameters from the CR + r = cf.GM(cr1, **{"grid_north_pole_latitude": -45}) + self.assertTrue(isinstance(r, cf.RotatedLatitudeLongitude)) + self.assertEqual(r.grid_mapping_name, "rotated_latitude_longitude") + self.assertEqual(r.grid_north_pole_latitude, -45) + self.assertEqual(r.grid_north_pole_longitude, 190.0) + + # TODOGM extend greatly - by creating CRs with more varied + # 'grid_mapping_name' attributes with corresponding valid and invalid + # parameters... + + @unittest.skipUnless(pyproj_imported, "Requires pyproj package.") + def test_grid_mapping_map_parameter_validation(self): + """Test the validation of map parameters to Grid Mapping classes.""" + g1 = cf.Mercator( + false_easting=10.0, + false_northing=cf.Data(-20, units="cm"), + standard_parallel=(None, 50), + longitude_of_projection_origin=cf.Data( + -40.0, units="degrees_east" + ), + scale_factor_at_projection_origin=3.0, + prime_meridian_name="brussels", + ) + self.assertEqual(g1.false_easting, cf.Data(10.0, "m")) + self.assertEqual(g1.false_northing, cf.Data(-0.2, "m")) + self.assertEqual( + g1.standard_parallel, (None, cf.Data(50, "degrees_north")) + ) + self.assertEqual( + g1.longitude_of_projection_origin, cf.Data(-40.0, "degrees_east") + ) + self.assertEqual(g1.scale_factor_at_projection_origin, cf.Data(3.0, 1)) + self.assertEqual(g1.prime_meridian_name, "brussels") + + # TODOGM extend this test a lot with testing like the above and + # with systematic coverage over valid inputs + + @unittest.skipUnless(pyproj_imported, "Requires pyproj package.") + def test_grid_mapping__get_cf_grid_mapping_from_name(self): + """Test the '_get_cf_grid_mapping_from_name' function.""" + for gm_name, cf_gm_class in { + "vertical_perspective": cf.VerticalPerspective, + "oblique_mercator": cf.ObliqueMercator, + "albers_conical_equal_area": cf.AlbersEqualArea, + "lambert_conformal_conic": cf.LambertConformalConic, + "some_unsupported_name": None, + }.items(): + pass # TODO UPDATE with class + # self.assertEqual( + # cf._get_cf_grid_mapping_from_name(gm_name), cf_gm_class + # ) + + def test_grid_mapping_convert_proj_angular_data_to_cf(self): + """Test the 'convert_proj_angular_data_to_cf' function.""" + + # Check representative valid inputs + for input_with_correct_output in [ + # Check float value and no suffix + (("45.0", None), cf.Data(45.0, units="degrees")), + (("45.0", "lat"), cf.Data(45.0, units="degrees_north")), + (("45.0", "lon"), cf.Data(45.0, units="degrees_east")), + # Check integer and no suffix + (("100", None), cf.Data(100, units="degrees")), + (("100", "lat"), cf.Data(100, units="degrees_north")), + (("100", "lon"), cf.Data(100, units="degrees_east")), + # Check >360 degrees (over a full revolution) and "r" suffix + ((f"{3.0 * np.pi}r", None), cf.Data(3.0 * np.pi, units="radians")), + ( + (f"{3.0 * np.pi}r", "lat"), + cf.Data(540.0, units="degrees_north"), + ), + ((f"{3.0 * np.pi}r", "lon"), cf.Data(540.0, units="degrees_east")), + # Check "R" suffix + ((f"{0.5 * np.pi}R", None), cf.Data(0.5 * np.pi, units="radians")), + ((f"{0.5 * np.pi}R", "lat"), cf.Data(90.0, units="degrees_north")), + ((f"{0.5 * np.pi}R", "lon"), cf.Data(90.0, units="degrees_east")), + # Check integer value and "d" suffix + (("10d", None), cf.Data(10, units="degrees")), + (("10d", "lat"), cf.Data(10, units="degrees_north")), + (("10d", "lon"), cf.Data(10, units="degrees_east")), + # Check >180 float and "D" suffix + (("200.123D", None), cf.Data(200.123, units="degrees")), + (("200.123D", "lat"), cf.Data(200.123, units="degrees_north")), + (("200.123D", "lon"), cf.Data(200.123, units="degrees_east")), + # Check negative numeric value and "°" suffix + (("-70.5°", None), cf.Data(-70.5, units="degrees")), + (("-70.5°", "lat"), cf.Data(-70.5, units="degrees_north")), + (("-70.5°", "lon"), cf.Data(-70.5, units="degrees_east")), + # Check zero and lack of digits after point edge cases + (("0", None), cf.Data(0, units="degrees")), + (("0.0", "lat"), cf.Data(0.0, units="degrees_north")), + (("-0.", "lon"), cf.Data(0.0, units="degrees_east")), + ]: + _input, correct_output = input_with_correct_output + d = cf.convert_proj_angular_data_to_cf(*_input) + self.assertTrue(d.equals(correct_output, verbose=2)) + + def test_grid_mapping_convert_cf_angular_data_to_proj(self): + """Test the 'convert_cf_angular_data_to_proj' function.""" + # Check representative valid inputs + for input_with_correct_output in [ + # Check float and basic lat/lon-context free degree unit + (cf.Data(45.0, units="degrees"), "45.0"), + # Check integer and various units possible for lat/lon + (cf.Data(45, units="degrees_north"), "45"), + (cf.Data(45, units="degrees_N"), "45"), + (cf.Data(45, units="degreeN"), "45"), + (cf.Data(45, units="degrees_east"), "45"), + (cf.Data(45, units="degrees_E"), "45"), + (cf.Data(45, units="degreeE"), "45"), + (cf.Data(45, units="degree"), "45"), + # Check negative + (cf.Data(-0.1, units="degrees"), "-0.1"), + (cf.Data(-10, units="degrees"), "-10"), + # Check zero case + (cf.Data(0, units="degrees_north"), "0"), + (cf.Data(0.0, units="degrees_north"), "0.0"), + # Check radians units cases and >180 + (cf.Data(190, units="radians"), "190R"), + (cf.Data(190.0, units="radians"), "190.0R"), + # Check flot with superfluous 0 + (cf.Data(120.100, units="degrees"), "120.1"), + ]: + _input, correct_output = input_with_correct_output + p = cf.convert_cf_angular_data_to_proj(_input) + self.assertEqual(p, correct_output) + + # Check representative invalid inputs error correctly + for bad_input in [ + cf.Data([1, 2, 3]), # not singular (size 1) + cf.Data(45), # no units + cf.Data(45, "m"), # non-angular units + cf.Data(2, "elephants"), # bad/non-CF units + ]: + with self.assertRaises(ValueError): + cf.convert_cf_angular_data_to_proj(bad_input) + with self.assertRaises(TypeError): + # Non-numeric value + cf.convert_cf_angular_data_to_proj(cf.Data("N", "radians")) + + # Note that 'convert_cf_angular_data_to_proj' and + # 'convert_proj_angular_data_to_cf' are not strict inverse + # functions, since the former will convert to the *simplest* + # way to specify the PROJ input, namely with no suffix for + # degrees(_X) units and the 'R' suffix for radians, whereas + # the input might have 'D' or 'r' etc. instead. + # + # However, we check that inputs that are expected to be + # undone to their original form by operation of both + # functions, namely those with this 'simplest' PROJ form, + # do indeed get re-generated with operation of both: + for p in ("10", "-10", "10.11", "0", "1R", "0.2R", "-0.1R", "0R"): + p2 = cf.convert_cf_angular_data_to_proj( + cf.convert_proj_angular_data_to_cf(p) + ) + self.assertEqual(p, p2) + + # With a lat or lon 'context'. Only non-radians inputs will + # be re-generated since degrees_X gets converted back to the + # default degrees, so skip those in these test cases. + if not p.endswith("R"): + p3 = cf.convert_cf_angular_data_to_proj( + cf.convert_proj_angular_data_to_cf(p, context="lat") + ) + self.assertEqual(p, p3) + + p4 = cf.convert_cf_angular_data_to_proj( + cf.convert_proj_angular_data_to_cf(p, context="lon") + ) + self.assertEqual(p, p4) + + +if __name__ == "__main__": + print("Run date:", datetime.datetime.now()) + cf.environment() + print() + unittest.main(verbosity=2) diff --git a/setup.py b/setup.py index 597f938d40..ca0401725f 100755 --- a/setup.py +++ b/setup.py @@ -317,6 +317,8 @@ def compile(): "cf.regrid", "cf.umread_lib", "cf.test", + "cf.gridmappings", + "cf.gridmappings.abstract", ], package_data={"cf": package_data}, scripts=["scripts/cfa"],