import logging

from dataclasses import dataclass

from functools import reduce

from operator import mul as operator_mul

import cfdm

import numpy as np

from cfdm import is_log_level_debug, is_log_level_detail, is_log_level_info

from . import (

AuxiliaryCoordinate,

Bounds,

CellMeasure,

CellMethod,

Constructs,

Count,

DimensionCoordinate,

Domain,

DomainAncillary,

DomainAxis,

FieldList,

Flags,

Index,

List,

Quantization,

mixin,

)

from .constants import masked as cf_masked

from .data import Data

from .data.array import (

GatheredArray,

RaggedContiguousArray,

RaggedIndexedArray,

RaggedIndexedContiguousArray,

)

from .decorators import (

_deprecated_kwarg_check,

_inplace_enabled,

_inplace_enabled_define_and_cleanup,

_manage_log_level_via_verbosity,

)

from .formula_terms import FormulaTerms

from .functions import (

_DEPRECATION_ERROR,

_DEPRECATION_ERROR_ARG,

_DEPRECATION_ERROR_ATTRIBUTE,

_DEPRECATION_ERROR_KWARG_VALUE,

_DEPRECATION_ERROR_KWARGS,

_DEPRECATION_ERROR_METHOD,

DeprecationError,

_section,

flat,

parse_indices,

)

from .functions import relaxed_identities as cf_relaxed_identities

from .functions import size as cf_size

from .query import Query, eq, ge, gt, le, lt

from .subspacefield import SubspaceField

from .timeduration import TimeDuration

from .units import Units

logger = logging.getLogger(__name__)

# --------------------------------------------------------------------

# Commonly used units

# --------------------------------------------------------------------

# _units_degrees = Units("degrees")

_units_radians = Units("radians")

_units_metres = Units("m")

_units_1 = Units("1")

# --------------------------------------------------------------------

# Map each allowed input collapse method name to its corresponding

# Data method. Input collapse methods not in this sictionary are

# assumed to have a corresponding Data method with the same name.

# --------------------------------------------------------------------

_collapse_methods = {

**{

name: name

for name in [

"mean", # results in 'mean': 'mean' entry, etc.

"mean_absolute_value",

"mean_of_upper_decile",

"max",

"maximum_absolute_value",

"min",

"max",

"minimum_absolute_value",

"mid_range",

"range",

"median",

"sd",

"sum",

"sum_of_squares",

"integral",

"root_mean_square",

"var",

"sample_size",

"sum_of_weights",

"sum_of_weights2",

]

},

**{ # non-identical mapped names:

"avg": "mean",

"average": "mean",

"maximum": "max",

"minimum": "min",

"standard_deviation": "sd",

"variance": "var",

},

}

# --------------------------------------------------------------------

# Map each allowed input collapse method name to its corresponding CF

# cell method.

# --------------------------------------------------------------------

_collapse_cell_methods = {

**{

name: name

for name in [

"point",

"mean",

"mean_absolute_value",

"mean_of_upper_decile",

"maximum",

"maximum_absolute_value",

"minimum",

"minimum_absolute_value",

"mid_range",

"range",

"median",

"standard_deviation",

"sum",

"root_mean_square",

"sum_of_squares",

"variance",

]

},

**{ # non-identical mapped names:

"avg": "mean",

"average": "mean",

"max": "maximum",

"min": "minimum",

"sd": "standard_deviation",

"integral": "sum",

"var": "variance",

"sample_size": "point",

"sum_of_weights": "sum",

"sum_of_weights2": "sum",

},

}

# --------------------------------------------------------------------

# These Data methods may be weighted

# --------------------------------------------------------------------

_collapse_weighted_methods = set(

(

"mean",

"mean_absolute_value",

"mean_of_upper_decile",

"avg",

"average",

"sd",

"standard_deviation",

"sum",

"var",

"variance",

"sum_of_weights",

"sum_of_weights2",

"integral",

"root_mean_square",

)

)

# --------------------------------------------------------------------

# These Data methods may specify a number of degrees of freedom

# --------------------------------------------------------------------

_collapse_ddof_methods = set(("sd", "var"))

_earth_radius = 6371229.0

_relational_methods = (

"__eq__",

"__ne__",

"__lt__",

"__le__",

"__gt__",

"__ge__",

)

@dataclass()

class _Axis_characterisation:

"""Characterise a domain axis.

Used by `_binary_operation` to help with ascertaining if there is

a common axis in two fields.

.. versionaddedd:: 3.16.3

"""

# The size of the axis, an integer.

