diff --git a/meteoinfo-chart/src/main/java/org/meteoinfo/chart/geo/MapGridLine.java b/meteoinfo-chart/src/main/java/org/meteoinfo/chart/geo/MapGridLine.java index e86cb0d0..1d1a7668 100644 --- a/meteoinfo-chart/src/main/java/org/meteoinfo/chart/geo/MapGridLine.java +++ b/meteoinfo-chart/src/main/java/org/meteoinfo/chart/geo/MapGridLine.java @@ -24,8 +24,8 @@ import java.util.stream.Collectors; public class MapGridLine extends GridLine { protected ProjectionInfo projInfo; - protected Extent extent; - protected Extent lonLatExtent; + protected Extent extent = new Extent(-100, 100, -100, 100); + protected Extent lonLatExtent = new Extent(-180, 180, -90, 90); protected List longitudeLocations; protected List latitudeLocations; protected boolean fixLocations = false; @@ -476,7 +476,7 @@ public class MapGridLine extends GridLine { for (GridLabel aGL : tLabels) { if (!aGL.isBorder()) { if (aGL.isLongitude()) { - if (Math.abs(aGL.getCoord().X) < 1000 && Math.abs(aGL.getCoord().Y) < 1000) { + if (Math.abs(aGL.getCoord().X) < 100000 && Math.abs(aGL.getCoord().Y) < 100000) { continue; } diff --git a/meteoinfo-chart/src/main/java/org/meteoinfo/chart/geo/MapPlot.java b/meteoinfo-chart/src/main/java/org/meteoinfo/chart/geo/MapPlot.java index 70ca813d..67558c2f 100644 --- a/meteoinfo-chart/src/main/java/org/meteoinfo/chart/geo/MapPlot.java +++ b/meteoinfo-chart/src/main/java/org/meteoinfo/chart/geo/MapPlot.java @@ -757,12 +757,10 @@ public class MapPlot extends Plot2D implements IWebMapPanel { */ @Override public void setDrawExtent(Extent extent) { - if (!this.fixDrawExtent) { - super.setDrawExtent(extent); + super.setDrawExtent(extent); - if (!this.isLonLatMap()) { - ((MapGridLine) this.gridLine).setExtent(extent); - } + if (!this.isLonLatMap()) { + ((MapGridLine) this.gridLine).setExtent(extent); } } diff --git a/meteoinfo-common/src/main/java/org/meteoinfo/common/util/GlobalUtil.java b/meteoinfo-common/src/main/java/org/meteoinfo/common/util/GlobalUtil.java index 9951cd38..aa3c2bc5 100644 --- a/meteoinfo-common/src/main/java/org/meteoinfo/common/util/GlobalUtil.java +++ b/meteoinfo-common/src/main/java/org/meteoinfo/common/util/GlobalUtil.java @@ -67,7 +67,7 @@ import java.util.zip.ZipInputStream; public static String getVersion(){ String version = GlobalUtil.class.getPackage().getImplementationVersion(); if (version == null || version.equals("")) { - version = "3.8.1"; + version = "3.8.2"; } return version; } diff --git a/meteoinfo-lab/milconfig.xml b/meteoinfo-lab/milconfig.xml index 1b040270..6fd07f0d 100644 --- a/meteoinfo-lab/milconfig.xml +++ b/meteoinfo-lab/milconfig.xml @@ -1,33 +1,37 @@ - - - - - - - - + + + + + + + + + + - + + + - + @@ -36,5 +40,5 @@
- + diff --git a/meteoinfo-lab/pylib/mipylib/meteolib/__init__$py.class b/meteoinfo-lab/pylib/mipylib/meteolib/__init__$py.class index ed5111b9..90febd02 100644 Binary files a/meteoinfo-lab/pylib/mipylib/meteolib/__init__$py.class and b/meteoinfo-lab/pylib/mipylib/meteolib/__init__$py.class differ diff --git a/meteoinfo-lab/pylib/mipylib/meteolib/__init__.py b/meteoinfo-lab/pylib/mipylib/meteolib/__init__.py index 8c3b3e41..2af2a97b 100644 --- a/meteoinfo-lab/pylib/mipylib/meteolib/__init__.py +++ b/meteoinfo-lab/pylib/mipylib/meteolib/__init__.py @@ -3,8 +3,10 @@ from .wrf import * from . import constants from .calc import * from .interpolate import * +from ._eof import * __all__ = ['wrf','constants','meteo','calc','interpolate'] __all__ += meteo.__all__ __all__ += calc.__all__ -__all__ += interpolate.__all__ \ No newline at end of file +__all__ += interpolate.__all__ +__all__ += _eof.__all__ \ No newline at end of file diff --git a/meteoinfo-lab/pylib/mipylib/meteolib/_eof.py b/meteoinfo-lab/pylib/mipylib/meteolib/_eof.py new file mode 100644 index 00000000..b8333d1d --- /dev/null +++ b/meteoinfo-lab/pylib/mipylib/meteolib/_eof.py @@ -0,0 +1,740 @@ +import mipylib.numeric as np +import warnings + +__all__ = ['EOF','varimax'] + + +class EOF(object): + """EOF analysis""" + + def __init__(self, dataset, weights=None, center=True, ddof=1): + """ + The EOF solution is computed at initialization time. Method + calls are used to retrieve computed quantities. + + **Arguments:** + + *dataset* + A `NDArray` with two or more dimensions containing the data to be analysed. + The first dimension is assumed to represent time. Missing + values are permitted, in the form of `nan` values. Missing values must be constant + with time (e.g., values of an oceanographic field over land). + + **Optional arguments:** + + *weights* + An array of weights whose shape is compatible with those of + the input array *dataset*. The weights can have the same + shape as *dataset* or a shape compatible with an array + broadcast (i.e., the shape of the weights can can match the + rightmost parts of the shape of the input array *dataset*). + If the input array *dataset* does not require weighting then + the value *None* may be used. Defaults to *None* (no + weighting). + + *center* + If *True*, the mean along the first axis of *dataset* (the + time-mean) will be removed prior to analysis. If *False*, + the mean along the first axis will not be removed. Defaults + to *True* (mean is removed). + + The covariance interpretation relies on the input data being + anomaly data with a time-mean of 0. Therefore, this option + should usually be set to *True*. Setting this option to + *True* has the useful side effect of propagating missing + values along the time dimension, ensuring that a solution + can be found even if missing values occur in different + locations at different times. + + *ddof* + 'Delta degrees of freedom'. The divisor used to normalize + the covariance matrix is *N - ddof* where *N* is the + number of samples. Defaults to *1*. + + **Returns:** + + *solver* + An `Eof` instance. + """ + # Store the input data in an instance variable. + if dataset.ndim < 2: + raise ValueError('the input data set must be ' + 'at least two dimensional') + self._data = dataset.copy() + + # Store information about the shape/size of the input data. + self._records = self._data.shape[0] + self._originalshape = self._data.shape[1:] + channels = np.prod(self._originalshape) + + # Weight the data set according to weighting argument. + if weights is not None: + try: + # The broadcast_arrays call returns a list, so the second index + # is retained, but also we want to remove the time dimension + # from the weights so the first index from the broadcast + # array is taken. + self._weights = np.broadcast_arrays( + self._data[0:1], weights)[1][0] + self._data = self._data * self._weights + except ValueError: + raise ValueError('weight array dimensions are incompatible') + except TypeError: + raise TypeError('weights are not a valid type') + else: + self._weights = None + + # Remove the time mean of the input data unless explicitly told + # not to by the "center" argument. + self._centered = center + if center: + self._data = self._center(self._data) + + # Reshape to two dimensions (time, space) creating the design matrix. + self._data = self._data.reshape([self._records, channels]) + + # Find the indices of values that are not missing in one row. All the + # rows will have missing values in the same places provided the + # array was centered. If it wasn't then it is possible that some + # missing values will be missed and the singular value decomposition + # will produce not a number for everything. + if not self._valid_nan(self._data): + raise ValueError('missing values detected in different ' + 'locations at different times') + nonMissingIndex = np.where(np.logical_not(np.isnan(self._data[0])))[0] + + # Remove missing values from the design matrix. + dataNoMissing = self._data[:, nonMissingIndex] + if dataNoMissing.size == 0: + raise ValueError('all input data is missing') + + # Compute the singular value decomposition of the design matrix. + try: + A, Lh, E = np.linalg.svd(dataNoMissing, full_matrices=False) + except (np.linalg.LinAlgError, ValueError): + raise ValueError('error encountered in SVD, check that missing ' + 'values are in the same places at each time and ' + 'that all the values are not missing') + + # Singular values are the square-root of the eigenvalues of the + # covariance matrix. Construct the eigenvalues appropriately and + # normalize by N-ddof where N is the number of observations. This + # corresponds to the eigenvalues of the normalized covariance matrix. + self._ddof = ddof + normfactor = float(self._records - self._ddof) + self._L = Lh * Lh / normfactor + + # Store the number of eigenvalues (and hence EOFs) that were actually + # computed. + self.neofs = len(self._L) + + # Re-introduce missing values into the eigenvectors in the same places + # as they exist in the input maps. Create an array of not-a-numbers + # and then introduce data values where required. We have to use the + # astype method to ensure the eigenvectors are the same type as the + # input dataset since multiplication by np.NaN will promote to 64-bit. + self._flatE = np.ones([self.neofs, channels], + dtype=self._data.dtype) * np.nan + self._flatE = self._flatE.astype(self._data.dtype) + self._flatE[:, nonMissingIndex] = E + + # Remove the scaling on the principal component time-series that is + # implicitly introduced by using SVD instead of eigen-decomposition. + # The PCs may be re-scaled later if required. + self._P = A * Lh + + def _center(self, in_array): + """Remove the mean of an array along the first dimension.""" + # Compute the mean along the first dimension. + mean = in_array.mean(axis=0) + + # Return the input array with its mean along the first dimension + # removed. + return (in_array - mean) + + def _valid_nan(self, in_array): + inan = np.isnan(in_array) + return (inan.any(axis=0) == inan.all(axis=0)).all() + + def pcs(self, pcscaling=0, npcs=None): + """Principal component time series (PCs). + + **Optional arguments:** + + *pcscaling* + Set the scaling of the retrieved PCs. The following + values are accepted: + + * *0* : Un-scaled PCs (default). + * *1* : PCs are scaled to unit variance (divided by the + square-root of their eigenvalue). + * *2* : PCs are multiplied by the square-root of their + eigenvalue. + + *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. + + **Returns:** + + *pcs* + An array where the columns are the ordered PCs. + + **Examples:** + + All un-scaled PCs:: + + pcs = solver.pcs() + + First 3 PCs scaled to unit variance:: + + pcs = solver.pcs(npcs=3, pcscaling=1) + + """ + slicer = slice(0, npcs) + if pcscaling == 0: + # Do not scale. + return self._P[:, slicer].copy() + elif pcscaling == 1: + # Divide by the square-root of the eigenvalue. + return self._P[:, slicer] / np.sqrt(self._L[slicer]) + elif pcscaling == 2: + # Multiply by the square root of the eigenvalue. + return self._P[:, slicer] * np.sqrt(self._L[slicer]) + else: + raise ValueError('invalid PC scaling option: ' + '{!s}'.format(pcscaling)) + + def eofs(self, eofscaling=0, neofs=None): + """Empirical orthogonal functions (EOFs). + + **Optional arguments:** + + *eofscaling* + Sets the scaling of the EOFs. The following values are + accepted: + + * *0* : Un-scaled EOFs (default). + * *1* : EOFs are divided by the square-root of their + eigenvalues. + * *2* : EOFs are multiplied by the square-root of their + eigenvalues. + + *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 with no scaling:: + + eofs = solver.eofs() + + The leading EOF with scaling applied:: + + eof1 = solver.eofs(neofs=1, eofscaling=1) + + """ + if neofs is None or neofs > self.neofs: + neofs = self.neofs + slicer = slice(0, neofs) + neofs = neofs or self.neofs + flat_eofs = self._flatE[slicer].copy() + if eofscaling == 0: + # No modification. A copy needs to be returned in case it is + # modified. If no copy is made the internally stored eigenvectors + # could be modified unintentionally. + rval = flat_eofs + elif eofscaling == 1: + # Divide by the square-root of the eigenvalues. + rval = flat_eofs / np.sqrt(self._L[slicer])[:, np.newaxis] + elif eofscaling == 2: + # Multiply by the square-root of the eigenvalues. + rval = flat_eofs * np.sqrt(self._L[slicer])[:, np.newaxis] + else: + raise ValueError('invalid eof scaling option: ' + '{!s}'.format(eofscaling)) + + return rval.reshape((neofs,) + self._originalshape) + + def eofs_correlation(self, neofs=None): + """Correlation map EOFs. + + Empirical orthogonal functions (EOFs) expressed as the + correlation between the 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 = solver.eofsAsCorrelation() + + The leading EOF:: + + eof1 = solver.eofsAsCorrelation(neofs=1) + + """ + # Retrieve the specified number of PCs. + pcs = self.pcs(npcs=neofs, pcscaling=1) + # Compute the correlation of the PCs with the input field. + c = correlation_map( + pcs, + self._data.reshape((self._records,) + self._originalshape)) + + return c + + def eofs_covariance(self, neofs=None, pcscaling=1): + """Covariance map EOFs. + + Empirical orthogonal functions (EOFs) expressed as the + covariance between the 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: + + * *0* : Un-scaled PCs. + * *1* : PCs are scaled to unit variance (divided by the + square-root of their eigenvalue) (default). + * *2* : PCs are multiplied by the square-root of their + eigenvalue. + + The default is to divide PCs by the square-root of their + eigenvalue so that the PCs are scaled to unit variance + (option 1). + + **Returns:** + + *eofs* + An array with the ordered EOFs along the first dimension. + + **Examples:** + + All EOFs:: + + eofs = solver.eofsAsCovariance() + + The leading EOF:: + + eof1 = solver.eofsAsCovariance(neofs=1) + + The leading EOF using un-scaled PCs:: + + eof1 = solver.eofsAsCovariance(neofs=1, pcscaling=0) + + """ + pcs = self.pcs(npcs=neofs, pcscaling=pcscaling) + # Divide the input data by the weighting (if any) before computing + # the covariance maps. + data = self._data.reshape((self._records,) + self._originalshape) + if self._weights is not None: + with warnings.catch_warnings(): + warnings.simplefilter('ignore', RuntimeWarning) + data /= self._weights + c = covariance_map(pcs, data, ddof=self._ddof) + + return c + + def eigenvalues(self, neigs=None): + """Eigenvalues (decreasing variances) associated with each EOF. + + **Optional argument:** + + *neigs* + Number of eigenvalues to return. Defaults to all + eigenvalues. If the number of eigenvalues requested is more + than the number that are available, then all available + eigenvalues will be returned. + + **Returns:** + + *eigenvalues* + An array containing the eigenvalues arranged largest to + smallest. + + **Examples:** + + All eigenvalues:: + + eigenvalues = solver.eigenvalues() + + The first eigenvalue:: + + eigenvalue1 = solver.eigenvalues(neigs=1) + + """ + # Create a slicer and use it on the eigenvalue array. A copy must be + # returned in case the slicer takes all elements, in which case a + # reference to the eigenvalue array is returned. If this is modified + # then the internal eigenvalues array would then be modified as well. + slicer = slice(0, neigs) + return self._L[slicer].copy() + + def variance_fraction(self, neigs=None): + """Fractional EOF mode variances. + + The fraction of the total variance explained by each EOF mode, + values between 0 and 1 inclusive. + + **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 EOF mode:: + + variance_fractions = solver.varianceFraction() + + The fractional variance represented by the first EOF mode:: + + variance_fraction_mode_1 = solver.VarianceFraction(neigs=1) + + """ + # Return the array of eigenvalues divided by the sum of the + # eigenvalues. + slicer = slice(0, neigs) + return self._L[slicer] / self._L.sum() + + def total_anomaly_variance(self): + """ + Total variance associated with the field of anomalies (the sum + of the eigenvalues). + + **Returns:** + + *total_variance* + A scalar value. + + **Example:** + + Get the total variance:: + + total_variance = solver.totalAnomalyVariance() + + """ + # Return the sum of the eigenvalues. + return self._L.sum() + + def northtest(self, neigs=None, vfscaled=False): + """Typical errors for eigenvalues. + + The method of North et al. (1982) is used to compute the typical + error for each eigenvalue. It is assumed that the number of + times in the input data set is the same as the number of + independent realizations. If this assumption is not valid then + the result may be inappropriate. + + **Optional arguments:** + + *neigs* + The number of eigenvalues to return typical errors for. + Defaults to typical errors for all eigenvalues. If the + number of eigenvalues requested is more than the number that + are available, then typical errors for all available + eigenvalues will be returned. + + *vfscaled* + If *True* scale the errors by the sum of the eigenvalues. + This yields typical errors with the same scale as the values + returned by `Eof.varianceFraction`. If *False* then no + scaling is done. Defaults to *False* (no scaling). + + **Returns:** + + *errors* + An array containing the typical errors. + + **References** + + North G.R., T.L. Bell, R.F. Cahalan, and F.J. Moeng (1982) + Sampling errors in the estimation of empirical orthogonal + functions. *Mon. Weather. Rev.*, **110**, pp 669-706. + + **Examples:** + + Typical errors for all eigenvalues:: + + errors = solver.northTest() + + Typical errors for the first 5 eigenvalues scaled by the sum of + the eigenvalues:: + + errors = solver.northTest(neigs=5, vfscaled=True) + + """ + slicer = slice(0, neigs) + # Compute the factor that multiplies the eigenvalues. The number of + # records is assumed to be the number of realizations of the field. + factor = np.sqrt(2.0 / self._records) + # If requested, allow for scaling of the eigenvalues by the total + # variance (sum of the eigenvalues). + if vfscaled: + factor /= self._L.sum() + # Return the typical errors. + return self._L[slicer] * factor + + def varimax(self, eofs, **kwargs): + """Rotation empirical orthogonal functions (REOFs). + + **Optional arguments:** + + *eofs* + EOF data array. + + *kwargs* + Parameters for varimax rotation computation. + + **Returns:** + + *reofs* + An array with the ordered REOFs along the first dimension. + """ + neofs = eofs.shape[0] + channels = np.prod(self._originalshape) + eofs = eofs.reshape(neofs, channels) + eofs = eofs.T + + reofs = varimax(eofs, kwargs)[0] + reofs = reofs.T + reofs = reofs.reshape((neofs,) + self._originalshape) + + return reofs + + +def _check_flat_center(pcs, field): + """ + Check PCs and a field for shape compatibility, flatten both to 2D, + and center along the first dimension. + + This set of operations is common to both covariance and correlation + calculations. + + """ + # Get the number of times in the field. + records = field.shape[0] + if records != pcs.shape[0]: + # Raise an error if the field has a different number of times to the + # PCs provided. + raise ValueError("PCs and field must have the same first dimension") + if len(pcs.shape) > 2: + # Raise an error if the PCs are more than 2D. + raise ValueError("PCs must be 1D or 2D") + # Check if the field is 1D. + if len(field.shape) == 1: + originalshape = tuple() + channels = 1 + else: + # Record the shape of the field and the number of spatial elements. + originalshape = field.shape[1:] + channels = np.prod(originalshape) + # Record the number of PCs. + if len(pcs.shape) == 1: + npcs = 1 + npcs_out = tuple() + else: + npcs = pcs.shape[1] + npcs_out = (npcs,) + # Create a flattened field so iterating over space is simple. Also do this + # for the PCs to ensure they are 2D. + field_flat = field.reshape([records, channels]) + pcs_flat = pcs.reshape([records, npcs]) + # Centre both the field and PCs in the time dimension. + field_flat = field_flat - field_flat.mean(axis=0) + pcs_flat = pcs_flat - pcs_flat.mean(axis=0) + return pcs_flat, field_flat, npcs_out + originalshape + + +def correlation_map(pcs, field): + """Correlation maps for a set of PCs and a spatial-temporal field. + + Given an array where the columns are PCs (e.g., as output from + `~eofs.standard.Eof.pcs`) and an array containing spatial-temporal + data where the first dimension represents time, one correlation map + per PC is computed. + + The field must have the same temporal dimension as the PCs. Any + number of spatial dimensions (including zero) are allowed in the + field and there can be any number of PCs. + + **Arguments:** + + *pcs* + PCs as the columns of an array. + + *field* + Spatial-temporal field with time as the first dimension. + + **Returns:** + + *correlation_maps* + An array with the correlation maps along the first dimension. + + **Example:** + + Compute correlation maps for each PC:: + + pcs = solver.pcs(pcscaling=1) + correlation_maps = correlation_maps(pcs, field) + + """ + # Check PCs and fields for validity, flatten the arrays ready for the + # computation and remove the mean along the leading dimension. + pcs_cent, field_cent, out_shape = _check_flat_center(pcs, field) + # Compute the standard deviation of the PCs and the fields along the time + # dimension (the leading dimension). + pcs_std = pcs_cent.std(axis=0) + field_std = field_cent.std(axis=0) + # Set the divisor. + div = float(pcs_cent.shape[0]) + # Compute the correlation map. + cor = np.dot(field_cent.T, pcs_cent).T / div + cor /= np.outer(pcs_std, field_std) + # Return the correlation with the appropriate shape. + return cor.reshape(out_shape) + +def covariance_map(pcs, field, ddof=1): + """Covariance maps for a set of PCs and a spatial-temporal field. + + Given an array where the columns are PCs (e.g., as output from + `eofs.standard.Eof.pcs`) and an array containing spatial-temporal + data where the first dimension represents time, one covariance map + per PC is computed. + + The field must have the same temporal dimension as the PCs. Any + number of spatial dimensions (including zero) are allowed in the + field and there can be any number of PCs. + + **Arguments:** + + *pcs* + PCs as the columns of an array. + + *field* + Spatial-temporal field with time as the first dimension. + + **Optional arguments:** + + *ddof* + 'Delta degrees of freedom'. The divisor used to normalize + the covariance matrix is *N - ddof* where *N* is the + number of samples. Defaults to *1*. + + **Returns:** + + *covariance_maps* + An array with the covariance maps along the first dimension. + + **Example:** + + Compute covariance maps for each PC:: + + pcs = solver.pcs(pcscaling=1) + covariance_maps = covariance_maps(pcs, field) + + """ + # Check PCs and fields for validity, flatten the arrays ready for the + # computation and remove the mean along the leading dimension. + pcs_cent, field_cent, out_shape = _check_flat_center(pcs, field) + # Set the divisor according to the specified delta-degrees of freedom. + div = float(pcs_cent.shape[0] - ddof) + # Compute the covariance map, making sure it has the appropriate shape. + cov = (np.dot(field_cent.T, pcs_cent).T / div).reshape(out_shape) + return cov + +def varimax(x, norm=True, tol=1e-10, it_max=1000): + """ + Rotate EOFs according to varimax algorithm + + :param x: (*array_like*) Input 2-D array. + :param norm: (*boolean*) Determines whether to do Kaiser normalization the rows + of the loadings before performing the rotation. Default is `True`. + :param tol: (*float*) Tolerance. + :param it_max: (*int*) Specifies the maximum number of iterations to do. + + :returns: Rotated EOFs and rotate matrix. + """ + has_nan = False + if x.contains_nan(): #Has NaN value + mask = np.isnan(x).sum(axis=1) + valid_idx = np.where(mask==0)[0] + xx = x[valid_idx,:] + has_nan = True + else: + xx = x.copy() + + if norm: + h = np.sqrt(np.sum(xx**2, axis=1)) + xx = xx / h[:, None] + + p, nc = xx.shape + TT = np.eye(nc) + d = 0 + for _ in range(it_max): + z = np.dot(xx, TT) + B = np.dot(xx.T, (z**3 - np.dot(z, np.diag(np.squeeze(np.dot(np.ones((1,p)), (z**2))))) / p)) + U, S, Vh = np.linalg.svd(B) + TT = np.dot(U, Vh) + d2 = d + d = np.sum(S) + # End if exceeded tolerance. + if d < d2 * (1 + tol): + break + + # Final matrix. + r = np.dot(xx, TT) + + if norm: + r = r * h[:,None] + + if has_nan: + rr = np.ones(x.shape) * np.nan + rr[valid_idx,:] = r + r = rr + + return r, TT diff --git a/meteoinfo-lab/pylib/mipylib/meteolib/meteo$py.class b/meteoinfo-lab/pylib/mipylib/meteolib/meteo$py.class index 659c4198..bb252a00 100644 Binary files a/meteoinfo-lab/pylib/mipylib/meteolib/meteo$py.class and b/meteoinfo-lab/pylib/mipylib/meteolib/meteo$py.class differ diff --git a/meteoinfo-lab/pylib/mipylib/meteolib/meteo.py b/meteoinfo-lab/pylib/mipylib/meteolib/meteo.py index cb449912..45f951c9 100644 --- a/meteoinfo-lab/pylib/mipylib/meteolib/meteo.py +++ b/meteoinfo-lab/pylib/mipylib/meteolib/meteo.py @@ -19,7 +19,7 @@ __all__ = [ 'moist_lapse','p2h','qair2rh','rh2dewpoint', 'sigma_to_pressure','tc2tf', 'tf2tc','uv2ds','pressure_to_height_std', - 'height_to_pressure_std','eof','vapor_pressure','varimax' + 'height_to_pressure_std','eof','vapor_pressure' ] def uv2ds(u, v): diff --git a/meteoinfo-lab/pylib/mipylib/numeric/core/_ndarray.py b/meteoinfo-lab/pylib/mipylib/numeric/core/_ndarray.py index f5a96919..97d95ef5 100644 --- a/meteoinfo-lab/pylib/mipylib/numeric/core/_ndarray.py +++ b/meteoinfo-lab/pylib/mipylib/numeric/core/_ndarray.py @@ -100,8 +100,12 @@ class NDArray(object): #deal with Ellipsis if Ellipsis in indices: + n = 0 + for ii in indices: + if ii is not None: + n += 1; + indices1 = [] - n = self.ndim - len(indices) + 1 for ii in indices: if ii is Ellipsis: for _ in range(n): @@ -110,7 +114,7 @@ class NDArray(object): indices1.append(ii) indices = indices1 -#for all int indices + #for all int indices if len(indices) == self.ndim: allint = True aindex = self._array.getIndex() @@ -978,6 +982,40 @@ class NDArray(object): r = ArrayMath.round(self._array, decimals) return NDArray(r) + def clip(self, min=None, max=None): + """ + Clip (limit) the values in an array. + + Given an interval, values outside the interval are clipped to the interval edges. For example, + if an interval of [0, 1] is specified, values smaller than 0 become 0, and values larger than 1 + become 1. + + Parameters + ---------- + min, max : array_like or None + Minimum and maximum value. If ``None``, clipping is not performed on + the corresponding edge. Only one of `min` and `max` may be + ``None``. Both are broadcast against this array. + + Returns + ------- + clipped_array : NDArray + An array with the elements of this array, but where values + < `min` are replaced with `min`, and those > `max` + with `max`. + """ + if min is None: + if max is None: + raise ValueError("Only one of min and max my be None!") + + r = ArrayMath.clipMax(self._array, max) + elif max is None: + r = ArrayMath.clipMin(self._array, min) + else: + r = ArrayMath.clip(self._array, min, max) + + return self.array_wrap(r) + def rot90(self, k=1): """ Rotate an array by 90 degrees in the counter-clockwise direction. The first two dimensions diff --git a/meteoinfo-lab/pylib/mipylib/numeric/core/numeric$py.class b/meteoinfo-lab/pylib/mipylib/numeric/core/numeric$py.class index 7725279e..fc8775d7 100644 Binary files a/meteoinfo-lab/pylib/mipylib/numeric/core/numeric$py.class and b/meteoinfo-lab/pylib/mipylib/numeric/core/numeric$py.class differ diff --git a/meteoinfo-lab/pylib/mipylib/numeric/core/numeric.py b/meteoinfo-lab/pylib/mipylib/numeric/core/numeric.py index 6738532c..03ffab08 100644 --- a/meteoinfo-lab/pylib/mipylib/numeric/core/numeric.py +++ b/meteoinfo-lab/pylib/mipylib/numeric/core/numeric.py @@ -616,6 +616,37 @@ def round(x, decimals=0): else: return __builtin__.round(x, decimals) +def clip(a, a_min, a_max): + """ + Clip (limit) the values in an array. + + Given an interval, values outside the interval are clipped to + the interval edges. For example, if an interval of ``[0, 1]`` + is specified, values smaller than 0 become 0, and values larger + than 1 become 1. + + Equivalent to but faster than ``np.minimum(a_max, np.maximum(a, a_min))``. + + No check is performed to ensure ``a_min < a_max``. + + Parameters + ---------- + a : array_like + Array containing elements to clip. + a_min, a_max : array_like or None + Minimum and maximum value. If ``None``, clipping is not performed on + the corresponding edge. Only one of `a_min` and `a_max` may be + ``None``. Both are broadcast against `a`. + + Returns + ------- + clipped_array : ndarray + An array with the elements of `a`, but where values + < `a_min` are replaced with `a_min`, and those > `a_max` + with `a_max`. + """ + return a.clip(a_min, a_max) + def square(x): """ Return the element-wise square of the input. diff --git a/meteoinfo-lab/pylib/mipylib/plotlib/_axes$py.class b/meteoinfo-lab/pylib/mipylib/plotlib/_axes$py.class index 8fe2a976..0b65cf29 100644 Binary files a/meteoinfo-lab/pylib/mipylib/plotlib/_axes$py.class and b/meteoinfo-lab/pylib/mipylib/plotlib/_axes$py.class differ diff --git a/meteoinfo-lab/pylib/mipylib/plotlib/_axes.py b/meteoinfo-lab/pylib/mipylib/plotlib/_axes.py index d64d1e03..4e315284 100644 --- a/meteoinfo-lab/pylib/mipylib/plotlib/_axes.py +++ b/meteoinfo-lab/pylib/mipylib/plotlib/_axes.py @@ -747,7 +747,6 @@ class Axes(object): extent = self._axes.getDrawExtent() extent.minX = xmin extent.maxX = xmax - self._axes.setFixDrawExtent(False) self._axes.setDrawExtent(extent) self._axes.setExtent(extent.clone()) self._axes.setFixDrawExtent(True) @@ -776,11 +775,20 @@ class Axes(object): extent = self._axes.getDrawExtent() extent.minY = ymin extent.maxY = ymax - self._axes.setFixDrawExtent(False) self._axes.setDrawExtent(extent) self._axes.setExtent(extent.clone()) self._axes.setFixDrawExtent(True) + def set_draw_extent(self, extent): + """ + Set axes draw extent. + + :param extent: The extent + """ + if not self._axes.isFixDrawExtent(): + self._axes.setDrawExtent(extent.clone()) + self._axes.setExtent(extent.clone()) + def twinx(self): """ Make a second axes that shares the x-axis. The new axes will overlay *ax*. The ticks @@ -1507,6 +1515,7 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(graphics, zorder=zorder) + return graphics def scatter(self, *args, **kwargs): @@ -1626,7 +1635,6 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(graphics, zorder=zorder) - self._axes.setAutoExtent() return graphics @@ -1879,7 +1887,6 @@ class Axes(object): yerrU, line, eline, capsize) zorder = kwargs.pop('zorder', None) self.add_graphic(graphics, zorder=zorder) - self._axes.setAutoExtent() return graphics @@ -2008,7 +2015,6 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(graphics, zorder=zorder) - self._axes.setAutoExtent() if autowidth: barswidth = kwargs.pop('barswidth', 0.8) self._axes.setBarsWidth(barswidth) @@ -2135,7 +2141,7 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(graphics, zorder=zorder) - self._axes.setAutoExtent() + if autoheight: barsheight = kwargs.pop('barsheight', 0.8) self._axes.setBarsWidth(barsheight) @@ -2300,7 +2306,6 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(graphics, zorder=zorder) - self._axes.setAutoExtent() return linefmt @@ -2378,8 +2383,7 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(graphics, zorder=zorder) - self._axes.setExtent(graphics.getExtent()) - self._axes.setDrawExtent(graphics.getExtent()) + self.set_draw_extent(graphics.getExtent()) return graphics @@ -2538,8 +2542,7 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(graphics, zorder=zorder) - self._axes.setExtent(graphics.getExtent()) - self._axes.setDrawExtent(graphics.getExtent()) + self.set_draw_extent(graphics.getExtent()) return graphics @@ -2640,6 +2643,7 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(igraphic, zorder=zorder) + self.set_draw_extent(igraphic.getExtent()) gridline = self._axes.getGridLine() gridline.setTop(True) @@ -2698,8 +2702,7 @@ class Axes(object): if visible: zorder = kwargs.pop('zorder', None) self.add_graphic(graphics, zorder=zorder) - self._axes.setExtent(graphics.getExtent()) - self._axes.setDrawExtent(graphics.getExtent()) + self.set_draw_extent(graphics.getExtent()) return graphics @@ -2746,8 +2749,8 @@ class Axes(object): if visible: zorder = kwargs.pop('zorder', None) self.add_graphic(graphics, zorder=zorder) - self._axes.setExtent(graphics.getExtent()) - self._axes.setDrawExtent(graphics.getExtent()) + self.set_draw_extent(graphics.getExtent()) + return graphics def text(self, x, y, s, **kwargs): @@ -2790,7 +2793,6 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(graphic, zorder=zorder) - self._axes.setAutoExtent() return graphic @@ -2820,7 +2822,6 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(graphic, zorder=zorder) - self._axes.setAutoExtent() return graphic @@ -2865,7 +2866,6 @@ class Axes(object): dy = dy - sy * 2 graphic = GraphicFactory.createArrowLine(x0, y0, dx, dy, alb) self.add_graphic(graphic) - self._axes.setAutoExtent() return ctext, graphic @@ -2911,7 +2911,7 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(graphics, zorder=zorder) - self._axes.setAutoExtent() + return graphics def patch(self, x, y=None, **kwargs): @@ -2932,7 +2932,7 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(graphics, zorder=zorder) - self._axes.setAutoExtent() + return graphics def rectangle(self, position, curvature=None, **kwargs): @@ -2949,7 +2949,7 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(graphic, zorder=zorder) - self._axes.setAutoExtent() + return graphic def fill_between(self, x, y1, y2=0, where=None, **kwargs): @@ -2997,7 +2997,6 @@ class Axes(object): graphics = GraphicFactory.createFillBetweenPolygons(xdata, y1, y2, where, pb) zorder = kwargs.pop('zorder', None) self.add_graphic(graphics, zorder=zorder) - self._axes.setAutoExtent() return pb @@ -3046,7 +3045,6 @@ class Axes(object): graphics = GraphicFactory.createFillBetweenPolygonsX(ydata, x1, x2, where, pb) zorder = kwargs.pop('zorder', None) self.add_graphic(graphics, zorder=zorder) - self._axes.setAutoExtent() return pb @@ -3127,7 +3125,6 @@ class Axes(object): for graphic in graphics: self.add_graphic(graphic) - self._axes.setAutoExtent() self._axes.setAspectType(AspectType.EQUAL) self._axes.getAxis(Location.BOTTOM).setVisible(False) self._axes.getAxis(Location.LEFT).setVisible(False) @@ -3273,7 +3270,6 @@ class Axes(object): graphics.setAntiAlias(antialias) self.add_graphic(graphics) - self._axes.setAutoExtent() return graphics @@ -3569,7 +3565,6 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(graphics, zorder=zorder) - self._axes.setAutoExtent() return graphics @@ -3675,7 +3670,6 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(igraphic, zorder=zorder) - self._axes.setAutoExtent() return igraphic @@ -3830,7 +3824,6 @@ class Axes(object): zorder = kwargs.pop('zorder', None) self.add_graphic(igraphic, zorder=zorder) - self._axes.setAutoExtent() return igraphic diff --git a/meteoinfo-lab/pylib/mipylib/plotlib/_mapaxes$py.class b/meteoinfo-lab/pylib/mipylib/plotlib/_mapaxes$py.class index 29e6cae3..ae19d051 100644 Binary files a/meteoinfo-lab/pylib/mipylib/plotlib/_mapaxes$py.class and b/meteoinfo-lab/pylib/mipylib/plotlib/_mapaxes$py.class differ diff --git a/meteoinfo-lab/pylib/mipylib/plotlib/_mapaxes.py b/meteoinfo-lab/pylib/mipylib/plotlib/_mapaxes.py index 8c31158a..1066cc81 100644 --- a/meteoinfo-lab/pylib/mipylib/plotlib/_mapaxes.py +++ b/meteoinfo-lab/pylib/mipylib/plotlib/_mapaxes.py @@ -850,8 +850,7 @@ class MapAxes(Axes): if visible: zorder = kwargs.pop('zorder', None) graphics = self.add_graphic(graphics, projection=proj, zorder=zorder) - self._axes.setExtent(graphics.getExtent()) - self._axes.setDrawExtent(graphics.getExtent()) + self.set_draw_extent(graphics.getExtent()) return graphics @@ -931,8 +930,7 @@ class MapAxes(Axes): if visible: zorder = kwargs.pop('zorder', None) contours = self.add_graphic(graphics, projection=proj, zorder=zorder) - self._axes.setDrawExtent(graphics.getExtent()) - self._axes.setExtent(graphics.getExtent()) + self.set_draw_extent(graphics.getExtent()) return graphics @@ -996,8 +994,7 @@ class MapAxes(Axes): if zorder is None: zorder = 0 graphics = self.add_graphic(graphics, projection=proj, zorder=zorder) - self._axes.setDrawExtent(graphics.getExtent()) - self._axes.setExtent(graphics.getExtent()) + self.set_draw_extent(graphics.getExtent()) return graphics @@ -1133,8 +1130,7 @@ class MapAxes(Axes): if zorder is None: zorder = 0 igraphic = self.add_graphic(igraphic, zorder=zorder) - self._axes.setDrawExtent(igraphic.getExtent().clone()) - self._axes.setExtent(igraphic.getExtent().clone()) + self.set_draw_extent(igraphic.getExtent()) gridline = self._axes.getGridLine() gridline.setTop(True) @@ -1152,7 +1148,7 @@ class MapAxes(Axes): :param cmap: (*string*) Color map string. :param colors: (*list*) If None (default), the colormap specified by cmap will be used. If a string, like ‘r’ or ‘red’, all levels will be plotted in this color. If a tuple of matplotlib - color args (string, float, rgb, etc), different levels will be plotted in different colors in + color args (string, float, rgb, etc.), different levels will be plotted in different colors in the order specified. :param fill_value: (*float*) Fill_value. Default is ``-9999.0``. :param proj: (*ProjectionInfo*) Map projection of the data. Default is None. @@ -1199,8 +1195,7 @@ class MapAxes(Axes): if visible: zorder = kwargs.pop('zorder', None) graphics = self.add_graphic(graphics, projection=proj, zorder=zorder) - self._axes.setExtent(graphics.getExtent()) - self._axes.setDrawExtent(graphics.getExtent()) + self.set_draw_extent(graphics.getExtent()) return graphics @@ -1216,7 +1211,7 @@ class MapAxes(Axes): :param cmap: (*string*) Color map string. :param colors: (*list*) If None (default), the colormap specified by cmap will be used. If a string, like ‘r’ or ‘red’, all levels will be plotted in this color. If a tuple of matplotlib - color args (string, float, rgb, etc), different levels will be plotted in different colors in + color args (string, float, rgb, etc.), different levels will be plotted in different colors in the order specified. :param fill_value: (*float*) Fill_value. Default is ``-9999.0``. :param proj: (*ProjectionInfo*) Map projection of the data. Default is None. @@ -1224,7 +1219,7 @@ class MapAxes(Axes): :param zorder: (*int*) Z-order of created layer for display. :param select: (*boolean*) Set the return layer as selected layer or not. - :returns: (*VectoryLayer*) Polygon VectoryLayer created from array data. + :returns: (*graphics*) Polygon graphics created from array data. """ n = len(args) if n <= 2: @@ -1256,8 +1251,7 @@ class MapAxes(Axes): if zorder is None: zorder = 0 graphics = self.add_graphic(graphics, projection=proj, zorder=zorder) - self._axes.setDrawExtent(graphics.getExtent().clone()) - self._axes.setExtent(graphics.getExtent().clone()) + self.set_draw_extent(graphics.getExtent()) return graphics @@ -1354,9 +1348,8 @@ class MapAxes(Axes): visible = kwargs.pop('visible', True) if visible: zorder = kwargs.pop('zorder', None) - rGraphics = self.add_graphic(graphics, projection=proj, zorder=zorder) - self._axes.setDrawExtent(rGraphics.getExtent().clone()) - self._axes.setExtent(rGraphics.getExtent().clone()) + graphics = self.add_graphic(graphics, projection=proj, zorder=zorder) + self.set_draw_extent(graphics.getExtent()) return graphics @@ -1452,8 +1445,7 @@ class MapAxes(Axes): if visible: zorder = kwargs.pop('zorder', None) graphics = self.add_graphic(graphics, projection=proj, zorder=zorder) - self._axes.setDrawExtent(graphics.getExtent().clone()) - self._axes.setExtent(graphics.getExtent().clone()) + self.set_draw_extent(graphics.getExtent()) return graphics @@ -1531,8 +1523,7 @@ class MapAxes(Axes): if visible: zorder = kwargs.pop('zorder', None) graphics = self.add_graphic(graphics, projection=proj, zorder=zorder) - self._axes.setDrawExtent(graphics.getExtent().clone()) - self._axes.setExtent(graphics.getExtent().clone()) + self.set_draw_extent(graphics.getExtent()) return graphics @@ -1570,8 +1561,7 @@ class MapAxes(Axes): if visible: zorder = kwargs.pop('zorder', None) graphics = self.add_graphic(graphics, projection=proj, zorder=zorder) - self._axes.setDrawExtent(graphics.getExtent().clone()) - self._axes.setExtent(graphics.getExtent().clone()) + self.set_draw_extent(graphics.getExtent()) return graphics diff --git a/meteoinfo-map/src/main/java/org/meteoinfo/map/forms/FrmMeteoData.java b/meteoinfo-map/src/main/java/org/meteoinfo/map/forms/FrmMeteoData.java index 19ac85a3..a48728c0 100644 --- a/meteoinfo-map/src/main/java/org/meteoinfo/map/forms/FrmMeteoData.java +++ b/meteoinfo-map/src/main/java/org/meteoinfo/map/forms/FrmMeteoData.java @@ -1044,9 +1044,11 @@ public class FrmMeteoData extends javax.swing.JDialog { viewGridData(); } else if (_meteoDataInfo.isStationData()) { viewStationData(); - if (_gridData.getXArray() != null && _gridData.getYArray() != null) { - if (_gridData.getXNum() > 0 && _gridData.getYNum() > 0) { - viewGridData(); + if (_gridData != null) { + if (_gridData.getXArray() != null && _gridData.getYArray() != null) { + if (_gridData.getXNum() > 0 && _gridData.getYNum() > 0) { + viewGridData(); + } } } } diff --git a/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/math/ArrayMath.java b/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/math/ArrayMath.java index b81815b2..f30f8d41 100644 --- a/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/math/ArrayMath.java +++ b/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/math/ArrayMath.java @@ -3548,6 +3548,110 @@ public class ArrayMath { return r; } + /** + * Clip array by minimum and maximum values + * + * @param a Input array + * @param min Minimum + * @param max Maximum + * @return Output array + */ + public static Array clip(Array a, double min, double max) { + Array r = Array.factory(a.getDataType(), a.getShape()); + double v; + if (a.getIndexPrivate().isFastIterator()) { + for (int i = 0; i < a.getSize(); i++) { + v = a.getDouble(i); + if (v < min) { + v = min; + } else if (v > max) { + v = max; + } + r.setDouble(i, v); + } + } else { + IndexIterator iterA = a.getIndexIterator(); + IndexIterator iterR = r.getIndexIterator(); + while (iterA.hasNext()) { + v = iterA.getDoubleNext(); + if (v < min) { + v = min; + } else if (v > max) { + v = max; + } + iterR.setDoubleNext(v); + } + } + + return r; + } + + /** + * Clip array by minimum value + * + * @param a Input array + * @param min Minimum + * @return Output array + */ + public static Array clipMin(Array a, double min) { + Array r = Array.factory(a.getDataType(), a.getShape()); + double v; + if (a.getIndexPrivate().isFastIterator()) { + for (int i = 0; i < a.getSize(); i++) { + v = a.getDouble(i); + if (v < min) { + v = min; + } + r.setDouble(i, v); + } + } else { + IndexIterator iterA = a.getIndexIterator(); + IndexIterator iterR = r.getIndexIterator(); + while (iterA.hasNext()) { + v = iterA.getDoubleNext(); + if (v < min) { + v = min; + } + iterR.setDoubleNext(v); + } + } + + return r; + } + + /** + * Clip array by maximum value + * + * @param a Input array + * @param max Maximum + * @return Output array + */ + public static Array clipMax(Array a, double max) { + Array r = Array.factory(a.getDataType(), a.getShape()); + double v; + if (a.getIndexPrivate().isFastIterator()) { + for (int i = 0; i < a.getSize(); i++) { + v = a.getDouble(i); + if (v > max) { + v = max; + } + r.setDouble(i, v); + } + } else { + IndexIterator iterA = a.getIndexIterator(); + IndexIterator iterR = r.getIndexIterator(); + while (iterA.hasNext()) { + v = iterA.getDoubleNext(); + if (v > max) { + v = max; + } + iterR.setDoubleNext(v); + } + } + + return r; + } + /** * Return the complex conjugate, element-wise. * The complex conjugate of a complex number is obtained by changing the sign of its imaginary part.