diff --git a/lib/eofs/experimental/__init__.py b/lib/eofs/experimental/__init__.py
new file mode 100644
index 0000000..efa4e5c
--- /dev/null
+++ b/lib/eofs/experimental/__init__.py
@@ -0,0 +1,17 @@
+"""Experimental features."""
+# (c) Copyright 2013 Andrew Dawson. All Rights Reserved.
+#
+# This file is part of eofs.
+#
+# eofs is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# eofs is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+# for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with eofs. If not, see .
diff --git a/lib/eofs/experimental/rotation/__init__.py b/lib/eofs/experimental/rotation/__init__.py
new file mode 100644
index 0000000..cbd7d24
--- /dev/null
+++ b/lib/eofs/experimental/rotation/__init__.py
@@ -0,0 +1,37 @@
+"""Rotation of EOFs."""
+# (c) Copyright 2013 Andrew Dawson. All Rights Reserved.
+#
+# This file is part of eofs.
+#
+# eofs is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# eofs is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+# for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with eofs. If not, see .
+from __future__ import absolute_import
+
+from . import standard
+
+
+# Define the objects imported by imports of the form:
+# from eofs.experimental.rotation import *
+__all__ = ['standard']
+
+try:
+ from . import cdms
+ __all__.append('cdms')
+except ImportError:
+ pass
+
+try:
+ from . import iris
+ __all__.append('iris')
+except ImportError:
+ pass
diff --git a/lib/eofs/experimental/rotation/cdms.py b/lib/eofs/experimental/rotation/cdms.py
new file mode 100644
index 0000000..3efd308
--- /dev/null
+++ b/lib/eofs/experimental/rotation/cdms.py
@@ -0,0 +1,252 @@
+# (c) Copyright 2013 Andrew Dawson. All Rights Reserved.
+#
+# This file is part of eofs.
+#
+# eofs is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# eofs is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+# for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with eofs. If not, see .
+from __future__ import absolute_import
+import numpy as np
+import numpy.ma as ma
+from cdms2 import createVariable
+
+from .standard import Rotator as StandardRotator
+
+
+class Rotator(StandardRotator):
+ """Rotate EOFs from the `cdms2` interface."""
+
+ def __init__(self, *args, **kwargs):
+ """
+ Rotator(solver, neofs, method='varimax', scaled=True)
+
+ Create an EOF rotator.
+
+ **Arguments:**
+
+ *solver*
+ An `~eofs.cdms.Eof` instance that can generate the EOFs to
+ be rotated.
+
+ *neofs*
+ Number of EOFs to use for the rotation.
+
+ **Keyword arguments:**
+
+ *method*
+ The method to use for rotation. Currently the only accepted
+ value is 'varimax'.
+
+ *scaled*
+ If *True* the EOFs are multiplied by the square-root of
+ their eigenvalue before the rotation. If *False* the
+ orthonormal EOFs are used for the rotation. Defaults to
+ *True*.
+
+ **Examples:**
+
+ A varimax rotator based on the first 10 scaled EOFs::
+
+ solver = Eof(data, weights='area')
+ rotator = Rotator(solver, 10)
+
+ A varimax rotator based on the first 5 un-scaled EOFs::
+
+ solver = Eof(data, weights='area')
+ rotator = Rotator(solver, 5, scaled=False)
+
+ """
+ super(Rotator, self).__init__(*args, **kwargs)
+
+ def eofs(self, *args, **kwargs):
+ """
+ eofs(neofs=None, renormalize=True)
+
+ Rotated empirical orthogonal functions (EOFs).
+
+ **Optional arguments:**
+
+ *neofs*
+ Number of EOFs to return. Defaults to all EOFs. If the
+ number of EOFs requested is more than the number that are
+ available, then all available EOFs will be returned.
+
+ *renormalize*
+ If *True* and the rotation was based on scaled EOFs then the
+ scaling will be removed and the returned rotated EOFs will
+ be scaled to unit variance. If *False* the returned rotated
+ EOFs will retain the scaling.
+
+ **Returns:**
+
+ *eofs*
+ A `cdms2` variable containing the ordered rotated EOFs. The
+ EOFs are numbered from 0 to *neofs* - 1.
+
+ **Examples:**
+
+ All rotated EOFs with scaling::
+
+ rotator = Rotator(solver, 10, scaled=True)
+ rotated_eofs = rotator.eofs()
+
+ All rotated EOFs with scaling removed::
+
+ rotator = Rotator(solver, 10, scaled=True)
+ rotated_eofs = rotator.eofs(renormalize=True)
+
+ The leading rotated EOF with scaling::
+
+ rotator = Rotator(solver, 10, scaled=True)
+ rotated_eofs = rotator.eofs(neofs=1)
+
+ """
+ return super(Rotator, self).eofs(*args, **kwargs)
+
+ def varianceFraction(self, *args, **kwargs):
+ """
+ varianceFraction(neigs=None)
+
+ Fractional rotated EOF mode variances.
+
+ The fraction of the total variance explained by each rotated EOF
+ mode, values between 0 and 1 inclusive. Only meaningful if the
+ *scaled* option was set to *True* when creating the `Rotator`.
+
+ **Optional argument:**
+
+ *neigs*
+ Number of eigenvalues to return the fractional variance for.
+ Defaults to all eigenvalues. If the number of eigenvalues
+ requested is more than the number that are available, then
+ fractional variances for all available eigenvalues will be
+ returned.
+
+ **Returns:**
+
+ *variance_fractions*
+ A `cdms2` variable containing the fractional variances for
+ each rotated EOF mode. The EOF modes are numbered from 0 to
+ *neigs* - 1.
+
+ **Examples:**
+
+ The fractional variance represented by each rotated EOF mode::
+
+ rotator = Rotator(solver, 10, scaled=True)
+ variance_fractions = rotator.varianceFraction()
+
+ The fractional variance represented by the first rotated EOF mode::
+
+ rotator = Rotator(solver, 10, scaled=True)
+ variance_fractions = rotator.varianceFraction(neigs=1)
+
+ """
+ return super(Rotator, self).varianceFraction(*args, **kwargs)
+
+ def pcs(self, *args, **kwargs):
+ """
+ pcs(npcs=None, normalized=False)
+
+ Principal component time series (PCs).
+
+ The PC time series associated with the rotated EOFs.
+
+ **Optional arguments:**
+
+ *npcs*
+ Number of PCs to retrieve. Defaults to all the PCs. If the
+ number of PCs requested is more than the number that are
+ available, then all available PCs will be returned.
+
+ *normalized*
+ If *True* the returned PCs are scaled to unit variance. If
+ *False* no scaling is done. Defaults to *False*.
+
+ **Returns:**
+
+ *pcs*
+ A `cdms2` variable containing the PCs. The PCs are numbered
+ from 0 to *npcs* - 1.
+
+ **Examples:**
+
+ All un-scaled PCs::
+
+ pcs = rotator.pcs()
+
+ First 3 PCs scaled to unit variance::
+
+ pcs = rotator.pcs(npcs=3, normalized=True)
+
+ """
+ return super(Rotator, self).pcs(*args, **kwargs)
+
+ def _solverdata(self):
+ """Get the raw data from the EOF solver."""
+ return self._solver._solver._data
+
+ def strip_metadata(self, var):
+ """Strip basic metadata from a `cdms2` variable.
+
+ **Argument:**
+
+ *var*
+ A `cdms2` variable to strip.
+
+ **Returns:**
+
+ *data*
+ An array of the raw data contained in *var*.
+
+ *metadata*
+ A dictionary containing basic metadata from *var*. The
+ dictionary has entries:
+ 'dimensions': list of dimension coordinates
+ 'id': variable's id
+ 'standard_name': if *var* has a standard name
+ 'long_name': if *var* has a long name
+
+ """
+ metadata = {}
+ metadata['dimensions'] = var.getAxisList()
+ metadata['id'] = var.id
+ for name in ('standard_name', 'long_name'):
+ var_name = getattr(var, name, None)
+ if var_name is not None:
+ metadata[name] = var_name
+ return var.asma(), metadata
+
+ def apply_metadata(self, data, metadata):
+ """Construct a cube from raw data and a metadata dictionary.
+
+ **Arguments:**
+
+ *data*
+ Raw data array.
+
+ *metadata*
+ A dictionary of metadata returned from
+ `~Rotator.strip_metadata`.
+
+ **Returns:**
+
+ *var*
+ A `cdms2` variable.
+
+ """
+ dims = metadata['dimensions']
+ var = createVariable(data, axes=dims, id=metadata['id'])
+ for name in ('standard_name', 'long_name'):
+ if name in metadata:
+ setattr(var, name, metadata[name])
+ return var
diff --git a/lib/eofs/experimental/rotation/iris.py b/lib/eofs/experimental/rotation/iris.py
new file mode 100644
index 0000000..2eef840
--- /dev/null
+++ b/lib/eofs/experimental/rotation/iris.py
@@ -0,0 +1,253 @@
+"""Rotation of EOFs for the `iris` interface."""
+# (c) Copyright 2013 Andrew Dawson. All Rights Reserved.
+#
+# This file is part of eofs.
+#
+# eofs is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# eofs is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+# for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with eofs. If not, see .
+from __future__ import absolute_import
+import numpy as np
+import numpy.ma as ma
+from iris.cube import Cube
+
+from .standard import Rotator as StandardRotator
+
+
+class Rotator(StandardRotator):
+ """Rotate EOFs from the `iris` interface."""
+
+ def __init__(self, *args, **kwargs):
+ """
+ Rotator(solver, neofs, method='varimax', scaled=True)
+
+ Create an EOF rotator.
+
+ **Arguments:**
+
+ *solver*
+ An `~eofs.iris.Eof` instance that can generate the EOFs to
+ be rotated.
+
+ *neofs*
+ Number of EOFs to use for the rotation.
+
+ **Keyword arguments:**
+
+ *method*
+ The method to use for rotation. Currently the only accepted
+ value is 'varimax'.
+
+ *scaled*
+ If *True* the EOFs are multiplied by the square-root of
+ their eigenvalue before the rotation. If *False* the
+ orthonormal EOFs are used for the rotation. Defaults to
+ *True*.
+
+ **Examples:**
+
+ A varimax rotator based on the first 10 scaled EOFs::
+
+ solver = Eof(data, weights='area')
+ rotator = Rotator(solver, 10)
+
+ A varimax rotator based on the first 5 un-scaled EOFs::
+
+ solver = Eof(data, weights='area')
+ rotator = Rotator(solver, 5, scaled=False)
+
+ """
+ super(Rotator, self).__init__(*args, **kwargs)
+
+ def eofs(self, *args, **kwargs):
+ """
+ eofs(neofs=None, renormalize=True)
+
+ Rotated empirical orthogonal functions (EOFs).
+
+ **Optional arguments:**
+
+ *neofs*
+ Number of EOFs to return. Defaults to all EOFs. If the
+ number of EOFs requested is more than the number that are
+ available, then all available EOFs will be returned.
+
+ *renormalize*
+ If *True* and the rotation was based on scaled EOFs then the
+ scaling will be removed and the returned rotated EOFs will
+ be scaled to unit variance. If *False* the returned rotated
+ EOFs will retain the scaling.
+
+ **Returns:**
+
+ *eofs*
+ A `~iris.cube.Cube` containing the ordered rotated EOFs. The
+ EOFs are numbered from 0 to *neofs* - 1.
+
+ **Examples:**
+
+ All rotated EOFs with scaling::
+
+ rotator = Rotator(solver, 10, scaled=True)
+ rotated_eofs = rotator.eofs()
+
+ All rotated EOFs with scaling removed::
+
+ rotator = Rotator(solver, 10, scaled=True)
+ rotated_eofs = rotator.eofs(renormalize=True)
+
+ The leading rotated EOF with scaling::
+
+ rotator = Rotator(solver, 10, scaled=True)
+ rotated_eofs = rotator.eofs(neofs=1)
+
+ """
+ return super(Rotator, self).eofs(*args, **kwargs)
+
+ def varianceFraction(self, *args, **kwargs):
+ """
+ varianceFraction(neigs=None)
+
+ Fractional rotated EOF mode variances.
+
+ The fraction of the total variance explained by each rotated EOF
+ mode, values between 0 and 1 inclusive. Only meaningful if the
+ *scaled* option was set to *True* when creating the `Rotator`.
+
+ **Optional argument:**
+
+ *neigs*
+ Number of eigenvalues to return the fractional variance for.
+ Defaults to all eigenvalues. If the number of eigenvalues
+ requested is more than the number that are available, then
+ fractional variances for all available eigenvalues will be
+ returned.
+
+ **Returns:**
+
+ *variance_fractions*
+ A `~iris.cube.Cube` containing the fractional variances for
+ each rotated EOF mode. The EOF modes are numbered from 0 to
+ *neigs* - 1.
+
+ **Examples:**
+
+ The fractional variance represented by each rotated EOF mode::
+
+ rotator = Rotator(solver, 10, scaled=True)
+ variance_fractions = rotator.varianceFraction()
+
+ The fractional variance represented by the first rotated EOF mode::
+
+ rotator = Rotator(solver, 10, scaled=True)
+ variance_fractions = rotator.varianceFraction(neigs=1)
+
+ """
+ return super(Rotator, self).varianceFraction(*args, **kwargs)
+
+ def pcs(self, *args, **kwargs):
+ """
+ pcs(npcs=None, normalized=False)
+
+ Principal component time series (PCs).
+
+ The PC time series associated with the rotated EOFs.
+
+ **Optional arguments:**
+
+ *npcs*
+ Number of PCs to retrieve. Defaults to all the PCs. If the
+ number of PCs requested is more than the number that are
+ available, then all available PCs will be returned.
+
+ *normalized*
+ If *True* the returned PCs are scaled to unit variance. If
+ *False* no scaling is done. Defaults to *False*.
+
+ **Returns:**
+
+ *pcs*
+ A `~iris.cube.Cube` containing the PCs. The PCs are numbered
+ from 0 to *npcs* - 1.
+
+ **Examples:**
+
+ All un-scaled PCs::
+
+ pcs = rotator.pcs()
+
+ First 3 PCs scaled to unit variance::
+
+ pcs = rotator.pcs(npcs=3, normalized=True)
+
+ """
+ return super(Rotator, self).pcs(*args, **kwargs)
+
+ def _solverdata(self):
+ """Get the raw data from the EOF solver."""
+ return self._solver._solver._data
+
+ def strip_metadata(self, cube):
+ """Strip basic metadata from a `~iris.cube.Cube`.
+
+ **Argument:**
+
+ *cube*
+ A `~iris.cube.Cube` to strip.
+
+ **Returns:**
+
+ *data*
+ An array of the raw data contained in *cube*.
+
+ *metadata*
+ A dictionary containing basic metadata from *cube*. The
+ dictionary has entries:
+ 'dimensions': tuple of dimension coordinates
+ 'standard_name': if *cube* has a standard name
+ 'long_name': if *cube* has a long name
+ 'var_name': if cube has a var_name
+
+ """
+ metadata = {}
+ metadata['dimensions'] = cube.dim_coords
+ for name in ('standard_name', 'long_name', 'var_name'):
+ cube_name = getattr(cube, name, None)
+ if cube_name is not None:
+ metadata[name] = cube_name
+ return cube.data, metadata
+
+ def apply_metadata(self, data, metadata):
+ """Construct a cube from raw data and a metadata dictionary.
+
+ **Arguments:**
+
+ *data*
+ Raw data array.
+
+ *metadata*
+ A dictionary of metadata returned from
+ `~Rotator.strip_metadata`.
+
+ **Returns:**
+
+ *cube*
+ A `~iris.cube.Cube`.
+
+ """
+ dims = metadata['dimensions']
+ cube = Cube(data,
+ dim_coords_and_dims=zip(dims, range(len(dims))))
+ for name in ('standard_name', 'long_name', 'var_name'):
+ if name in metadata:
+ setattr(cube, name, metadata[name])
+ return cube
diff --git a/lib/eofs/experimental/rotation/kernels.py b/lib/eofs/experimental/rotation/kernels.py
new file mode 100644
index 0000000..f823f9b
--- /dev/null
+++ b/lib/eofs/experimental/rotation/kernels.py
@@ -0,0 +1,82 @@
+"""Rotation kernels for EOF rotation."""
+# (c) Copyright 2013 Andrew Dawson. All Rights Reserved.
+#
+# This file is part of eofs.
+#
+# eofs is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# eofs is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+# for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with eofs. If not, see .
+from __future__ import absolute_import
+
+import numpy as np
+
+
+def _varimax_kernel(eofs, eps=1e-10, itermax=1000, kaisernorm=True):
+ """Rotation of EOFs according to the varimax criterion.
+
+ **Arguments:**
+
+ *eofs*
+ A 2-dimensional `~numpy.ndarray` with dimensions [neofs, nspace]
+ containing no missing values.
+
+ **Optional arguments:**
+
+ *epsilon*
+ Tolerance value used to determine convergence of the rotation
+ algorithm. Defaults to 1e-10.
+
+ *itermax*
+ Maximum number of iterations to allow in the rotation algorithm.
+
+ *kaisernorm*
+ If *True* uses Kaiser row normalization. If *False* no
+ normalization is used in the kernel. Defaults to *True*.
+
+ """
+ try:
+ neofs, nspace = eofs.shape
+ except ValueError:
+ raise ValueError('kernel requires a 2-D input')
+ if neofs < 2:
+ raise ValueError('at least 2 EOFs are required for rotation')
+ if kaisernorm:
+ # Apply Kaiser row normalization.
+ scale = np.sqrt((eofs ** 2).sum(axis=0))
+ eofs /= scale
+ # Define the initial values of the rotation matrix and the convergence
+ # monitor.
+ rotation = np.eye(neofs, dtype=eofs.dtype)
+ delta = 0.
+ # Iteratively compute the rotation matrix.
+ for i in range(itermax):
+ z = np.dot(eofs.T, rotation)
+ b = np.dot(eofs,
+ z ** 3 - np.dot(z, np.diag((z ** 2).sum(axis=0)) / nspace))
+ u, s, v = np.linalg.svd(b)
+ rotation = np.dot(u, v)
+ delta_previous = delta
+ delta = s.sum()
+ if delta < delta_previous * (1. + eps):
+ # Convergence is reached, stop the iteration.
+ break
+ # Apply the rotation to the input EOFs.
+ reofs = np.dot(eofs.T, rotation).T
+ if kaisernorm:
+ # Remove the normalization.
+ reofs *= scale
+ return reofs, rotation
+
+
+KERNEL_MAPPING = {
+ 'varimax': _varimax_kernel,
+}
diff --git a/lib/eofs/experimental/rotation/standard.py b/lib/eofs/experimental/rotation/standard.py
new file mode 100644
index 0000000..cce12aa
--- /dev/null
+++ b/lib/eofs/experimental/rotation/standard.py
@@ -0,0 +1,635 @@
+"""Rotation of EOFs."""
+# (c) Copyright 2013 Andrew Dawson. All Rights Reserved.
+#
+# This file is part of eofs.
+#
+# eofs is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# eofs is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+# for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with eofs. If not, see .
+from __future__ import absolute_import
+import collections
+import warnings
+
+import numpy as np
+import numpy.ma as ma
+
+from .kernels import KERNEL_MAPPING
+from ...tools.standard import correlation_map, covariance_map
+
+
+class Rotator(object):
+ """Rotate EOFs from the standard `numpy` interface."""
+
+ def __init__(self, solver, neofs, method='varimax', scaled=True,
+ kernelargs=None):
+ """Create an EOF rotator.
+
+ **Arguments:**
+
+ *solver*
+ An `~eofs.standard.Eof` instance that can generate the EOFs
+ to be rotated.
+
+ *neofs*
+ Number of EOFs to use for the rotation.
+
+ **Keyword arguments:**
+
+ *method*
+ The method to use for rotation. Currently the only accepted
+ value is 'varimax'.
+
+ *scaled*
+ If *True* the EOFs are multiplied by the square-root of
+ their eigenvalue before the rotation. If *False* the
+ orthonormal EOFs are used for the rotation. Defaults to
+ *True*.
+
+ **Examples:**
+
+ A varimax rotator based on the first 10 scaled EOFs::
+
+ solver = Eof(data, weights=weights_array)
+ rotator = Rotator(solver, 10)
+
+ A varimax rotator based on the first 5 un-scaled EOFs::
+
+ solver = Eof(data, weights=weights_array)
+ rotator = Rotator(solver, 5, scaled=False)
+
+ """
+ # This is just a reference, not a copy, so it should be lightweight.
+ self._solver = solver
+ self._scaled = scaled
+ self.neofs = neofs
+ # Retrieve the required quantities from the solver.
+ eofscaling = 2 if scaled else 0
+ eofs = self._solver.eofs(eofscaling=eofscaling, neofs=neofs)
+ # Remove metadata and store it for later use.
+ eofs, self._eofmetadata = self.strip_metadata(eofs)
+ # Remove missing values and metadata and reshape to [neofs, nspace].
+ eofs, self._eofinfo = self._to2d(eofs)
+ # Do the rotation using the appropriate kernel.
+ kwargs = {}
+ if kernelargs is not None:
+ try:
+ kwargs.update(kernelargs)
+ except TypeError:
+ raise TypeError('kernel arguments must be a '
+ 'dictionary of keyword arguments')
+ try:
+ self._eofs_rot, R = KERNEL_MAPPING[method.lower()](eofs, **kwargs)
+ self._rotation_matrix = R
+ except KeyError:
+ raise ValueError("unknown rotation method: '{!s}'".format(method))
+ # Compute variances of the rotated EOFs as these are used by several
+ # methods.
+ self._eofs_rot_var = (self._eofs_rot ** 2).sum(axis=1)
+ self._var_idx = np.argsort(self._eofs_rot_var)[::-1]
+ # Reorder rotated EOFs according to their variance.
+ self._eofs_rot = self._eofs_rot[self._var_idx]
+ self._eofs_rot_var = self._eofs_rot_var[self._var_idx]
+
+ def eofs(self, neofs=None, renormalize=False):
+ """Rotated empirical orthogonal functions (EOFs).
+
+ **Optional arguments:**
+
+ *neofs*
+ Number of EOFs to return. Defaults to all EOFs. If the
+ number of EOFs requested is more than the number that are
+ available, then all available EOFs will be returned.
+
+ *renormalize*
+ If *True* and the rotation was based on scaled EOFs then the
+ scaling will be removed and the returned rotated EOFs will
+ be scaled to unit variance. If *False* the returned rotated
+ EOFs will retain the scaling.
+
+ **Returns:**
+
+ *eofs*
+ An array with the ordered rotated EOFs along the first
+ dimension.
+
+ **Examples:**
+
+ All rotated EOFs with scaling::
+
+ rotator = Rotator(solver, 10, scaled=True)
+ rotated_eofs = rotator.eofs()
+
+ All rotated EOFs with scaling removed::
+
+ rotator = Rotator(solver, 10, scaled=True)
+ rotated_eofs = rotator.eofs(renormalize=True)
+
+ The leading rotated EOF with scaling::
+
+ rotator = Rotator(solver, 10, scaled=True)
+ rotated_eofs = rotator.eofs(neofs=1)
+
+ """
+ # Determine the correct slice.
+ if (neofs is None) or neofs > self.neofs:
+ neofs = self.neofs
+ slicer = slice(0, neofs)
+ # Optionally renormalize.
+ if renormalize and self._scaled:
+ eofs_rot = self._eofs_rot / \
+ np.sqrt(self._eofs_rot_var)[:, np.newaxis]
+ else:
+ eofs_rot = self._eofs_rot
+ eofs_rot = self._from2d(eofs_rot, self._eofinfo)[slicer]
+ eofs_rot = self.apply_metadata(eofs_rot, self._eofmetadata)
+ return eofs_rot
+
+ def varianceFraction(self, neigs=None):
+ """Fractional rotated EOF mode variances.
+
+ The fraction of the total variance explained by each rotated EOF
+ mode, values between 0 and 1 inclusive. Only meaningful if the
+ *scaled* option was set to *True* when creating the `Rotator`.
+
+ **Optional argument:**
+
+ *neigs*
+ Number of eigenvalues to return the fractional variance for.
+ Defaults to all eigenvalues. If the number of eigenvalues
+ requested is more than the number that are available, then
+ fractional variances for all available eigenvalues will be
+ returned.
+
+ **Returns:**
+
+ *variance_fractions*
+ An array containing the fractional variances.
+
+ **Examples:**
+
+ The fractional variance represented by each rotated EOF mode::
+
+ rotator = Rotator(solver, 10, scaled=True)
+ variance_fractions = rotator.varianceFraction()
+
+ The fractional variance represented by the first rotated EOF mode::
+
+ rotator = Rotator(solver, 10, scaled=True)
+ variance_fractions = rotator.varianceFraction(neigs=1)
+
+ """
+ if (neigs is None) or (neigs > self.neofs):
+ neigs = self.neofs
+ slicer = slice(0, neigs)
+ # Compute fractions of variance accounted for by each rotated mode.
+ eigenvalues = self._solver.eigenvalues(neigs=neigs)
+ variance_fractions = self._solver.varianceFraction(neigs=neigs)
+ ev, ev_metadata = self.strip_metadata(eigenvalues)
+ vf, vf_metadata = self.strip_metadata(variance_fractions)
+ if self._scaled:
+ ratio = vf[0] / ev[0]
+ vf_rot = self._eofs_rot_var[slicer] * ratio
+ else:
+ vf_rot = np.array([1. / float(self._eofs_rot.shape[1])] * neigs)
+ vf_rot = self.apply_metadata(vf_rot, vf_metadata)
+ return vf_rot
+
+ def eigenvalues(self, neigs=None):
+ """Variances of the rotated EOFs."""
+ if neigs > self.neofs or neigs is None:
+ neigs = self.neofs
+ slicer = slice(0, neigs)
+ eigenvalues = self._solver.eigenvalues(neigs=neigs)
+ ev, ev_metadata = self.strip_metadata(eigenvalues)
+ variances = self._eofs_rot_var[slicer]
+ variances = self.apply_metadata(variances, ev_metadata)
+ return variances
+
+ def pcs(self, npcs=None, normalized=False):
+ """Principal component time series (PCs).
+
+ The PC time series associated with the rotated EOFs.
+
+ **Optional arguments:**
+
+ *npcs*
+ Number of PCs to retrieve. Defaults to all the PCs. If the
+ number of PCs requested is more than the number that are
+ available, then all available PCs will be returned.
+
+ *normalized*
+ If *True* the returned PCs are scaled to unit variance. If
+ *False* no scaling is done. Defaults to *False*.
+
+ **Returns:**
+
+ *pcs*
+ An array where the columns are the ordered PCs.
+
+ **Examples:**
+
+ All un-scaled PCs::
+
+ pcs = rotator.pcs()
+
+ First 3 PCs scaled to unit variance::
+
+ pcs = rotator.pcs(npcs=3, normalized=True)
+
+ """
+ # Determine the correct slice.
+ if (npcs is None) or (npcs > self.neofs):
+ npcs = self.neofs
+ slicer = slice(0, npcs)
+ # Compute the PCs:
+ # 1. Obtain the non-rotated PCs.
+ pcs = self._solver.pcs(npcs=self.neofs, pcscaling=1)
+ # 2. Apply the rotation matrix to standardized PCs.
+ pcs = np.dot(pcs, self._rotation_matrix)
+ # 3. Reorder according to variance.
+ pcs = pcs[:, self._var_idx]
+ if not normalized:
+ # Optionally scale by square root of variance of rotated EOFs.
+ pcs *= np.sqrt(self._eofs_rot_var)
+ # Select only the required PCs.
+ pcs = pcs[:, slicer]
+ # Collect the metadata used for PCs by the solver and apply it to
+ # these PCs.
+ _, pcs_metadata = self.strip_metadata(self._solver.pcs(npcs=npcs))
+ pcs = self.apply_metadata(pcs, pcs_metadata)
+ return pcs
+
+ def eofsAsCorrelation(self, neofs=None):
+ """Correlation map rotated EOFs.
+
+ Rotated empirical orthogonal functions (EOFs) expressed as the
+ correlation between the rotated principal component time series (PCs)
+ and the time series of the `Eof` input *dataset* at each grid
+ point.
+
+ .. note::
+
+ These are not related to the EOFs computed from the
+ correlation matrix.
+
+ **Optional argument:**
+
+ *neofs*
+ Number of EOFs to return. Defaults to all EOFs. If the
+ number of EOFs requested is more than the number that are
+ available, then all available EOFs will be returned.
+
+ **Returns:**
+
+ *eofs*
+ An array with the ordered EOFs along the first dimension.
+
+ **Examples:**
+
+ All EOFs::
+
+ eofs = rotator.eofsAsCorrelation()
+
+ The leading EOF::
+
+ eof1 = rotator.eofsAsCorrelation(neofs=1)
+
+ """
+ # Get original dimensions of data
+ records = self._solver._records
+ originalshape = self._solver._originalshape
+ # Retrieve the specified number of PCs.
+ pcs = self.pcs(npcs=neofs, normalized=True)
+ # Compute the correlation of the PCs with the input field.
+ c = correlation_map(
+ pcs,
+ self._solverdata().reshape((records,) + originalshape))
+ # The results of the correlation_map function will be a masked array.
+ # For consistency with other return values, this is converted to a
+ # numpy array filled with numpy.nan.
+ if not self._solver._filled:
+ c = c.filled(fill_value=np.nan)
+ return c
+
+ def eofsAsCovariance(self, neofs=None, normalized=True):
+ """Covariance map rotated EOFs.
+
+ Rotated empirical orthogonal functions (EOFs) expressed as the
+ covariance between the rotated principal component time series (PCs)
+ and the time series of the `Eof` input *dataset* at each grid
+ point.
+
+ **Optional arguments:**
+
+ *neofs*
+ Number of EOFs to return. Defaults to all EOFs. If the
+ number of EOFs requested is more than the number that are
+ available, then all available EOFs will be returned.
+
+ *pcscaling*
+ Set the scaling of the PCs used to compute covariance. The
+ following values are accepted:
+ *normalized*
+ If *True* the PCs used to compute covariance are scaled to
+ unit variance. If *False* no scaling is done.
+ Defaults to *True* which is the same as scaling option *1*
+ for non-rotated covariance maps.
+
+ **Returns:**
+
+ *eofs*
+ An array with the ordered EOFs along the first dimension.
+
+ **Examples:**
+
+ All EOFs::
+
+ eofs = rotator.eofsAsCovariance()
+
+ The leading EOF::
+
+ eof1 = rotator.eofsAsCovariance(neofs=1)
+
+ The leading EOF using un-scaled PCs::
+
+ eof1 = rotator.eofsAsCovariance(neofs=1, normalized=False)
+
+ """
+ # Get original dimensions of data
+ records = self._solver._records
+ originalshape = self._solver._originalshape
+ # Retrieve the specified number of PCs.
+ pcs = self.pcs(npcs=neofs, normalized=normalized)
+ # Divide the input data by the weighting (if any) before computing
+ # the covariance maps.
+ data = self._solverdata().reshape((records,) + originalshape)
+ if self._solver._weights is not None:
+ with warnings.catch_warnings():
+ warnings.simplefilter('ignore', RuntimeWarning)
+ data /= self._solver._weights
+ c = covariance_map(pcs, data, ddof=self._solver._ddof)
+ # The results of the covariance_map function will be a masked array.
+ # For consitsency with other return values, this is converted to a
+ # numpy array filled with numpy.nan.
+ if not self._solver._filled:
+ c = c.filled(fill_value=np.nan)
+ return c
+
+ def reconstuctedField(self, neofs):
+ """Reconstructed data field based on a subset of rotated EOFs.
+
+ If weights were passed to the `Eof` instance the returned
+ reconstructed field will automatically have this weighting
+ removed. Otherwise the returned field will have the same
+ weighting as the `Eof` input *dataset*.
+
+ **Argument:**
+
+ *neofs*
+ Number of EOFs to use for the reconstruction. If the
+ number of EOFs requested is more than the number that are
+ available, then all available EOFs will be used for the
+ reconstruction. Alternatively this argument can be an
+ iterable of mode numbers (where the first mode is 1) in
+ order to facilitate reconstruction with arbitrary modes.
+
+ **Returns:**
+
+ *reconstruction*
+ An array the same shape as the `Eof` input *dataset*
+ contaning the reconstruction using *neofs* EOFs.
+
+ **Examples:**
+
+ Reconstruct the input field using 3 rotated EOFs::
+
+ reconstruction = rotator.reconstructedField(3)
+
+ Reconstruct the input field using rotated EOFs 1, 2 and 5::
+
+ reconstruction = rotator.reconstuctedField([1, 2, 5])
+
+ """
+ # Determine how the PCs and EOFs will be selected.
+ if isinstance(neofs, collections.Iterable):
+ modes = [m - 1 for m in neofs]
+ else:
+ modes = slice(0, neofs)
+ # Create array containing rotated EOFs including not a number entries
+ # of original input data.
+ originalshape = self._solver._originalshape
+ channels = np.product(originalshape)
+ nan_idx = np.isnan(self._solver._flatE).all(axis=0)
+ L = self._eofs_rot[modes]
+ neofs = L.shape[0]
+ loadings = np.zeros((neofs, channels)) * np.nan
+ loadings[:, ~nan_idx] = L
+ # Project principal components onto the rotated EOFs to
+ # compute the reconstructed field.
+ P = self.pcs(npcs=None, normalized=True)[:, modes]
+ rval = np.dot(P, loadings)
+ # Reshape the reconstructed field so it has the same shape as the
+ # input data set.
+ records = self._solver._records
+ rval = rval.reshape((records,) + originalshape)
+ # Un-weight the reconstructed field if weighting was performed on
+ # the input data set.
+ if self._solver._weights is not None:
+ rval = rval / self._solver._weights
+ # Return the reconstructed field.
+ if self._solver._filled:
+ rval = ma.array(rval, mask=np.where(np.isnan(rval), True, False))
+ return rval
+
+ def projectField(self, field, neofs=None, eofscaling=0, weighted=True):
+ """Project a field onto the rotated EOFs.
+
+ Given a data set, projects it onto the rotated EOFs to generate a
+ corresponding set of pseudo-PCs.
+
+ **Argument:**
+
+ *field*
+ A `numpy.ndarray` or `numpy.ma.MaskedArray` with two or more
+ dimensions containing the data to be projected onto the
+ EOFs. It must have the same corresponding spatial dimensions
+ (including missing values in the same places) as the `Eof`
+ input *dataset*. *field* may have a different length time
+ dimension to the `Eof` input *dataset* or no time dimension
+ at all.
+
+ **Optional arguments:**
+
+ *neofs*
+ Number of EOFs to project onto. Defaults to all EOFs. If the
+ number of EOFs requested is more than the number that are
+ available, then the field will be projected onto all
+ available EOFs.
+
+ *eofscaling*
+ Set the scaling of the EOFs that are projected onto. The
+ following values are accepted:
+
+ * *0* : Un-scaled EOFs (default).
+ * *1* : EOFs are divided by the square-root of their
+ eigenvalue.
+ * *2* : EOFs are multiplied by the square-root of their
+ eigenvalue.
+
+ *weighted*
+ If *True* then *field* is weighted using the same weights
+ used for the EOF analysis prior to projection. If *False*
+ then no weighting is applied. Defaults to *True* (weighting
+ is applied). Generally only the default setting should be
+ used.
+
+ **Returns:**
+
+ *pseudo_pcs*
+ An array where the columns are the ordered pseudo-PCs.
+
+ **Examples:**
+
+ Project a data set onto all EOFs::
+
+ pseudo_pcs = solver.projectField(data)
+
+ Project a data set onto the four leading EOFs::
+
+ pseudo_pcs = solver.projectField(data, neofs=4)
+
+ """
+ # Check that the shape/dimension of the data set is compatible with
+ # the EOFs.
+ solver = self._solver
+ solver._verify_projection_shape(field, solver._originalshape)
+ input_ndim = field.ndim
+ eof_ndim = len(solver._originalshape) + 1
+ # Create a slice object for truncating the EOFs.
+ slicer = slice(0, neofs)
+ # If required, weight the data set with the same weighting that was
+ # used to compute the EOFs.
+ field = field.copy()
+ if weighted:
+ wts = solver.getWeights()
+ if wts is not None:
+ field = field * wts
+ # Fill missing values with NaN if this is a masked array.
+ try:
+ field = field.filled(fill_value=np.nan)
+ except AttributeError:
+ pass
+ # Flatten the data set into [time, space] dimensionality.
+ if eof_ndim > input_ndim:
+ field = field.reshape((1,) + field.shape)
+ records = field.shape[0]
+ channels = np.product(field.shape[1:])
+ field_flat = field.reshape([records, channels])
+ # Locate the non-missing values and isolate them.
+ if not solver._valid_nan(field_flat):
+ raise ValueError('missing values detected in different '
+ 'locations at different times')
+ nonMissingIndex = np.where(np.logical_not(np.isnan(field_flat[0])))[0]
+ field_flat = field_flat[:, nonMissingIndex]
+ # Locate the non-missing values in the EOFs and check they match those
+ # in the data set, then isolate the non-missing values.
+ eofNonMissingIndex = np.where(
+ np.logical_not(np.isnan(solver._flatE[0])))[0]
+ if eofNonMissingIndex.shape != nonMissingIndex.shape or \
+ (eofNonMissingIndex != nonMissingIndex).any():
+ raise ValueError('field and EOFs have different '
+ 'missing value locations')
+ # The correct projection of a new data field on the rotated EOFs
+ # follows the following equation:
+ # PC_new = x_new @ E @ L^(1/2) @ R
+ # where
+ # PC_new : standardized pseudo-PC for new input data field
+ # x_new : new input data field
+ # E : non-rotated EOFs (eigenvectors)
+ # L^(1/2) : Square root of diagonal matrix containing the eigenvalues
+ # R : rotation matrix
+ eofs_flat = solver._flatE[:self.neofs, eofNonMissingIndex]
+ projected_pcs = field_flat @ eofs_flat.T
+ projected_pcs /= np.sqrt(solver._L[:self.neofs])
+ projected_pcs = projected_pcs @ self._rotation_matrix
+ # Reorder the obtained (standardized) rotated EOFs according
+ # to their variance.
+ projected_pcs = projected_pcs[:, self._var_idx]
+ # Select desired PCs
+ projected_pcs = projected_pcs[:, slicer]
+ # PCs are standardized. In order to match the correct eofscaling
+ # we have to multiply the PCs with
+ # the square root of the rotated variance (eofscaling == 1)
+ # the rotated variance (eofscaling == 2)
+ if eofscaling == 0:
+ pass
+ elif eofscaling == 1:
+ projected_pcs *= np.sqrt(self._eofs_rot_var[slicer])
+ elif eofscaling == 2:
+ projected_pcs *= self._eofs_rot_var[slicer]
+ else:
+ raise ValueError('invalid PC scaling option: '
+ '{!s}'.format(eofscaling))
+ if eof_ndim > input_ndim:
+ # If an extra dimension was introduced, remove it before returning
+ # the projected PCs.
+ projected_pcs = projected_pcs[0]
+ return projected_pcs
+
+ def _solverdata(self):
+ """Get the raw data from the EOF solver."""
+ return self._solver._data
+
+ def _to2d(self, eofs):
+ """Re-shape EOFs to 2D and remove missing values."""
+ # Re-shape to 2D.
+ info = {}
+ neofs = eofs.shape[0]
+ channels = eofs.shape[1:]
+ nspace = np.prod(channels)
+ ev = eofs.reshape([neofs, nspace])
+ # Remove missing values.
+ try:
+ ev = ev.filled(fill_value=np.nan)
+ info['filled'] = True
+ except AttributeError:
+ info['filled'] = False
+ nmi = np.where(np.isnan(ev[0]) == False)[0]
+ ev_nm = ev[:, nmi]
+ # Record information about the transform.
+ info['neofs'] = neofs
+ info['original_shape'] = channels
+ info['dtype'] = eofs.dtype
+ info['non_missing_locations'] = nmi
+ return ev_nm, info
+
+ def _from2d(self, ev, info):
+ """
+ Re-shape a 2D array to full shape and replace missing values.
+
+ """
+ channels = np.prod(info['original_shape'])
+ eofs = np.ones([info['neofs'], channels],
+ dtype=info['dtype']) * np.nan
+ eofs[:, info['non_missing_locations']] = ev
+ if info['filled']:
+ eofs = ma.array(eofs, mask=np.where(np.isnan(eofs), True, False))
+ eofs = eofs.reshape((info['neofs'],) + info['original_shape'])
+ return eofs
+
+ def apply_metadata(self, var, metadata):
+ """Convert an array to a metadata holding instance."""
+ return var
+
+ def strip_metadata(self, var):
+ """Convert a metadata holding instance to an array."""
+ return var, None
diff --git a/lib/eofs/tests/data/latitude_r.npy b/lib/eofs/tests/data/latitude_r.npy
new file mode 100644
index 0000000..2eecf31
Binary files /dev/null and b/lib/eofs/tests/data/latitude_r.npy differ
diff --git a/lib/eofs/tests/data/longitude_r.npy b/lib/eofs/tests/data/longitude_r.npy
new file mode 100644
index 0000000..142c2d9
Binary files /dev/null and b/lib/eofs/tests/data/longitude_r.npy differ
diff --git a/lib/eofs/tests/data/rotated_eigenvalues.scaled.npy b/lib/eofs/tests/data/rotated_eigenvalues.scaled.npy
new file mode 100644
index 0000000..903cdd1
Binary files /dev/null and b/lib/eofs/tests/data/rotated_eigenvalues.scaled.npy differ
diff --git a/lib/eofs/tests/data/rotated_eigenvalues.unscaled.npy b/lib/eofs/tests/data/rotated_eigenvalues.unscaled.npy
new file mode 100644
index 0000000..59988bf
Binary files /dev/null and b/lib/eofs/tests/data/rotated_eigenvalues.unscaled.npy differ
diff --git a/lib/eofs/tests/data/rotated_eofs.scaled.npy b/lib/eofs/tests/data/rotated_eofs.scaled.npy
new file mode 100644
index 0000000..b946707
Binary files /dev/null and b/lib/eofs/tests/data/rotated_eofs.scaled.npy differ
diff --git a/lib/eofs/tests/data/rotated_eofs.unscaled.npy b/lib/eofs/tests/data/rotated_eofs.unscaled.npy
new file mode 100644
index 0000000..d4b9fd1
Binary files /dev/null and b/lib/eofs/tests/data/rotated_eofs.unscaled.npy differ
diff --git a/lib/eofs/tests/data/rotated_pcs.scaled.npy b/lib/eofs/tests/data/rotated_pcs.scaled.npy
new file mode 100644
index 0000000..c6b23c2
Binary files /dev/null and b/lib/eofs/tests/data/rotated_pcs.scaled.npy differ
diff --git a/lib/eofs/tests/data/rotated_pcs.unscaled.npy b/lib/eofs/tests/data/rotated_pcs.unscaled.npy
new file mode 100644
index 0000000..ad4a8d0
Binary files /dev/null and b/lib/eofs/tests/data/rotated_pcs.unscaled.npy differ
diff --git a/lib/eofs/tests/data/rotated_variance.scaled.npy b/lib/eofs/tests/data/rotated_variance.scaled.npy
new file mode 100644
index 0000000..692b538
Binary files /dev/null and b/lib/eofs/tests/data/rotated_variance.scaled.npy differ
diff --git a/lib/eofs/tests/data/rotated_variance.unscaled.npy b/lib/eofs/tests/data/rotated_variance.unscaled.npy
new file mode 100644
index 0000000..684d658
Binary files /dev/null and b/lib/eofs/tests/data/rotated_variance.unscaled.npy differ
diff --git a/lib/eofs/tests/data/sst_r.npy b/lib/eofs/tests/data/sst_r.npy
new file mode 100644
index 0000000..30fbafa
Binary files /dev/null and b/lib/eofs/tests/data/sst_r.npy differ
diff --git a/lib/eofs/tests/data/time_r.npy b/lib/eofs/tests/data/time_r.npy
new file mode 100644
index 0000000..2981d65
Binary files /dev/null and b/lib/eofs/tests/data/time_r.npy differ
diff --git a/lib/eofs/tests/reference.py b/lib/eofs/tests/reference.py
index 0a45b03..5952d81 100644
--- a/lib/eofs/tests/reference.py
+++ b/lib/eofs/tests/reference.py
@@ -271,3 +271,122 @@ def reference_multivariate_solution(container_type, weights):
except TypeError:
solution[var] = None, None
return solution
+
+
+def _read_reference_rotated_solution(scaled):
+ if scaled:
+ ident = 'scaled'
+ else:
+ ident = 'unscaled'
+ field_names = ['time_r',
+ 'latitude_r',
+ 'longitude_r',
+ 'sst_r',
+ 'rotated_eigenvalues.{!s}'.format(ident),
+ 'rotated_eofs.{!s}'.format(ident),
+ 'rotated_pcs.{!s}'.format(ident),
+ 'rotated_variance.{!s}'.format(ident),]
+ fields = {name.split('.')[0]: _tomasked(_retrieve_test_field(name))
+ for name in field_names}
+ fields['sst_r'] -= fields['sst_r'].mean(axis=0)
+ fields['rotated_pcs'] = fields['rotated_pcs'].transpose()
+ return fields
+
+
+def reference_rotated_solution(container_type, scaled):
+ """Obtain a reference rotated EOF analysis solution.
+
+ **Arguments:**
+
+ *container_type*
+ Either 'standard', 'cdms', or 'iris'.
+
+ *scaled*
+ If *True* use scaled EOFs, if *False* use un-scaled EOFs.
+
+ """
+ container_type = container_type.lower()
+ if container_type not in ('standard', 'iris', 'cdms'):
+ raise ValueError("unknown container type "
+ "'{!s}'".format(container_type))
+ solution = _read_reference_rotated_solution(scaled)
+ time_units = 'days since 1800-1-1 00:00:00'
+ neofs = len(solution['rotated_eigenvalues'])
+ if container_type == 'cdms':
+ try:
+ time_dim = cdms2.createAxis(solution['time_r'], id='time')
+ time_dim.designateTime()
+ time_dim.units = time_units
+ lat_dim = cdms2.createAxis(solution['latitude_r'], id='latitude')
+ lat_dim.designateLatitude()
+ lon_dim = cdms2.createAxis(solution['longitude_r'], id='longitude')
+ lon_dim.designateLongitude()
+ eof_dim = cdms2.createAxis(np.arange(1, neofs+1), id='eof')
+ eof_dim.long_name = 'eof_number'
+ solution['sst_r'] = cdms2.createVariable(
+ solution['sst_r'],
+ axes=[time_dim, lat_dim, lon_dim],
+ id='sst')
+ solution['rotated_eigenvalues'] = cdms2.createVariable(
+ solution['rotated_eigenvalues'],
+ axes=[eof_dim],
+ id='eigenvalues')
+ solution['rotated_eofs'] = cdms2.createVariable(
+ solution['rotated_eofs'],
+ axes=[eof_dim, lat_dim, lon_dim],
+ id='eofs')
+ solution['rotated_pcs'] = cdms2.createVariable(
+ solution['rotated_pcs'],
+ axes=[time_dim, eof_dim],
+ id='pcs')
+ solution['rotated_variance'] = cdms2.createVariable(
+ solution['rotated_variance'],
+ axes=[eof_dim],
+ id='variance')
+ except NameError:
+ raise ValueError("cannot use container 'cdms' without the "
+ "cdms2 module")
+ elif container_type == 'iris':
+ try:
+ time_dim = DimCoord(solution['time_r'],
+ standard_name='time',
+ units=Unit(time_units, 'gregorian'))
+ lat_dim = DimCoord(solution['latitude_r'],
+ standard_name='latitude',
+ units='degrees_north')
+ lat_dim.guess_bounds()
+ lon_dim = DimCoord(solution['longitude_r'],
+ standard_name='longitude',
+ units='degrees_east')
+ lon_dim.guess_bounds()
+ eof_dim = DimCoord(np.arange(1, neofs+1),
+ long_name='eof')
+ solution['sst_r']= Cube(
+ solution['sst_r'],
+ dim_coords_and_dims=zip((time_dim, lat_dim, lon_dim),
+ range(3)),
+ long_name='sst')
+ solution['rotated_eigenvalues']= Cube(
+ solution['rotated_eigenvalues'],
+ dim_coords_and_dims=zip((eof_dim,),
+ range(1)),
+ long_name='rotated_eigenvalues')
+ solution['rotated_eofs']= Cube(
+ solution['rotated_eofs'],
+ dim_coords_and_dims=zip((eof_dim, lat_dim, lon_dim),
+ range(3)),
+ long_name='eofs')
+ solution['rotated_pcs']= Cube(
+ solution['rotated_pcs'],
+ dim_coords_and_dims=zip((time_dim, eof_dim),
+ range(2)),
+ long_name='pcs')
+ solution['rotated_variance']= Cube(
+ solution['rotated_variance'],
+ dim_coords_and_dims=zip((eof_dim,),
+ range(1)),
+ long_name='variance')
+ except NameError:
+ raise ValueError("cannot use container 'iris' without the "
+ "iris module")
+ return solution
diff --git a/lib/eofs/tests/test_rotation.py b/lib/eofs/tests/test_rotation.py
new file mode 100644
index 0000000..34af78e
--- /dev/null
+++ b/lib/eofs/tests/test_rotation.py
@@ -0,0 +1,200 @@
+"""Test EOF rotations against reference solutions."""
+# (c) Copyright 2013 Andrew Dawson. All Rights Reserved.
+#
+# This file is part of eofs.
+#
+# eofs is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# eofs is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+# for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with eofs. If not, see .
+from nose import SkipTest
+import numpy as np
+try:
+ from iris.cube import Cube
+except:
+ pass
+
+import eofs
+import eofs.experimental.rotation as rotation
+from eofs.tests import EofsTest
+
+from utils import sign_adjustments, error
+from reference import reference_rotated_solution
+
+
+# Create a mapping from interface name to solver class.
+solvers = {'standard': eofs.standard.Eof}
+try:
+ solvers['cdms'] = eofs.cdms.Eof
+except AttributeError:
+ pass
+try:
+ solvers['iris'] = eofs.iris.Eof
+except AttributeError:
+ pass
+
+# Create a mapping from interface name to rotator class.
+rotators = {'standard': rotation.standard.Rotator}
+try:
+ rotators['cdms'] = rotation.cdms.Rotator
+except AttributeError:
+ pass
+try:
+ rotators['iris'] = rotation.iris.Rotator
+except AttributeError:
+ pass
+
+
+class RotatorTest(EofsTest):
+ """Base class for all rotation test classes."""
+ interface = None
+ scaled = None
+ method = None
+ neofs = 5
+
+ @classmethod
+ def setup_class(cls):
+ try:
+ cls.solution = reference_rotated_solution(cls.interface,
+ cls.scaled)
+ except ValueError:
+ raise SkipTest('library component not available '
+ 'for {!s} interface'.format(cls.interface))
+ try:
+ # use default kws for solver, already well tested for weights etc.
+ solver = solvers[cls.interface](cls.solution['sst_r'], ddof=0)
+ cls.rotator = rotators[cls.interface](solver,
+ cls.neofs,
+ method=cls.method,
+ scaled=cls.scaled)
+ except KeyError:
+ raise SkipTest('library component not available '
+ 'for {!s} interface'.format(cls.interface))
+
+ def test_eofs(self):
+ # generate EOF tests for normalized and non-normalized EOFs
+ for renormalize in (False, True):
+ yield self.check_eofs, renormalize
+
+ def check_eofs(self, renormalize):
+ # rotated EOFs should match the (possibly normalized) reference
+ # solution
+ eofs = self._tomasked(
+ self.rotator.eofs(neofs=self.neofs, renormalize=renormalize))
+ reofs = self._tomasked(self.solution['rotated_eofs']).copy()
+ eofs *= sign_adjustments(eofs, reofs)
+ if renormalize:
+ variance = (reofs ** 2.).sum(axis=1).sum(axis=1)
+ reofs = reofs / np.sqrt(variance)[:, np.newaxis, np.newaxis]
+ self.assert_almost_equal(error(eofs, reofs), 0, places=3)
+
+ def test_pcs(self):
+ # generate PC tests for normalized and non-normalized PCs
+ for normalized in (False, True):
+ yield self.check_pcs, normalized
+
+ def check_pcs(self, normalized):
+ # rotated PCs should math the (possibly normalized) reference solution
+ pcs = self._tomasked(
+ self.rotator.pcs(npcs=self.neofs, normalized=normalized))
+ rpcs = self._tomasked(self.solution['rotated_pcs']).copy()
+ pcs *= sign_adjustments(pcs.transpose(), rpcs.transpose()).transpose()
+ if normalized:
+ rpcs /= rpcs.std(axis=0, ddof=1)
+ self.assert_almost_equal(error(pcs, rpcs), 0, places=3)
+
+ def test_varianceFraction(self):
+ # rotated variance fractions should match the reference solution
+ if not self.scaled:
+ # variance fraction is not meaningful when normalized EOFs are
+ # rotated
+ return
+ variance = self._tomasked(
+ self.rotator.varianceFraction(neigs=self.neofs)) * 100.
+ rvariance = self._tomasked(self.solution['rotated_variance'])
+ self.assert_array_almost_equal(variance, rvariance, decimal=3)
+
+
+#-----------------------------------------------------------------------------
+# Tests for the standard interface
+
+
+class StandardRotatorTest(RotatorTest):
+ """Base class for all standard interface solution test cases."""
+ interface = 'standard'
+
+ def _tomasked(self, value):
+ return value
+
+
+class TestStandardScaledVarimax(StandardRotatorTest):
+ """Rotation of scaled EOFs."""
+ scaled = True
+ method = 'varimax'
+
+
+class TestStandardUnscaledVarimax(StandardRotatorTest):
+ """Rotation of un-scaled EOFs."""
+ scaled = False
+ method = 'varimax'
+
+
+#-----------------------------------------------------------------------------
+# Tests for the cdms interface
+
+
+class CDMSRotatorTest(RotatorTest):
+ """Base class for all cdms interface solution test cases."""
+ interface = 'cdms'
+
+ def _tomasked(self, value):
+ try:
+ return value.asma()
+ except AttributeError:
+ return value
+
+
+class TestCDMSRotatorScaledVarimax(CDMSRotatorTest):
+ """Rotation of scaled EOFs."""
+ scaled = True
+ method = 'varimax'
+
+
+class TestCDMSRotatorUnscaledVarimax(CDMSRotatorTest):
+ """Rotation of un-scaled EOFs."""
+ scaled = False
+ method = 'varimax'
+
+
+#-----------------------------------------------------------------------------
+# Tests for the cdms interface
+
+
+class IrisRotatorTest(RotatorTest):
+ """Base class for all iris interface solution test cases."""
+ interface = 'iris'
+
+ def _tomasked(self, value):
+ if type(value) is not Cube:
+ return value
+ return value.data
+
+
+class TestIrisRotatorScaledVarimax(IrisRotatorTest):
+ """Rotation of scaled EOFs."""
+ scaled = True
+ method = 'varimax'
+
+
+class TestIrisRotatorUnscaledVarimax(IrisRotatorTest):
+ """Rotation of un-scaled EOFs."""
+ scaled = False
+ method = 'varimax'
diff --git a/setup.py b/setup.py
index 3084a9a..1b31225 100644
--- a/setup.py
+++ b/setup.py
@@ -24,6 +24,8 @@
packages = ['eofs',
'eofs.tools',
'eofs.multivariate',
+ 'eofs.experimental',
+ 'eofs.experimental.rotation',
'eofs.examples',
'eofs.tests']
diff --git a/tests.py b/tests.py
new file mode 100644
index 0000000..70242f4
--- /dev/null
+++ b/tests.py
@@ -0,0 +1,38 @@
+"""Run the test suite for `eofs`."""
+# (c) Copyright 2013 Andrew Dawson. All Rights Reserved.
+#
+# This file is part of eofs.
+#
+# eofs is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# eofs is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with eofs. If not, see .
+import nose
+
+
+#: Test modules to run by default.
+default_test_modules = [
+ 'eofs.tests.test_solution',
+ 'eofs.tests.test_error_handling',
+ 'eofs.tests.test_tools',
+ 'eofs.tests.test_multivariate_solution',
+ 'eofs.tests.test_multivariate_error_handling',
+ 'eofs.tests.test_rotation',
+]
+
+
+def run():
+ """Run the :mod:`eofs` tests."""
+ nose.main(defaultTest=default_test_modules)
+
+
+if __name__ == '__main__':
+ run()