size: int = -1

# The domain axis identifier. E.g. 'domainaxis0'

axis: str = ""

# The coordinate constructs that characterize the axis

coords: tuple = ()

# The identifiers of the coordinate

# constructs. E.g. ('dimensioncoordinate1',)

keys: tuple = ()

# Whether or not the axis is spanned by the field's data array

axis_in_data_axes: bool = True

class Field(mixin.FieldDomain, mixin.PropertiesData, cfdm.Field):

"""A field construct of the CF data model.

The field construct is central to the CF data model, and includes

all the other constructs. A field corresponds to a CF-netCDF data

variable with all of its metadata. All CF-netCDF elements are

mapped to a field construct or some element of the CF field

construct. The field construct contains all the data and metadata

which can be extracted from the file using the CF conventions.

The field construct consists of a data array and the definition of

its domain (that describes the locations of each cell of the data

array), field ancillary constructs containing metadata defined

over the same domain, and cell method constructs to describe how

the cell values represent the variation of the physical quantity

within the cells of the domain. The domain is defined collectively

by the following constructs of the CF data model: domain axis,

dimension coordinate, auxiliary coordinate, cell measure,

coordinate reference and domain ancillary constructs.

The field construct also has optional properties to describe

aspects of the data that are independent of the domain. These

correspond to some netCDF attributes of variables (e.g. units,

long_name and standard_name), and some netCDF global file

attributes (e.g. history and institution).

**NetCDF interface**

{{netCDF variable}}

{{netCDF global attributes}}

{{netCDF group attributes}}

{{netCDF geometry group}}

Some components exist within multiple constructs, but when written

to a netCDF dataset the netCDF names associated with such

components will be arbitrarily taken from one of them. The netCDF

variable, dimension and sample dimension names and group

structures for such components may be set or removed consistently

across all such components with the `nc_del_component_variable`,

`nc_set_component_variable`, `nc_set_component_variable_groups`,

`nc_clear_component_variable_groups`,

`nc_del_component_dimension`, `nc_set_component_dimension`,

`nc_set_component_dimension_groups`,

`nc_clear_component_dimension_groups`,

`nc_del_component_sample_dimension`,

`nc_set_component_sample_dimension`,

`nc_set_component_sample_dimension_groups`,

`nc_clear_component_sample_dimension_groups` methods.

CF-compliance issues for field constructs read from a netCDF

dataset may be accessed with the `dataset_compliance` method.

"""

def __new__(cls, *args, **kwargs):

"""Store component classes."""

instance = super().__new__(cls)

instance._AuxiliaryCoordinate = AuxiliaryCoordinate

instance._DimensionCoordinate = DimensionCoordinate

instance._Bounds = Bounds

instance._Constructs = Constructs

instance._Domain = Domain

instance._DomainAncillary = DomainAncillary

instance._DomainAxis = DomainAxis

instance._Quantization = Quantization

instance._RaggedContiguousArray = RaggedContiguousArray

instance._RaggedIndexedArray = RaggedIndexedArray

instance._RaggedIndexedContiguousArray = RaggedIndexedContiguousArray

instance._GatheredArray = GatheredArray

instance._Count = Count

instance._Index = Index

instance._List = List

return instance

_special_properties = mixin.PropertiesData._special_properties

_special_properties += ("flag_values", "flag_masks", "flag_meanings")

def __init__(

self, properties=None, source=None, copy=True, _use_data=True

):

"""**Initialisation**

:Parameters:

properties: `dict`, optional

Set descriptive properties. The dictionary keys are

property names, with corresponding values. Ignored if the

*source* parameter is set.

*Parameter example:*

``properties={'standard_name': 'air_temperature'}``

Properties may also be set after initialisation with the

`set_properties` and `set_property` methods.

{{init source: optional}}

{{init copy: `bool`, optional}}

"""

super().__init__(

properties=properties,

source=source,

copy=copy,

_use_data=_use_data,

)

if source:

flags = getattr(source, "Flags", None)

if flags is not None:

self.Flags = flags.copy()

def __getitem__(self, indices):

