# Copyright Iris contributors
#
# This file is part of Iris and is released under the BSD license.
# See LICENSE in the root of the repository for full licensing details.
"""Various utilities and numeric transformations relevant to cartography."""
from collections import namedtuple
import copy
import warnings
import cartopy.crs as ccrs
import cartopy.img_transform
import cf_units
import dask.array as da
import numpy as np
import numpy.ma as ma
import iris.coord_systems
import iris.coords
import iris.exceptions
from iris.util import _meshgrid
import iris.warnings
from ._grid_angles import gridcell_angles, rotate_grid_vectors
# List of contents to control Sphinx autodocs.
# Unfortunately essential to get docs for the grid_angles functions.
__all__ = [
"DistanceDifferential",
"PartialDifferential",
"area_weights",
"cosine_latitude_weights",
"get_xy_contiguous_bounded_grids",
"get_xy_grids",
"gridcell_angles",
"project",
"rotate_grid_vectors",
"rotate_pole",
"rotate_winds",
"unrotate_pole",
"wrap_lons",
]
# This value is used as a fall-back if the cube does not define the earth
DEFAULT_SPHERICAL_EARTH_RADIUS = 6367470
# TODO: This should not be necessary, as CF is always in meters
DEFAULT_SPHERICAL_EARTH_RADIUS_UNIT = cf_units.Unit("m")
# Distance differentials for coordinate systems at specified locations
DistanceDifferential = namedtuple("DistanceDifferential", "dx1 dy1 dx2 dy2")
# Partial differentials between coordinate systems
PartialDifferential = namedtuple("PartialDifferential", "dx1 dy1")
[docs]
def wrap_lons(lons, base, period):
"""Wrap longitude values into the range between base and base+period.
Parameters
----------
lons :
base :
period :
Examples
--------
.. testsetup::
import numpy as np
from iris.analysis.cartography import wrap_lons
::
>>> print(wrap_lons(np.array([185, 30, -200, 75]), -180, 360))
[-175. 30. 160. 75.]
Notes
-----
This function maintains laziness when called; it does not realise data.
See more at :doc:`/userguide/real_and_lazy_data`.
"""
# It is important to use 64bit floating precision when changing a floats
# numbers range.
lons = lons.astype(np.float64)
return ((lons - base) % period) + base
[docs]
def unrotate_pole(rotated_lons, rotated_lats, pole_lon, pole_lat):
"""Convert rotated-pole to unrotated longitudes and latitudes.
``pole_lat`` should describe the location of the rotated pole that
describes the arrays of rotated-pole longitudes and latitudes.
Convert arrays of rotated-pole longitudes and latitudes to unrotated
arrays of longitudes and latitudes. The values of ``pole_lon`` and
``pole_lat`` should describe the location of the rotated pole that
describes the arrays of rotated-pole longitudes and latitudes.
As the arrays of rotated-pole longitudes and latitudes must describe a
rectilinear grid, the arrays of rotated-pole longitudes and latitudes must
be of the same shape as each other.
.. note:: Uses proj.4 to perform the conversion.
Parameters
----------
rotated_lons :
An array of rotated-pole longitude values.
rotated_lats :
An array of rotated-pole latitude values.
pole_lon :
The longitude of the rotated pole that describes the arrays of
rotated-pole longitudes and latitudes.
pole_lat :
The latitude of the rotated pole that describes the arrays of
rotated-pole longitudes and latitudes.
Returns
-------
An array of unrotated longitudes and an array of unrotated latitudes.
Examples
--------
::
lons, lats = unrotate_pole(rotated_lons, rotated_lats, pole_lon, pole_lat)
"""
src_proj = ccrs.RotatedGeodetic(pole_longitude=pole_lon, pole_latitude=pole_lat)
target_proj = ccrs.Geodetic()
res = target_proj.transform_points(x=rotated_lons, y=rotated_lats, src_crs=src_proj)
unrotated_lon = res[..., 0]
unrotated_lat = res[..., 1]
return unrotated_lon, unrotated_lat
[docs]
def rotate_pole(lons, lats, pole_lon, pole_lat):
"""Convert unrotated longitudes and latitudes to rotated-pole.
The values of ``pole_lon`` and ``pole_lat``
should describe the rotated pole that the arrays of longitudes and
latitudes are to be rotated onto.
As the arrays of longitudes and latitudes must describe a rectilinear grid,
the arrays of rotated-pole longitudes and latitudes must be of the same
shape as each other.
.. note:: Uses proj.4 to perform the conversion.
Parameters
----------
lons :
An array of longitude values.
lats :
An array of latitude values.
pole_lon :
The longitude of the rotated pole that the arrays of longitudes and
latitudes are to be rotated onto.
pole_lat :
The latitude of the rotated pole that the arrays of longitudes and
latitudes are to be rotated onto.
Returns
-------
An array of rotated-pole longitudes and an array of rotated-pole latitudes.
Examples
--------
::
rotated_lons, rotated_lats = rotate_pole(lons, lats, pole_lon, pole_lat)
"""
src_proj = ccrs.Geodetic()
target_proj = ccrs.RotatedGeodetic(pole_longitude=pole_lon, pole_latitude=pole_lat)
res = target_proj.transform_points(x=lons, y=lats, src_crs=src_proj)
rotated_lon = res[..., 0]
rotated_lat = res[..., 1]
return rotated_lon, rotated_lat
def _get_lon_lat_coords(cube):
def search_for_coord(coord_iterable, coord_name):
return [coord for coord in coord_iterable if coord_name in coord.name()]
lat_coords = search_for_coord(cube.dim_coords, "latitude") or search_for_coord(
cube.coords(), "latitude"
)
lon_coords = search_for_coord(cube.dim_coords, "longitude") or search_for_coord(
cube.coords(), "longitude"
)
if len(lat_coords) > 1 or len(lon_coords) > 1:
raise ValueError(
"Calling `_get_lon_lat_coords` with multiple same-type (i.e. dim/aux) lat or lon coords"
" is currently disallowed"
)
lat_coord = lat_coords[0]
lon_coord = lon_coords[0]
return lon_coord, lat_coord
def _xy_range(cube, mode=None):
"""Return the x & y range of this Cube.
Parameters
----------
cube :
The cube for which to calculate xy extents.
mode : optional
If the coordinate has bounds, set this to specify the
min/max calculation.
Set to iris.coords.POINT_MODE or iris.coords.BOUND_MODE.
"""
# Helpful error if we have an inappropriate CoordSystem
cs = cube.coord_system("CoordSystem")
cs_valid_types = (
iris.coord_systems.GeogCS,
iris.coord_systems.RotatedGeogCS,
)
if (cs is not None) and not isinstance(cs, cs_valid_types):
raise ValueError("Latlon coords cannot be found with {0}.".format(type(cs)))
x_coord, y_coord = cube.coord(axis="X"), cube.coord(axis="Y")
cs = cube.coord_system("CoordSystem")
if x_coord.has_bounds() != y_coord.has_bounds():
raise ValueError(
"Cannot get the range of the x and y coordinates if they do "
"not have the same presence of bounds."
)
if x_coord.has_bounds():
if mode not in [iris.coords.POINT_MODE, iris.coords.BOUND_MODE]:
raise ValueError('When the coordinate has bounds, please specify "mode".')
_mode = mode
else:
_mode = iris.coords.POINT_MODE
# Get the x and y grids
if isinstance(cs, iris.coord_systems.RotatedGeogCS):
if _mode == iris.coords.POINT_MODE:
x, y = get_xy_grids(cube)
else:
x, y = get_xy_contiguous_bounded_grids(cube)
else:
if _mode == iris.coords.POINT_MODE:
x = x_coord.points
y = y_coord.points
else:
x = x_coord.bounds
y = y_coord.bounds
# Get the x and y range
if getattr(x_coord, "circular", False):
x_range = (np.min(x), np.min(x) + x_coord.units.modulus)
else:
x_range = (np.min(x), np.max(x))
y_range = (np.min(y), np.max(y))
return (x_range, y_range)
[docs]
def get_xy_grids(cube):
"""Return 2D X and Y points for a given cube.
Parameters
----------
cube :
The cube for which to generate 2D X and Y points.
Examples
--------
::
x, y = get_xy_grids(cube)
Notes
-----
This function maintains laziness when called; it does not realise data.
See more at :doc:`/userguide/real_and_lazy_data`.
"""
x_coord, y_coord = cube.coord(axis="X"), cube.coord(axis="Y")
x = x_coord.points
y = y_coord.points
if x.ndim == y.ndim == 1:
# Convert to 2D.
x, y = _meshgrid(x, y)
elif x.ndim == y.ndim == 2:
# They are already in the correct shape.
pass
else:
raise ValueError("Expected 1D or 2D XY coords")
return (x, y)
[docs]
def get_xy_contiguous_bounded_grids(cube):
"""Return 2d arrays for x and y bounds.
Returns array of shape (n+1, m+1).
Parameters
----------
cube : :class:`iris.cube.Cube`
Examples
--------
::
xs, ys = get_xy_contiguous_bounded_grids(cube)
Notes
-----
This function maintains laziness when called; it does not realise data.
See more at :doc:`/userguide/real_and_lazy_data`.
"""
x_coord, y_coord = cube.coord(axis="X"), cube.coord(axis="Y")
x = x_coord.contiguous_bounds()
y = y_coord.contiguous_bounds()
x, y = _meshgrid(x, y)
return (x, y)
def _quadrant_area(radian_lat_bounds, radian_lon_bounds, radius_of_earth):
"""Calculate spherical segment areas.
Area weights are calculated for each lat/lon cell as:
.. math::
r^2 (lon_1 - lon_0) ( sin(lat_1) - sin(lat_0))
The resulting array will have a shape of
*(radian_lat_bounds.shape[0], radian_lon_bounds.shape[0])*
The calculations are done at 64 bit precision and the returned array
will be of type numpy.float64.
Parameters
----------
radian_lat_bounds :
[n,2] array of latitude bounds (radians).
radian_lon_bounds :
[n,2] array of longitude bounds (radians).
radius_of_earth :
Radius of the earth (currently assumed spherical).
"""
# ensure pairs of bounds
if (
radian_lat_bounds.shape[-1] != 2
or radian_lon_bounds.shape[-1] != 2
or radian_lat_bounds.ndim != 2
or radian_lon_bounds.ndim != 2
):
raise ValueError("Bounds must be [n,2] array")
# fill in a new array of areas
radius_sqr = radius_of_earth**2
radian_lat_64 = radian_lat_bounds.astype(np.float64)
radian_lon_64 = radian_lon_bounds.astype(np.float64)
ylen = np.sin(radian_lat_64[:, 1]) - np.sin(radian_lat_64[:, 0])
xlen = radian_lon_64[:, 1] - radian_lon_64[:, 0]
areas = radius_sqr * np.outer(ylen, xlen)
# we use abs because backwards bounds (min > max) give negative areas.
return np.abs(areas)
[docs]
def area_weights(cube, normalize=False, compute=True, chunks=None):
r"""Return an array of area weights, with the same dimensions as the cube.
This is a 2D lat/lon area weights array, repeated over the non lat/lon
dimensions.
The cube must have coordinates 'latitude' and 'longitude' with bounds.
Area weights are calculated for each lat/lon cell as:
.. math::
r^2 (lon_1 - lon_0) (\sin(lat_1) - \sin(lat_0))
Currently, only supports a spherical datum.
Uses earth radius from the cube, if present and spherical.
Defaults to iris.analysis.cartography.DEFAULT_SPHERICAL_EARTH_RADIUS.
Parameters
----------
cube : :class:`iris.cube.Cube`
The cube to calculate area weights for.
normalize : bool, default=False
If False, weights are grid cell areas. If True, weights are grid
cell areas divided by the total grid area.
compute : bool, default=True
If False, return a lazy dask array. If True, return a numpy array.
chunks : tuple, optional
If compute is False and a value is provided, then the result will use
these chunks instead of the same chunks as the cube data. The values
provided here will only be used along dimensions that are not latitude
or longitude.
Returns
-------
broad_weights :
"""
# Get the radius of the earth
cs = cube.coord_system("CoordSystem")
if isinstance(cs, iris.coord_systems.GeogCS):
if cs.inverse_flattening != 0.0:
warnings.warn(
"Assuming spherical earth from ellipsoid.",
category=iris.warnings.IrisDefaultingWarning,
)
radius_of_earth = cs.semi_major_axis
elif isinstance(cs, iris.coord_systems.RotatedGeogCS) and (
cs.ellipsoid is not None
):
if cs.ellipsoid.inverse_flattening != 0.0:
warnings.warn(
"Assuming spherical earth from ellipsoid.",
category=iris.warnings.IrisDefaultingWarning,
)
radius_of_earth = cs.ellipsoid.semi_major_axis
else:
warnings.warn(
"Using DEFAULT_SPHERICAL_EARTH_RADIUS.",
category=iris.warnings.IrisDefaultingWarning,
)
radius_of_earth = DEFAULT_SPHERICAL_EARTH_RADIUS
# Get the lon and lat coords and axes
try:
lon, lat = _get_lon_lat_coords(cube)
except IndexError:
raise ValueError(
"Cannot get latitude/longitude coordinates from cube {!r}.".format(
cube.name()
)
)
if lat.ndim > 1:
raise iris.exceptions.CoordinateMultiDimError(lat)
if lon.ndim > 1:
raise iris.exceptions.CoordinateMultiDimError(lon)
lat_dim = cube.coord_dims(lat)
lat_dim = lat_dim[0] if lat_dim else None
lon_dim = cube.coord_dims(lon)
lon_dim = lon_dim[0] if lon_dim else None
if not (lat.has_bounds() and lon.has_bounds()):
msg = (
"Coordinates {!r} and {!r} must have bounds to determine "
"the area weights.".format(lat.name(), lon.name())
)
raise ValueError(msg)
# Convert from degrees to radians
lat = lat.copy()
lon = lon.copy()
for coord in (lat, lon):
if coord.units in (cf_units.Unit("degrees"), cf_units.Unit("radians")):
coord.convert_units("radians")
else:
msg = (
"Units of degrees or radians required, coordinate "
"{!r} has units: {!r}".format(coord.name(), coord.units.name)
)
raise ValueError(msg)
# Create 2D weights from bounds.
# Use the geographical area as the weight for each cell
if compute:
lat_bounds = lat.bounds
lon_bounds = lon.bounds
else:
lat_bounds = lat.lazy_bounds()
lon_bounds = lon.lazy_bounds()
ll_weights = _quadrant_area(lat_bounds, lon_bounds, radius_of_earth)
# Normalize the weights if necessary.
if normalize:
ll_weights /= ll_weights.sum()
# Now we create an array of weights for each cell. This process will
# handle adding the required extra dimensions and also take care of
# the order of dimensions.
broadcast_dims = [x for x in (lat_dim, lon_dim) if x is not None]
wshape = []
for idim, dim in zip((0, 1), (lat_dim, lon_dim)):
if dim is not None:
wshape.append(ll_weights.shape[idim])
ll_weights = ll_weights.reshape(wshape)
broad_weights = iris.util.broadcast_to_shape(
ll_weights, cube.shape, broadcast_dims, chunks=chunks
)
return broad_weights
[docs]
def cosine_latitude_weights(cube):
r"""Calculate cosine latitude weights, with the same dimensions as the cube.
Return an array of latitude weights, with the same dimensions as
the cube. The weights are the cosine of latitude.
These are n-dimensional latitude weights repeated over the dimensions
not covered by the latitude coordinate.
The cube must have a coordinate with 'latitude' in the name. Out of
range values (greater than 90 degrees or less than -90 degrees) will
be clipped to the valid range.
Weights are calculated for each latitude as:
.. math::
w_l = \cos \phi_l
Parameters
----------
cube : :class:`iris.cube.Cube`
Examples
--------
Compute weights suitable for averaging type operations::
from iris.analysis.cartography import cosine_latitude_weights
cube = iris.load_cube(iris.sample_data_path('air_temp.pp'))
weights = cosine_latitude_weights(cube)
Compute weights suitable for EOF analysis (or other covariance type
analyses)::
import numpy as np
from iris.analysis.cartography import cosine_latitude_weights
cube = iris.load_cube(iris.sample_data_path('air_temp.pp'))
weights = np.sqrt(cosine_latitude_weights(cube))
Notes
-----
This function maintains laziness when called; it does not realise data.
See more at :doc:`/userguide/real_and_lazy_data`.
"""
# Find all latitude coordinates, we want one and only one.
lat_coords = [coord for coord in cube.coords() if "latitude" in coord.name()]
if len(lat_coords) > 1:
raise ValueError("Multiple latitude coords are currently disallowed.")
try:
lat = lat_coords[0]
except IndexError:
raise ValueError(
"Cannot get latitude coordinate from cube {!r}.".format(cube.name())
)
# Get the dimension position(s) of the latitude coordinate.
lat_dims = cube.coord_dims(lat)
# Convert to radians.
lat = lat.copy()
lat.convert_units("radians")
# Compute the weights as the cosine of latitude. In some cases,
# particularly when working in 32-bit precision, the latitude values can
# extend beyond the allowable range of [-pi/2, pi/2] due to numerical
# precision. We first check for genuinely out of range values, and issue a
# warning if these are found. Then the cosine is computed and clipped to
# the valid range [0, 1].
threshold = np.deg2rad(0.001) # small value for grid resolution
if np.any(lat.points < -np.pi / 2.0 - threshold) or np.any(
lat.points > np.pi / 2.0 + threshold
):
warnings.warn(
"Out of range latitude values will be clipped to the valid range.",
category=iris.warnings.IrisDefaultingWarning,
)
points = lat.points
l_weights = np.cos(points).clip(0.0, 1.0)
# Create weights for each grid point. This operation handles adding extra
# dimensions and also the order of the dimensions.
broadcast_dims = [x for x in lat_dims if x is not None]
wshape = []
for idim, dim in enumerate(lat_dims):
if dim is not None:
wshape.append(l_weights.shape[idim])
l_weights = l_weights.reshape(wshape)
broad_weights = iris.util.broadcast_to_shape(l_weights, cube.shape, broadcast_dims)
return broad_weights
[docs]
def project(cube, target_proj, nx=None, ny=None):
"""Nearest neighbour regrid to a specified target projection.
Return a new cube that is the result of projecting a cube with 1 or 2
dimensional latitude-longitude coordinates from its coordinate system into
a specified projection e.g. Robinson or Polar Stereographic.
This function is intended to be used in cases where the cube's coordinates
prevent one from directly visualising the data, e.g. when the longitude
and latitude are two dimensional and do not make up a regular grid.
Parameters
----------
cube : :class:`iris.cube.Cube`
An instance of :class:`iris.cube.Cube`.
target_proj : :class:`iris.coord_systems.CoordSystem`
An instance of the Cartopy Projection class, or an instance of
:class:`iris.coord_systems.CoordSystem` from which a projection
will be obtained.
nx : optional
Desired number of sample points in the x direction for a domain
covering the globe.
ny : optional
Desired number of sample points in the y direction for a domain
covering the globe.
Returns
-------
:class:`iris.cube.Cube`
An instance of :class:`iris.cube.Cube` and a list describing the
extent of the projection.
Notes
-----
.. note::
If there are both dim and aux latitude-longitude coordinates, only
the dim coordinates will be used.
.. note::
This function assumes global data and will if necessary extrapolate
beyond the geographical extent of the source cube using a nearest
neighbour approach. nx and ny then include those points which are
outside of the target projection.
.. note::
Masked arrays are handled by passing their masked status to the
resulting nearest neighbour values. If masked, the value in the
resulting cube is set to 0.
.. note::
This function does not maintain laziness when called; it realises data.
See more at :doc:`/userguide/real_and_lazy_data`.
.. warning::
This function uses a nearest neighbour approach rather than any form
of linear/non-linear interpolation to determine the data value of each
cell in the resulting cube. Consequently it may have an adverse effect
on the statistics of the data e.g. the mean and standard deviation
will not be preserved.
.. warning::
If the target projection is non-rectangular, e.g. Robinson, the target
grid may include points outside the boundary of the projection. The
latitude/longitude of such points may be unpredictable.
"""
try:
lon_coord, lat_coord = _get_lon_lat_coords(cube)
except IndexError:
raise ValueError(
"Cannot get latitude/longitude coordinates from cube {!r}.".format(
cube.name()
)
)
if lat_coord.coord_system != lon_coord.coord_system:
raise ValueError(
"latitude and longitude coords appear to have "
"different coordinates systems."
)
if lon_coord.units != "degrees":
lon_coord = lon_coord.copy()
lon_coord.convert_units("degrees")
if lat_coord.units != "degrees":
lat_coord = lat_coord.copy()
lat_coord.convert_units("degrees")
# Determine source coordinate system
if lat_coord.coord_system is None:
# Assume WGS84 latlon if unspecified
warnings.warn(
"Coordinate system of latitude and longitude "
"coordinates is not specified. Assuming WGS84 Geodetic.",
category=iris.warnings.IrisDefaultingWarning,
)
orig_cs = iris.coord_systems.GeogCS(
semi_major_axis=6378137.0, inverse_flattening=298.257223563
)
else:
orig_cs = lat_coord.coord_system
# Convert to cartopy crs
source_cs = orig_cs.as_cartopy_crs()
# Obtain coordinate arrays (ignoring bounds) and convert to 2d
# if not already.
source_x = lon_coord.points
source_y = lat_coord.points
if source_x.ndim != 2 or source_y.ndim != 2:
source_x, source_y = _meshgrid(source_x, source_y)
# Calculate target grid
target_cs = None
if isinstance(target_proj, iris.coord_systems.CoordSystem):
target_cs = target_proj
target_proj = target_proj.as_cartopy_projection()
# Resolution of new grid
if nx is None:
nx = source_x.shape[1]
if ny is None:
ny = source_x.shape[0]
target_x, target_y, extent = cartopy.img_transform.mesh_projection(
target_proj, nx, ny
)
# Determine dimension mappings - expect either 1d or 2d
if lat_coord.ndim != lon_coord.ndim:
raise ValueError(
"The latitude and longitude coordinates have different dimensionality."
)
latlon_ndim = lat_coord.ndim
lon_dims = cube.coord_dims(lon_coord)
lat_dims = cube.coord_dims(lat_coord)
if latlon_ndim == 1:
xdim = lon_dims[0]
ydim = lat_dims[0]
elif latlon_ndim == 2:
if lon_dims != lat_dims:
raise ValueError(
"The 2d latitude and longitude coordinates "
"correspond to different dimensions."
)
# If coords are 2d assume that grid is ordered such that x corresponds
# to the last dimension (shortest stride).
xdim = lon_dims[1]
ydim = lon_dims[0]
else:
raise ValueError(
"Expected the latitude and longitude coordinates "
"to have 1 or 2 dimensions, got {} and "
"{}.".format(lat_coord.ndim, lon_coord.ndim)
)
# Create array to store regridded data
new_shape = list(cube.shape)
new_shape[xdim] = nx
new_shape[ydim] = ny
new_data = ma.zeros(new_shape, cube.data.dtype)
# Create iterators to step through cube data in lat long slices
new_shape[xdim] = 1
new_shape[ydim] = 1
index_it = np.ndindex(*new_shape)
if lat_coord.ndim == 1 and lon_coord.ndim == 1:
slice_it = cube.slices([lat_coord, lon_coord])
elif lat_coord.ndim == 2 and lon_coord.ndim == 2:
slice_it = cube.slices(lat_coord)
else:
raise ValueError(
"Expected the latitude and longitude coordinates "
"to have 1 or 2 dimensions, got {} and "
"{}.".format(lat_coord.ndim, lon_coord.ndim)
)
# # Mask out points outside of extent in source_cs - disabled until
# # a way to specify global/limited extent is agreed upon and code
# # is generalised to handle -180 to +180, 0 to 360 and >360 longitudes.
# source_desired_xy = source_cs.transform_points(target_proj,
# target_x.flatten(),
# target_y.flatten())
# if np.any(source_x < 0.0) and np.any(source_x > 180.0):
# raise ValueError('Unable to handle range of longitude.')
# # This does not work in all cases e.g. lon > 360
# if np.any(source_x > 180.0):
# source_desired_x = (source_desired_xy[:, 0].reshape(ny, nx) +
# 360.0) % 360.0
# else:
# source_desired_x = source_desired_xy[:, 0].reshape(ny, nx)
# source_desired_y = source_desired_xy[:, 1].reshape(ny, nx)
# outof_extent_points = ((source_desired_x < source_x.min()) |
# (source_desired_x > source_x.max()) |
# (source_desired_y < source_y.min()) |
# (source_desired_y > source_y.max()))
# # Make array a mask by default (rather than a single bool) to allow mask
# # to be assigned to slices.
# new_data.mask = np.zeros(new_shape)
# Step through cube data, regrid onto desired projection and insert results
# in new_data array
for index, ll_slice in zip(index_it, slice_it):
# Regrid source data onto target grid
index = list(index)
index[xdim] = slice(None, None)
index[ydim] = slice(None, None)
index = tuple(index) # Numpy>=1.16 : index with tuple, *not* list.
new_data[index] = cartopy.img_transform.regrid(
ll_slice.data,
source_x,
source_y,
source_cs,
target_proj,
target_x,
target_y,
)
# # Mask out points beyond extent
# new_data[index].mask[outof_extent_points] = True
# Remove mask if it is unnecessary
if not np.any(new_data.mask):
new_data = new_data.data
# Create new cube
new_cube = iris.cube.Cube(new_data)
# Add new grid coords
x_coord = iris.coords.DimCoord(
target_x[0, :],
"projection_x_coordinate",
units="m",
coord_system=copy.copy(target_cs),
)
y_coord = iris.coords.DimCoord(
target_y[:, 0],
"projection_y_coordinate",
units="m",
coord_system=copy.copy(target_cs),
)
new_cube.add_dim_coord(x_coord, xdim)
new_cube.add_dim_coord(y_coord, ydim)
# Add resampled lat/lon in original coord system
source_desired_xy = source_cs.transform_points(
target_proj, target_x.flatten(), target_y.flatten()
)
new_lon_points = source_desired_xy[:, 0].reshape(ny, nx)
new_lat_points = source_desired_xy[:, 1].reshape(ny, nx)
new_lon_coord = iris.coords.AuxCoord(
new_lon_points,
standard_name="longitude",
units="degrees",
coord_system=orig_cs,
)
new_lat_coord = iris.coords.AuxCoord(
new_lat_points,
standard_name="latitude",
units="degrees",
coord_system=orig_cs,
)
new_cube.add_aux_coord(new_lon_coord, [ydim, xdim])
new_cube.add_aux_coord(new_lat_coord, [ydim, xdim])
coords_to_ignore = set()
coords_to_ignore.update(cube.coords(contains_dimension=xdim))
coords_to_ignore.update(cube.coords(contains_dimension=ydim))
for coord in cube.dim_coords:
if coord not in coords_to_ignore:
new_cube.add_dim_coord(coord.copy(), cube.coord_dims(coord))
for coord in cube.aux_coords:
if coord not in coords_to_ignore:
new_cube.add_aux_coord(coord.copy(), cube.coord_dims(coord))
discarded_coords = coords_to_ignore.difference([lat_coord, lon_coord])
if discarded_coords:
warnings.warn(
"Discarding coordinates that share dimensions with {} and {}: {}".format(
lat_coord.name(),
lon_coord.name(),
[coord.name() for coord in discarded_coords],
),
category=iris.warnings.IrisIgnoringWarning,
)
# TODO handle derived coords/aux_factories
# Copy metadata across
new_cube.metadata = cube.metadata
return new_cube, extent
def _transform_xy(crs_from, x, y, crs_to):
"""Shorthand function to transform 2d points between coordinate reference systems.
Parameters
----------
crs_from : :class:`cartopy.crs.Projection`
The coordinate reference systems.
x, y : array
Point locations defined in 'crs_from'.
crs_to : :class:`cartopy.crs.Projection`
The coordinate reference systems.
Returns
-------
x, y
Arrays of locations defined in 'crs_to'.
"""
pts = crs_to.transform_points(crs_from, x, y)
return pts[..., 0], pts[..., 1]
def _inter_crs_differentials(crs1, x, y, crs2):
"""Calculate coordinate partial differentials from crs1 to crs2.
Returns dx2/dx1, dy2/dx1, dx2/dy1 and dy2/dy1, at given locations.
Parameters
----------
crs1 : :class:`cartopy.crs.Projection`
The coordinate systems for "from".
x, y : array
Point locations defined in 'crs1'.
crs2 : :class:`cartopy.crs.Projection`
The coordinate systems for "to".
Returns
-------
arrays
(dx2/dx1, dy2/dx1, dx2/dy1, dy2/dy1) at given locations. Each
element of this tuple will be the same shape as the 'x' and 'y'
arrays and will be the partial differentials between the two systems.
"""
# Get locations in target crs.
crs2_x, crs2_y = _transform_xy(crs1, x, y, crs2)
# Define small x-deltas in the source crs.
VECTOR_DELTAS_FACTOR = 360000.0 # Empirical factor to obtain small delta.
delta_x = (crs1.x_limits[1] - crs1.x_limits[0]) / VECTOR_DELTAS_FACTOR
delta_x = delta_x * np.ones(x.shape)
eps = 1e-9
# Reverse deltas where we would otherwise step outside the valid range.
invalid_dx = x + delta_x > crs1.x_limits[1] - eps
delta_x[invalid_dx] = -delta_x[invalid_dx]
# Calculate the transformed point with x = x + dx.
crs2_x2, crs2_y2 = _transform_xy(crs1, x + delta_x, y, crs2)
# Form differentials wrt dx.
dx2_dx = (crs2_x2 - crs2_x) / delta_x
dy2_dx = (crs2_y2 - crs2_y) / delta_x
# Define small y-deltas in the source crs.
delta_y = (crs1.y_limits[1] - crs1.y_limits[0]) / VECTOR_DELTAS_FACTOR
delta_y = delta_y * np.ones(y.shape)
# Reverse deltas where we would otherwise step outside the valid range.
invalid_dy = y + delta_y > crs1.y_limits[1] - eps
delta_y[invalid_dy] = -delta_y[invalid_dy]
# Calculate the transformed point with y = y + dy.
crs2_x2, crs2_y2 = _transform_xy(crs1, x, y + delta_y, crs2)
# Form differentials wrt dy.
dx2_dy = (crs2_x2 - crs2_x) / delta_y
dy2_dy = (crs2_y2 - crs2_y) / delta_y
return dx2_dx, dy2_dx, dx2_dy, dy2_dy
def _crs_distance_differentials(crs, x, y):
"""Calculate d(distance) / d(x) and ... / d(y).
Calculate d(distance) / d(x) and ... / d(y) for a coordinate
reference system at specified locations.
Parameters
----------
crs : :class:`cartopy.crs.Projection`
The coordinate reference system.
x, y : array
Locations at which to calculate the differentials,
defined in 'crs' coordinate reference system.
Returns
-------
(abs(ds/dx), abs(ds/dy))
Numerically approximated partial differentials,
i.e. scaling factors between changes in distance and changes in
coordinate values.
"""
# Make a true-latlon coordinate system for distance calculations.
crs_latlon = ccrs.Geodetic(globe=crs.globe)
# Transform points to true-latlon (just to get the true latitudes).
_, true_lat = _transform_xy(crs, x, y, crs_latlon)
# Get coordinate differentials w.r.t. true-latlon.
dlon_dx, dlat_dx, dlon_dy, dlat_dy = _inter_crs_differentials(crs, x, y, crs_latlon)
# Calculate effective scalings of X and Y coordinates.
lat_factor = np.cos(np.deg2rad(true_lat)) ** 2
ds_dx = np.sqrt(dlat_dx * dlat_dx + dlon_dx * dlon_dx * lat_factor)
ds_dy = np.sqrt(dlat_dy * dlat_dy + dlon_dy * dlon_dy * lat_factor)
return ds_dx, ds_dy
def _transform_distance_vectors(u_dist, v_dist, ds, dx2, dy2):
"""Transform distance vectors to another coordinate reference system.
Transform distance vectors from one coordinate reference system to
another, preserving magnitude and physical direction.
Parameters
----------
u_dist, v_dist : array
Components of each vector along the x and y directions of the source
crs at each location.
ds : `DistanceDifferential`
Distance differentials for the source and the target crs at specified
locations.
dx2, dy2 : `PartialDifferential`
Partial differentials from the source to the target crs.
Returns
-------
tuple
(ut_dist, vt_dist): Tuple of arrays containing the vector components
along the x and y directions of the target crs at each location.
"""
# Scale input distance vectors --> source-coordinate differentials.
u1, v1 = u_dist / ds.dx1, v_dist / ds.dy1
# Transform vectors into the target system.
u2 = dx2.dx1 * u1 + dx2.dy1 * v1
v2 = dy2.dx1 * u1 + dy2.dy1 * v1
# Rescale output coordinate vectors --> target distance vectors.
u2_dist, v2_dist = u2 * ds.dx2, v2 * ds.dy2
return u2_dist, v2_dist
def _transform_distance_vectors_tolerance_mask(src_crs, x, y, tgt_crs, ds, dx2, dy2):
"""Return a mask that can be applied to data array to mask elements.
Return a mask that can be applied to data array to mask elements
where the magnitude of vectors are not preserved due to numerical
errors introduced by the transformation between coordinate systems.
Parameters
----------
src_crs : `cartopy.crs.Projection`
The source coordinate reference systems.
x, y : array
Locations of each vector defined in 'src_crs'.
tgt_crs : `cartopy.crs.Projection`
The target coordinate reference systems.
ds : `DistanceDifferential`
Distance differentials for src_crs and tgt_crs at specified locations.
dx2, dy2 : `PartialDifferential`
Partial differentials from src_crs to tgt_crs.
Returns
-------
2d boolean array that is the same shape as x and y.
"""
if x.shape != y.shape:
raise ValueError(
"Arrays do not have matching shapes. "
"x.shape is {}, y.shape is {}.".format(x.shape, y.shape)
)
ones = np.ones(x.shape)
zeros = np.zeros(x.shape)
u_one_t, v_zero_t = _transform_distance_vectors(ones, zeros, ds, dx2, dy2)
u_zero_t, v_one_t = _transform_distance_vectors(zeros, ones, ds, dx2, dy2)
# Squared magnitudes should be equal to one within acceptable tolerance.
# A value of atol=2e-3 is used, which masks any magnitude changes >0.5%
# (approx percentage - based on experimenting).
sqmag_1_0 = u_one_t**2 + v_zero_t**2
sqmag_0_1 = u_zero_t**2 + v_one_t**2
mask = np.logical_not(
np.logical_and(
np.isclose(sqmag_1_0, ones, atol=2e-3),
np.isclose(sqmag_0_1, ones, atol=2e-3),
)
)
return mask
[docs]
def rotate_winds(u_cube, v_cube, target_cs):
r"""Transform wind vectors to a different coordinate system.
The input cubes contain U and V components parallel to the local X and Y
directions of the input grid at each point.
The output cubes contain the same winds, at the same locations, but
relative to the grid directions of a different coordinate system.
Thus in vector terms, the magnitudes will always be the same, but the
angles can be different.
The outputs retain the original horizontal dimension coordinates, but
also have two 2-dimensional auxiliary coordinates containing the X and
Y locations in the target coordinate system.
Parameters
----------
u_cube :
An instance of :class:`iris.cube.Cube` that contains the x-component
of the vector.
v_cube :
An instance of :class:`iris.cube.Cube` that contains the y-component
of the vector.
target_cs :
An instance of :class:`iris.coord_systems.CoordSystem` that specifies
the new grid directions.
Returns
-------
(u', v') tuple of :class:`iris.cube.Cube`
A (u', v') tuple of :class:`iris.cube.Cube` instances that are the u
and v components in the requested target coordinate system.
The units are the same as the inputs.
Notes
-----
.. note::
The U and V values relate to distance, with units such as 'm s-1'.
These are not the same as coordinate vectors, which transform in a
different manner.
.. note::
The names of the output cubes are those of the inputs, prefixed with
'transformed\_' (e.g. 'transformed_x_wind').
.. note::
This function does not maintain laziness when called; it realises data.
See more at :doc:`/userguide/real_and_lazy_data`.
.. warning::
Conversion between rotated-pole and non-rotated systems can be
expressed analytically. However, this function always uses a numerical
approach. In locations where this numerical approach does not preserve
magnitude to an accuracy of 0.1%, the corresponding elements of the
returned cubes will be masked.
"""
# Check u_cube and v_cube have the same shape. We iterate through
# the u and v cube slices which relies on the shapes matching.
if u_cube.shape != v_cube.shape:
msg = (
"Expected u and v cubes to have the same shape. "
"u cube has shape {}, v cube has shape {}."
)
raise ValueError(msg.format(u_cube.shape, v_cube.shape))
# Check the u_cube and v_cube have the same x and y coords.
msg = (
"Coordinates differ between u and v cubes. Coordinate {!r} from "
"u cube does not equal coordinate {!r} from v cube."
)
if u_cube.coord(axis="x") != v_cube.coord(axis="x"):
raise ValueError(
msg.format(u_cube.coord(axis="x").name(), v_cube.coord(axis="x").name())
)
if u_cube.coord(axis="y") != v_cube.coord(axis="y"):
raise ValueError(
msg.format(u_cube.coord(axis="y").name(), v_cube.coord(axis="y").name())
)
# Check x and y coords have the same coordinate system.
x_coord = u_cube.coord(axis="x")
y_coord = u_cube.coord(axis="y")
if x_coord.coord_system != y_coord.coord_system:
msg = (
"Coordinate systems of x and y coordinates differ. "
"Coordinate {!r} has a coord system of {!r}, but coordinate "
"{!r} has a coord system of {!r}."
)
raise ValueError(
msg.format(
x_coord.name(),
x_coord.coord_system,
y_coord.name(),
y_coord.coord_system,
)
)
# Convert from iris coord systems to cartopy CRSs to access
# transform functionality. Use projection as cartopy
# transform_vectors relies on x_limits and y_limits.
if x_coord.coord_system is not None:
src_crs = x_coord.coord_system.as_cartopy_projection()
else:
# Default to Geodetic (but actually use PlateCarree as a
# projection is needed).
src_crs = ccrs.PlateCarree()
target_crs = target_cs.as_cartopy_projection()
# Check the number of dimensions of the x and y coords is the same.
# Subsequent logic assumes either both 1d or both 2d.
x = x_coord.points
y = y_coord.points
if x.ndim != y.ndim or x.ndim > 2 or y.ndim > 2:
msg = (
"x and y coordinates must have the same number of dimensions "
"and be either 1D or 2D. The number of dimensions are {} and "
"{}, respectively.".format(x.ndim, y.ndim)
)
raise ValueError(msg)
# Check the dimension mappings match between u_cube and v_cube.
if u_cube.coord_dims(x_coord) != v_cube.coord_dims(x_coord):
raise ValueError(
"Dimension mapping of x coordinate differs between u and v cubes."
)
if u_cube.coord_dims(y_coord) != v_cube.coord_dims(y_coord):
raise ValueError(
"Dimension mapping of y coordinate differs between u and v cubes."
)
x_dims = u_cube.coord_dims(x_coord)
y_dims = u_cube.coord_dims(y_coord)
# Convert points to 2D, if not already, and determine dims.
if x.ndim == y.ndim == 1:
x, y = _meshgrid(x, y)
dims = (y_dims[0], x_dims[0])
else:
dims = x_dims
# Transpose x, y 2d arrays to match the order in cube's data
# array so that x, y and the sliced data all line up.
if dims[0] > dims[1]:
x = x.transpose()
y = y.transpose()
# Create resulting cubes - produce lazy output data if at least
# one input cube has lazy data
lazy_output = u_cube.has_lazy_data() or v_cube.has_lazy_data()
if lazy_output:
ut_cube = u_cube.copy(data=da.empty_like(u_cube.lazy_data()))
vt_cube = v_cube.copy(data=da.empty_like(v_cube.lazy_data()))
else:
ut_cube = u_cube.copy()
vt_cube = v_cube.copy()
ut_cube.rename("transformed_{}".format(u_cube.name()))
vt_cube.rename("transformed_{}".format(v_cube.name()))
# Get distance scalings for source crs.
ds_dx1, ds_dy1 = _crs_distance_differentials(src_crs, x, y)
# Get distance scalings for target crs.
x2, y2 = _transform_xy(src_crs, x, y, target_crs)
ds_dx2, ds_dy2 = _crs_distance_differentials(target_crs, x2, y2)
ds = DistanceDifferential(ds_dx1, ds_dy1, ds_dx2, ds_dy2)
# Calculate coordinate partial differentials from source crs to target crs.
dx2_dx1, dy2_dx1, dx2_dy1, dy2_dy1 = _inter_crs_differentials(
src_crs, x, y, target_crs
)
dx2 = PartialDifferential(dx2_dx1, dx2_dy1)
dy2 = PartialDifferential(dy2_dx1, dy2_dy1)
# Calculate mask based on preservation of magnitude.
mask = _transform_distance_vectors_tolerance_mask(
src_crs, x, y, target_crs, ds, dx2, dy2
)
apply_mask = mask.any()
if apply_mask:
# Make masked arrays to accept masking.
if lazy_output:
ut_cube = ut_cube.copy(data=da.ma.empty_like(ut_cube.core_data()))
vt_cube = vt_cube.copy(data=da.ma.empty_like(vt_cube.core_data()))
else:
ut_cube.data = ma.asanyarray(ut_cube.data)
vt_cube.data = ma.asanyarray(vt_cube.data)
# Project vectors with u, v components one horiz slice at a time and
# insert into the resulting cubes.
shape = list(u_cube.shape)
for dim in dims:
shape[dim] = 1
ndindex = np.ndindex(*shape)
for index in ndindex:
index = list(index)
for dim in dims:
index[dim] = slice(None, None)
index = tuple(index)
u = u_cube.core_data()[index]
v = v_cube.core_data()[index]
ut, vt = _transform_distance_vectors(u, v, ds, dx2, dy2)
if apply_mask:
if lazy_output:
ut = da.ma.masked_array(ut, mask=mask)
vt = da.ma.masked_array(vt, mask=mask)
else:
ut = ma.asanyarray(ut)
ut[mask] = ma.masked
vt = ma.asanyarray(vt)
vt[mask] = ma.masked
ut_cube.core_data()[index] = ut
vt_cube.core_data()[index] = vt
# Calculate new coords of locations in target coordinate system.
xyz_tran = target_crs.transform_points(src_crs, x, y)
xt = xyz_tran[..., 0].reshape(x.shape)
yt = xyz_tran[..., 1].reshape(y.shape)
# Transpose xt, yt 2d arrays to match the dim order
# of the original x an y arrays - i.e. undo the earlier
# transpose (if applied).
if dims[0] > dims[1]:
xt = xt.transpose()
yt = yt.transpose()
xt_coord = iris.coords.AuxCoord(
xt, standard_name="projection_x_coordinate", coord_system=target_cs
)
yt_coord = iris.coords.AuxCoord(
yt, standard_name="projection_y_coordinate", coord_system=target_cs
)
# Set units based on coord_system.
if isinstance(
target_cs,
(iris.coord_systems.GeogCS, iris.coord_systems.RotatedGeogCS),
):
xt_coord.units = yt_coord.units = "degrees"
else:
xt_coord.units = yt_coord.units = "m"
ut_cube.add_aux_coord(xt_coord, dims)
ut_cube.add_aux_coord(yt_coord, dims)
vt_cube.add_aux_coord(xt_coord.copy(), dims)
vt_cube.add_aux_coord(yt_coord.copy(), dims)
return ut_cube, vt_cube
-3.10.0.dev18-gold?style=flat){kind=icon}