"""Return a subspace of the field construct defined by indices.

f.__getitem__(indices) <==> f[indices]

Subspacing by indexing uses rules that are very similar to the

numpy indexing rules, the only differences being:

* An integer index i specified for a dimension reduces the size of

this dimension to unity, taking just the i-th element, but keeps

the dimension itself, so that the rank of the array is not

reduced.

* When two or more dimensions’ indices are sequences of integers

then these indices work independently along each dimension

(similar to the way vector subscripts work in Fortran). This is

the same indexing behaviour as on a Variable object of the

netCDF4 package.

* For a dimension that is cyclic, a range of indices specified by

a `slice` that spans the edges of the data (such as ``-2:3`` or

``3:-2:-1``) is assumed to "wrap" around, rather then producing

a null result.

.. seealso:: `indices`, `squeeze`, `subspace`, `where`

**Examples**

>>> f.shape

(12, 73, 96)

>>> f[0].shape

(1, 73, 96)

>>> f[3, slice(10, 0, -2), 95:93:-1].shape

(1, 5, 2)

>>> f.shape

(12, 73, 96)

>>> f[:, [0, 72], [5, 4, 3]].shape

(12, 2, 3)

>>> f.shape

(12, 73, 96)

>>> f[...].shape

(12, 73, 96)

>>> f[slice(0, 12), :, 10:0:-2].shape

(12, 73, 5)

>>> f[[True, True, False, True, True, False, False, True, True, True,

... True, True]].shape

(9, 64, 128)

>>> f[..., :6, 9:1:-2, [1, 3, 4]].shape

(6, 4, 3)

"""

if indices is Ellipsis:

return self.copy()

data = self.data

shape = data.shape

# Parse the index

if not isinstance(indices, tuple):

indices = (indices,)

if isinstance(indices[0], str) and indices[0] == "mask":

ancillary_mask = indices[:2]

indices2 = indices[2:]

else:

ancillary_mask = None

indices2 = indices

indices, roll = parse_indices(shape, indices2, cyclic=True)

if roll:

new = self

axes = data._axes

cyclic_axes = data._cyclic

for iaxis, shift in roll.items():

axis = axes[iaxis]

if axis not in cyclic_axes:

_ = self.get_data_axes()[iaxis]

raise IndexError(

"Can't take a cyclic slice from non-cyclic "

f"{self.constructs.domain_axis_identity(_)!r} axis"

)

new = new.roll(axis=iaxis, shift=shift)

else:

new = self.copy()

data = new.data

# ------------------------------------------------------------

# Subspace the field construct's data

# ------------------------------------------------------------

if ancillary_mask:

ancillary_mask = list(ancillary_mask)

findices = ancillary_mask + indices

else:

findices = indices

new_data = data[tuple(findices)]

if 0 in new_data.shape:

raise IndexError(

f"Indices {findices!r} result in a subspaced shape of "

f"{new_data.shape}, but can't create a subspace of "

f"{self.__class__.__name__} that has a size 0 axis"

)

# Set sizes of domain axes

data_axes = new.get_data_axes()

domain_axes = new.domain_axes(todict=True)

for axis, size in zip(data_axes, new_data.shape):

domain_axes[axis].set_size(size)

# Record which axes were cyclic before the subspace

org_cyclic = [data_axes.index(axis) for axis in new.cyclic()]

# Set the subspaced data

new.set_data(new_data, axes=data_axes, copy=False)

# Update axis cylcicity. Note that this can only entail

# setting an originally cyclic axis to be non-cyclic. Doing

# this now enables us to disable the (possibly very slow)

# automatic check for cyclicity on the 'set_construct' calls

# below.

if org_cyclic:

new_cyclic = new_data.cyclic()

[

new.cyclic(i, iscyclic=False)

for i in org_cyclic

if i not in new_cyclic

]

# ------------------------------------------------------------

# Subspace constructs with data

# ------------------------------------------------------------

if data_axes:

construct_data_axes = new.constructs.data_axes()

for key, construct in new.constructs.filter_by_axis(

*data_axes, axis_mode="or", todict=True

).items():

construct_axes = construct_data_axes[key]

dice = []

needs_slicing = False

for axis in construct_axes:

if axis in data_axes:

needs_slicing = True

dice.append(indices[data_axes.index(axis)])

else:

dice.append(slice(None))

# Generally we do not apply an ancillary mask to the

# metadata items, but for DSGs we do.

if ancillary_mask and new.DSG:

item_mask = []

for mask in ancillary_mask[1]:

iaxes = [

data_axes.index(axis)

for axis in construct_axes

if axis in data_axes

]

for i, (axis, size) in enumerate(

zip(data_axes, mask.shape)

):

if axis not in construct_axes:

if size > 1:

iaxes = None

break

mask = mask.squeeze(i)

if iaxes is None:

item_mask = None

break

else:

mask1 = mask.transpose(iaxes)

for i, axis in enumerate(construct_axes):

if axis not in data_axes:

mask1.inset_dimension(i)

item_mask.append(mask1)

if item_mask:

needs_slicing = True

dice = [ancillary_mask[0], item_mask] + dice

# Replace existing construct with its subspace

if needs_slicing:

new.set_construct(

construct[tuple(dice)],

key=key,

axes=construct_axes,

copy=False,

autocyclic={"no-op": True},

)

return new

def __setitem__(self, indices, value):

"""Called to implement assignment to x[indices]=value.

x.__setitem__(indices, value) <==> x[indices]=value

.. versionadded:: 2.0

"""

if isinstance(value, self.__class__):

value = self._conform_for_assignment(value)

try:

data = value.get_data(None, _fill_value=False)

except AttributeError:

pass

else:

if data is None:

raise ValueError(

f"Can't assign to a {self.__class__.__name__} from a "

f"{value.__class__.__name__} with no data"

)

value = data

data = self.get_data(_fill_value=False)

data[indices] = value

# @property

# def _cyclic(self):

# """Storage for axis cyclicity.

#

# Do not change the value in-place.

#

# """

# return self._custom.get("_cyclic", _empty_set)

#

# @_cyclic.setter

# def _cyclic(self, value):

# """value must be a set.

#

# Do not change the value in-place.

#

# """

# self._custom["_cyclic"] = value

#

# @_cyclic.deleter

# def _cyclic(self):

# self._custom["_cyclic"] = _empty_set

def analyse_items(self, relaxed_identities=None):

"""Analyse a domain.

:Returns:

`dict`

A description of the domain.

**Examples**

>>> print(f)

Axes : time(3) = [1979-05-01 12:00:00, ..., 1979-05-03 12:00:00] gregorian

: air_pressure(5) = [850.000061035, ..., 50.0000038147] hPa

: grid_longitude(106) = [-20.5400109887, ..., 25.6599887609] degrees

: grid_latitude(110) = [23.3200002313, ..., -24.6399995089] degrees

Aux coords : latitude(grid_latitude(110), grid_longitude(106)) = [[67.1246607722, ..., 22.8886948065]] degrees_N

: longitude(grid_latitude(110), grid_longitude(106)) = [[-45.98136251, ..., 35.2925499052]] degrees_E

Coord refs : <CF CoordinateReference: rotated_latitude_longitude>

>>> f.analyse_items()

{

'dim_coords': {'dim0': <CF Dim ....>,

'aux_coords': {'N-d': {'aux0': <CF AuxiliaryCoordinate: latitude(110, 106) degrees_N>,

'aux1': <CF AuxiliaryCoordinate: longitude(110, 106) degrees_E>},

'dim0': {'1-d': {},

'N-d': {},},

'dim1': {'1-d': {},

'N-d': {},},

'dim2': {'1-d': {},

'N-d': {'aux0': <CF AuxiliaryCoordinate: latitude(110, 106) degrees_N>,

'aux1': <CF AuxiliaryCoordinate: longitude(110, 106) degrees_E>},},

'dim3': {'1-d': {},

'N-d': {'aux0': <CF AuxiliaryCoordinate: latitude(110, 106) degrees_N>,

'aux1': <CF AuxiliaryCoordinate: longitude(110, 106) degrees_E>},},},

'axis_to_coord': {'dim0': <CF DimensionCoordinate: time(3) gregorian>,

'dim1': <CF DimensionCoordinate: air_pressure(5) hPa>,

'dim2': <CF DimensionCoordinate: grid_latitude(110) degrees>,

'dim3': <CF DimensionCoordinate: grid_longitude(106) degrees>},

'axis_to_id': {'dim0': 'time',

'dim1': 'air_pressure',

'dim2': 'grid_latitude',

'dim3': 'grid_longitude'},

'cell_measures': {'N-d': {},

'dim0': {'1-d': {},

'N-d': {},},

'dim1': {'1-d': {},

'N-d': {},},

'dim2': {'1-d': {},

'N-d': {},},

'dim3': {'1-d': {},

'N-d': {},},

},

'id_to_aux': {},

'id_to_axis': {'air_pressure': 'dim1',

'grid_latitude': 'dim2',

'grid_longitude': 'dim3',

'time': 'dim0'},

'id_to_coord': {'air_pressure': <CF DimensionCoordinate: air_pressure(5) hPa>,

'grid_latitude': <CF DimensionCoordinate: grid_latitude(110) degrees>,

'grid_longitude': <CF DimensionCoordinate: grid_longitude(106) degrees>,

'time': <CF DimensionCoordinate: time(3) gregorian>},

'id_to_key': {'air_pressure': 'dim1',

'grid_latitude': 'dim2',

'grid_longitude': 'dim3',

'time': 'dim0'},

'undefined_axes': [],

'warnings': [],

}

"""

# ------------------------------------------------------------

# Map each axis identity to its identifier, if such a mapping

# exists.

#

# For example:

# >>> id_to_axis

# {'time': 'dim0', 'height': dim1'}

# ------------------------------------------------------------

id_to_axis = {}

# ------------------------------------------------------------

# For each dimension that is identified by a 1-d auxiliary

# coordinate, map its dimension's its identifier.

#

# For example:

# >>> id_to_aux

# {'region': 'aux0'}

# ------------------------------------------------------------

id_to_aux = {}

# ------------------------------------------------------------

# The keys of the coordinate items which provide axis

# identities

#

# For example:

# >>> id_to_key

# {'region': 'aux0'}

# ------------------------------------------------------------

# id_to_key = {}

axis_to_id = {}

# ------------------------------------------------------------

# Map each dimension's identity to the coordinate which

# provides that identity.

#

# For example:

# >>> id_to_coord

# {'time': <CF Coordinate: time(12)>}

# ------------------------------------------------------------

id_to_coord = {}

axis_to_coord = {}

# ------------------------------------------------------------

# List the dimensions which are undefined, in that no unique

# identity can be assigned to them.

#

# For example:

# >>> undefined_axes

# ['dim2']

# ------------------------------------------------------------

undefined_axes = []

# ------------------------------------------------------------

#

# ------------------------------------------------------------

warnings = []

id_to_dim = {}

axis_to_aux = {}

axis_to_dim = {}

if relaxed_identities is None:

relaxed_identities = cf_relaxed_identities()

# dimension_coordinates = self.dimension_coordinates(view=True)

# auxiliary_coordinates = self.auxiliary_coordinates(view=True)

for axis in self.domain_axes(todict=True):

# dims = self.constructs.chain(

# "filter_by_type",

# ("dimension_coordinate",), "filter_by_axis", (axis,)

# mode="and", todict=True

# )

key, dim = self.dimension_coordinate(

item=True, default=(None, None), filter_by_axis=(axis,)

)

if dim is not None:

# This axis of the domain has a dimension coordinate

identity = dim.identity(strict=True, default=None)

if identity is None:

# Dimension coordinate has no identity, but it may

# have a recognised axis.

for ctype in ("T", "X", "Y", "Z"):

if getattr(dim, ctype, False):

identity = ctype

break

if identity is None and relaxed_identities:

identity = dim.identity(relaxed=True, default=None)

if identity:

if identity in id_to_axis:

warnings.append("Field has multiple {identity!r} axes")

axis_to_id[axis] = identity

id_to_axis[identity] = axis

axis_to_coord[axis] = key

id_to_coord[identity] = key

axis_to_dim[axis] = key

id_to_dim[identity] = key

continue

else:

key, aux = self.auxiliary_coordinate(

filter_by_axis=(axis,),

axis_mode="and", # TODO check this "and"

item=True,

default=(None, None),

)

if aux is not None:

# This axis of the domain does not have a

# dimension coordinate but it does have exactly

# one 1-d auxiliary coordinate, so that will do.

identity = aux.identity(strict=True, default=None)

if identity is None and relaxed_identities:

identity = aux.identity(relaxed=True, default=None)

if identity and aux.has_data():

if identity in id_to_axis:

warnings.append(

f"Field has multiple {identity!r} axes"

)

axis_to_id[axis] = identity

id_to_axis[identity] = axis

axis_to_coord[axis] = key

id_to_coord[identity] = key

axis_to_aux[axis] = key

id_to_aux[identity] = key

continue

# Still here? Then this axis is undefined

undefined_axes.append(axis)

return {

"axis_to_id": axis_to_id,

"id_to_axis": id_to_axis,

"axis_to_coord": axis_to_coord,

"axis_to_dim": axis_to_dim,

"axis_to_aux": axis_to_aux,

"id_to_coord": id_to_coord,

"id_to_dim": id_to_dim,

"id_to_aux": id_to_aux,

"undefined_axes": undefined_axes,

"warnings": warnings,

}

def _is_broadcastable(self, shape):

"""Checks the field's data array is broadcastable to a shape.

:Parameters:

shape: sequence of `int`

:Returns:

`bool`

"""

shape0 = getattr(self, "shape", None)

if shape is None:

return False

shape1 = shape

if tuple(shape1) == tuple(shape0):

# Same shape

return True

ndim0 = len(shape0)

ndim1 = len(shape1)

if not ndim0 or not ndim1:

# Either or both is scalar

return True

for setN in set(shape0), set(shape1):

if setN == {1}:

return True

if ndim1 > ndim0:

return False

for n, m in zip(shape1[::-1], shape0[::-1]):

if n != m and n != 1:

return False

return True

def _axis_positions(self, axes, parse=True, return_axes=False):

"""Convert the given axes to their positions in the data.

Any domain axes that are not spanned by the data are ignored.

If there is no data then an empty list is returned.

.. versionadded:: 3.9.0

:Parameters:

axes: (sequence of) `str` or `int`

The axes to be converted. TODO domain axis selection

parse: `bool`, optional

If False then do not parse the *axes*. Parsing should

always occur unless the given *axes* are the output of

a previous call to `parse_axes`. By default *axes* is

parsed by `_parse_axes`.

return_axes: `bool`, optional

If True then also return the domain axis identifiers

corresponding to the positions.

:Returns:

`list` [, `list`]

The domain axis identifiers. If *return_axes* is True

then also return the corresponding domain axis

identifiers.

"""

data_axes = self.get_data_axes(default=None)

if data_axes is None:

return []

if parse:

axes = self._parse_axes(axes)

axes = [a for a in axes if a in data_axes]

positions = [data_axes.index(a) for a in axes]

if return_axes:

return positions, axes

return positions

def _binary_operation(self, other, method):

"""Implement binary arithmetic and comparison operations on the

master data array with metadata-aware broadcasting.

It is intended to be called by the binary arithmetic and

comparison methods, such as `__sub__`, `__imul__`, `__rdiv__`,

`__lt__`, etc.

:Parameters:

other: `Field` or `Query` or any object that broadcasts to the field construct's data

method: `str`

The binary arithmetic or comparison method name (such as

``'__idiv__'`` or ``'__ge__'``).

:Returns:

`Field`

The new field, or the same field if the operation was an

in place augmented arithmetic assignment.

**Examples**

>>> h = f._binary_operation(g, '__add__')

>>> h = f._binary_operation(g, '__ge__')

>>> f._binary_operation(g, '__isub__')

>>> f._binary_operation(g, '__rdiv__')

"""

debug = is_log_level_debug(logger)

if isinstance(other, Query):

# --------------------------------------------------------

# Combine the field with a Query object

# --------------------------------------------------------

return NotImplemented

if not isinstance(other, self.__class__):

# --------------------------------------------------------

# Combine the field with anything other than a Query

# object or another field construct

# --------------------------------------------------------

if cf_size(other) == 1:

# ----------------------------------------------------

# No changes to the field metadata constructs are

# required so can use the metadata-unaware parent

# method

# ----------------------------------------------------

if other is None:

other = np.array(None, dtype=object)

other = Data(other)

if other.ndim > 0:

other.squeeze(inplace=True)

return super()._binary_operation(other, method)

if self._is_broadcastable(np.shape(other)):

return super()._binary_operation(other, method)

raise ValueError(

f"Can't combine {self.__class__.__name__!r} with "

f"{other.__class__.__name__!r} due to incompatible data "

f"shapes: {self.shape}, {np.shape(other)})"

)

# ============================================================

# Still here? Then combine the field with another field

# ============================================================

relaxed_identities = cf_relaxed_identities()

units = self.Units

sn = self.get_property("standard_name", None)

ln = self.get_property("long_name", None)

other_sn = other.get_property("standard_name", None)

other_ln = other.get_property("long_name", None)

field1 = other.copy()

inplace = method[2] == "i"

if not inplace:

field0 = self.copy()

else:

field0 = self

# Analyse the two fields' data array dimensions

out0 = {}

out1 = {}

for i, (f, out) in enumerate(zip((field0, field1), (out0, out1))):

data_axes = f.get_data_axes()

for axis in f.domain_axes(todict=True):

identity = None

if f.is_discrete_axis(axis):

# This is a discrete axis whose identity is

# inferred from all of its auxiliary coordinates

x = {}

for key, aux_coord in f.auxiliary_coordinates(

filter_by_axis=(axis,),

axis_mode="and",

todict=True,

).items():

identity = aux_coord.identity(

strict=True, default=None

)

if identity is None and relaxed_identities:

identity = aux_coord.identity(

relaxed=True, default=None

)