2025-12-05 11:42:54 +08:00

1138 lines
38 KiB
Python

#-----------------------------------------------------
# Author: Yaqiang Wang
# Date: 2014-12-27
# Purpose: MeteoInfo dimarray module
# Note: Jython
#-----------------------------------------------------
from org.meteoinfo.projection import KnownCoordinateSystems, Reproject, ProjectionInfo
from org.meteoinfo.data import GridData, GridArray
from org.meteoinfo.ndarray.math import ArrayMath, ArrayUtil
from org.meteoinfo.geometry.geoprocess import GeometryUtil
from org.meteoinfo.common import ResampleMethods
from org.meteoinfo.common import PointD
from org.meteoinfo.ndarray import Array, Range, MAMath, DataType
from org.meteoinfo.data.dimarray import Dimension, DimensionType
from java.lang import Double
from java.util import ArrayList
import math
import datetime
import numbers
import mipylib.miutil as miutil
import mipylib.numeric as np
from mipylib.numeric.core._ndarray import NDArray
from .accessor_dt import DateTimeAccessor
nan = Double.NaN
# Dimension array
class DimArray(NDArray):
def __init__(self, array, dims=None, proj=ProjectionInfo.LONG_LAT):
if isinstance(array, NDArray):
array = array._array
super(DimArray, self).__init__(array)
self.dims = None
if dims is None or len(dims) == 0:
dims = []
for i in range(self.ndim):
dim = Dimension()
dim.setDimValues(range(self.shape[i]))
dims.append(dim)
for dim in dims:
self.adddim(dim)
self.proj = proj
def __getitem__(self, indices):
if not isinstance(indices, tuple):
indices = [indices]
#deal with Ellipsis
if Ellipsis in indices:
indices1 = []
n = self.ndim - len(indices) + 1
for ii in indices:
if ii is Ellipsis:
for _ in range(n):
indices1.append(slice(None))
else:
indices1.append(ii)
indices = indices1
#for all int indices
if len(indices) == self.ndim:
allint = True
aindex = self._array.getIndex()
i = 0
for ii in indices:
if isinstance(ii, int):
if ii < 0:
ii = self.shape[i] + ii
aindex.setDim(i, ii)
else:
allint = False
break
i += 1
if allint:
return self._array.getObject(aindex)
#for None index - np.newaxis
newaxis = []
if len(indices) > self.ndim:
nindices = []
i = 0
for ii in indices:
if ii is None:
newaxis.append(i)
else:
nindices.append(ii)
i += 1
indices = nindices
#add slice(None) to end
if len(indices) < self.ndim:
if isinstance(indices, tuple):
indices = list(indices)
for i in range(self.ndim - len(indices)):
indices.append(slice(None))
if len(indices) != self.ndim:
print('indices must be ' + str(self.ndim) + ' dimensions!')
raise IndexError()
if not self.proj is None and not self.proj.isLonLat():
xlim = None
ylim = None
xidx = -1
yidx = -1
for i in range(0, self.ndim):
dim = self.dims[i]
if dim.getDimType() == DimensionType.X:
k = indices[i]
#if isinstance(k, (tuple, list)):
if isinstance(k, basestring):
xlims = k.split(':')
xlim = [float(xlims[0]), float(xlims[1])]
xidx = i
elif dim.getDimType() == DimensionType.Y:
k = indices[i]
#if isinstance(k, (tuple, list)):
if isinstance(k, basestring):
ylims = k.split(':')
ylim = [float(ylims[0]), float(ylims[1])]
yidx = i
if not xlim is None and not ylim is None:
fromproj = KnownCoordinateSystems.geographic.world.WGS1984
inpt = PointD(xlim[0], ylim[0])
outpt1 = Reproject.reprojectPoint(inpt, fromproj, self.proj)
inpt = PointD(xlim[1], ylim[1])
outpt2 = Reproject.reprojectPoint(inpt, fromproj, self.proj)
xlim = [outpt1.X, outpt2.X]
ylim = [outpt1.Y, outpt2.Y]
indices1 = []
for i in range(0, self.ndim):
if i == xidx:
#indices1.append(xlim
indices1.append(str(xlim[0]) + ':' + str(xlim[1]))
elif i == yidx:
#indices1.append(ylim)
indices1.append(str(ylim[0]) + ':' + str(ylim[1]))
else:
indices1.append(indices[i])
indices = indices1
ndims = []
ranges = []
flips = []
iszerodim = True
onlyrange = True
alllist = True
isempty = False
nshape = []
squeeze = False
for i in range(0, self.ndim):
isrange = True
k = indices[i]
if isinstance(k, int):
if k < 0:
k = self.dims[i].getLength() + k
sidx = k
eidx = k
step = 1
alllist = False
squeeze = True
elif isinstance(k, slice):
if isinstance(k.start, basestring):
sv = float(k.start)
sidx = self.dims[i].getValueIndex(sv)
elif isinstance(k.start, datetime.datetime):
sv = miutil.date2num(k.start)
sidx = self.dims[i].getValueIndex(sv)
else:
sidx = 0 if k.start is None else k.start
if sidx < 0:
sidx = self.dims[i].getLength() + sidx
if isinstance(k.stop, basestring):
ev = float(k.stop)
eidx = self.dims[i].getValueIndex(ev)
elif isinstance(k.stop, datetime.datetime):
ev = miutil.date2num(k.stop)
eidx = self.dims[i].getValueIndex(ev)
else:
eidx = self.dims[i].getLength() if k.stop is None else k.stop
if eidx < 0:
eidx = self.dims[i].getLength() + eidx
eidx -= 1
if isinstance(k.step, basestring):
nv = float(k.step) + self.dims[i].getDimValue()[0]
nidx = self.dims[i].getValueIndex(nv)
step = nidx - sidx
elif isinstance(k.step, datetime.timedelta):
nv = miutil.date2num(k.start + k.step)
nidx = self.dims[i].getValueIndex(nv)
step = nidx - sidx
else:
step = 1 if k.step is None else k.step
alllist = False
elif isinstance(k, (list, tuple, NDArray)):
onlyrange = False
isrange = False
if not isinstance(k[0], datetime.datetime):
if isinstance(k, (list, tuple)):
k = NDArray(k)
# if isinstance(k, NDArray):
# k = k.aslist()
# ranges.append(k)
else:
tlist = []
for tt in k:
sv = miutil.date2num(tt)
idx = self.dims[i].getValueIndex(sv)
tlist.append(idx)
# ranges.append(tlist)
k = NDArray(tlist)
ranges.append(k.asarray())
elif isinstance(k, basestring):
dim = self.dims[i]
kvalues = k.split(':')
sidx = dim.getValueIndex(float(kvalues[0]))
if len(kvalues) == 1:
eidx = sidx
step = 1
squeeze = True
else:
eidx = dim.getValueIndex(float(kvalues[1]))
if len(kvalues) == 2:
step = 1
else:
step = int(float(kvalues[2]) / dim.getDeltaValue())
if sidx > eidx and step > 0:
step = -step
alllist = False
else:
print(k)
raise IndexError()
if isrange:
if sidx >= self.shape[i] or eidx >= self.shape[i]:
raise IndexError()
n = abs(eidx - sidx) + 1
if n > 1:
dim = self.dims[i]
ndims.append(dim.extract(sidx, eidx, step))
if step < 0:
#step = abs(step)
flips.append(i)
if sidx > eidx:
iidx = eidx
eidx = sidx
sidx = iidx
rr = Range(sidx, eidx, abs(step))
ranges.append(rr)
nshape.append(eidx - sidx + 1 if eidx - sidx >= 0 else 0)
else:
if len(k) > 1:
dim = self.dims[i]
if isinstance(k, NDArray):
k = k.asarray()
ndims.append(dim.extract(k))
if isempty:
r = ArrayUtil.zeros(nshape, 'int')
return NDArray(r)
if onlyrange:
r = ArrayMath.section(self._array, ranges, False)
else:
if alllist:
r = ArrayMath.takeValues(self._array, ranges)
return NDArray(r)
else:
r = ArrayMath.take(self._array, ranges)
if newaxis:
for i in flips:
r = r.flip(i)
rr = Array.factory(r.getDataType(), r.getShape())
MAMath.copy(rr, r)
rr = NDArray(rr)
newshape = list(rr.shape)
for i in newaxis:
newshape.insert(i, 1)
rr = rr.reshape(newshape)
if onlyrange:
rr.base = self.get_base()
return rr
for i in flips:
r = r.flip(i)
r = r.reduce()
data = DimArray(r, ndims, self.proj)
if onlyrange:
data.base = self.get_base()
return data
@property
def array(self):
"""
Get NDArray object
:return: NDArray object
"""
return NDArray(self._array)
@array.setter
def array(self, value):
"""
Set NDArray object
:param value: NDArray object
"""
self._array = value._array
@property
def dt(self):
a = self.array
if a.dtype == np.dtype.datetime:
return DateTimeAccessor(a)
else:
return None
def member_names(self):
"""
Get member names. Only valid for Structure data type.
:returns: (*list*) Member names
"""
if self._array.getDataType() != DataType.STRUCTURE:
print('This method is only valid for structure array!')
return None
ms = self._array.getStructureMemberNames()
return list(ms)
def member_array(self, member, indices=None):
"""
Extract member array. Only valid for Structure data type.
:param member: (*string*) Member name.
:param indices: (*slice*) Indices.
:returns: (*array*) Extracted member array.
"""
if self._array.getDataType() != DataType.STRUCTURE:
print('This method is only valid for structure array!')
return None
a = self._array.getArrayObject()
m = a.findMember(member)
if m is None:
raise KeyError('The member %s not exists!' % member)
a = a.extractMemberArray(m)
r = DimArray(a, self.dims, self.proj)
if not indices is None:
r = r.__getitem__(indices)
return r
def value(self, indices):
#print type(indices)
if not isinstance(indices, tuple):
inds = []
inds.append(indices)
indices = inds
if len(indices) != self.ndim:
print('indices must be ' + str(self.ndim) + ' dimensions!')
return None
#origin = []
#size = []
#stride = []
dims = []
ranges = []
flips = []
for i in range(0, self.ndim):
k = indices[i]
if isinstance(indices[i], int):
sidx = k
eidx = k
step = 1
elif isinstance(k, slice):
sidx = 0 if k.start is None else k.start
eidx = self.dims[i].getLength()-1 if k.stop is None else k.stop
step = 1 if k.step is None else k.step
elif isinstance(k, tuple) or isinstance(k, list):
dim = self.dims[i]
sidx = dim.getValueIndex(k[0])
if len(k) == 1:
eidx = sidx
step = 1
else:
eidx = dim.getValueIndex(k[1])
if len(k) == 2:
step = 1
else:
step = int(k[2] / dim.getDeltaValue)
else:
print(k)
return None
if step < 0:
step = abs(step)
flips.append(i)
rr = Range(sidx, eidx, step)
ranges.append(rr)
#origin.append(sidx)
n = eidx - sidx + 1
#size.append(n)
#stride.append(step)
if n > 1:
dim = self.dims[i]
dims.append(dim.extract(sidx, eidx, step))
#r = ArrayMath.section(self._array, origin, size, stride)
r = ArrayMath.section(self._array, ranges)
for i in flips:
r = r.flip(i)
rr = Array.factory(r.getDataType(), r.getShape());
MAMath.copy(rr, r);
array = NDArray(rr)
data = DimArray(array, dims, self.proj)
return data
def _ufunc_finalize(self, obj, axis=None):
"""
Return a new array after universal function finalized.
:param obj: The object output from the universal function.
:param axis: (*int*) The axis for ufunc compute along. Default is `None`, means not consider.
:return: New array object.
"""
if isinstance(obj, (Array, NDArray)):
if axis is None:
return DimArray(obj, self.dims, self.proj)
else:
dims = []
for i in range(0, self.ndim):
if i != axis:
dims.append(self.dims[i])
return DimArray(obj, dims, self.proj)
else:
return obj
def array_wrap(self, arr, axis=None):
"""
Return a new array wrapped as self class object.
:param arr: The array to be wrapped.
:param axis: (*int*) The axis for ufunc compute along. Default is `None`, means not consider.
:return: New array object.
"""
if isinstance(arr, (Array, NDArray)):
if axis is None:
return DimArray(arr, self.dims, self.proj)
else:
dims = []
for i in range(0, self.ndim):
if isinstance(axis, (list, tuple)):
if not i in axis:
dims.append(self.dims[i])
else:
if i != axis:
dims.append(self.dims[i])
return DimArray(arr, dims, self.proj)
else:
return arr
# get dimension length
def dimlen(self, idx=0):
"""
Get dimension length.
:param idx: (*int*) Dimension index.
:returns: (*int*) Dimension length.
"""
return self.dims[idx].getLength()
def dimvalue(self, idx=0, convert=False):
"""
Get dimension values.
:param idx: (*int*) Dimension index.
:param convert: (*boolean*) If convert to real values (i.e. datetime). Default
is ``False``.
:returns: (*array_like*) Dimension values
"""
dim = self.dims[idx]
if convert:
if dim.getDimType() == DimensionType.T:
return miutil.nums2dates(dim.getDimValue())
else:
return NDArray(self.dims[idx].getDimValue())
else:
return NDArray(self.dims[idx].getDimValue())
def setdimvalue(self, idx, dimvalue):
"""
Set dimension value.
:param idx: (*int*) Dimension index.
:param dimvalue: (*array_like*) Dimension value.
"""
if isinstance(dimvalue, NDArray):
self.dims[idx].setDimValue(dimvalue._array)
else:
self.dims[idx].setDimValues(dimvalue)
def setdimtype(self, idx, dimtype):
"""
Set dimension type.
:param idx: (*int*) Dimension index.
:param dimtype: (*string*) Dimension type. [X | Y | Z | T].
"""
dtype = DimensionType.Other
if dimtype.upper() == 'X':
dtype = DimensionType.X
elif dimtype.upper() == 'Y':
dtype = DimensionType.Y
elif dimtype.upper() == 'Z':
dtype = DimensionType.Z
elif dimtype.upper() == 'T':
dtype = DimensionType.T
self.dims[idx].setDimType(dtype)
def adddim(self, dimvalue, dimtype=None, index=None):
"""
Add a dimension.
:param dimvalue: (*array_like*) Dimension value.
:param dimtype: (*string*) Dimension type.
:param index: (*int*) Index to be inserted.
"""
if isinstance(dimvalue, Dimension):
dim = dimvalue
else:
if isinstance(dimvalue, (NDArray, DimArray)):
dimvalue = dimvalue.aslist()
dtype = DimensionType.OTHER
if not dimtype is None:
if dimtype.upper() == 'X':
dtype = DimensionType.X
elif dimtype.upper() == 'Y':
dtype = DimensionType.Y
elif dimtype.upper() == 'Z':
dtype = DimensionType.Z
elif dimtype.upper() == 'T':
dtype = DimensionType.T
dim = Dimension(dtype)
dim.setDimValues(dimvalue)
if self.dims is None:
self.dims = [dim]
else:
if index is None:
self.dims.append(dim)
else:
self.dims.insert(index, dim)
self.ndim = len(self.dims)
def addtdim(self, t):
"""
Add a time dimension as first dimension.
:param t: (*array_like*) datetime array.
"""
if self.tdim() is None:
dim = Dimension(DimensionType.T)
t = miutil.date2num(t)
dim.setDimValues([t])
if self.dims is None:
self.dims = [dim]
else:
self.dims.insert(0, dim)
self.ndim = len(self.dims)
ss = list(self.shape)
ss.insert(0, 1)
ss = tuple(ss)
self._array = self._array.reshape(ss)
#self.shape = self._array.shape
def xdim(self):
"""
Get x dimension.
"""
for dim in self.dims:
if dim.getDimType() == DimensionType.X:
return dim
return None
def ydim(self):
"""
Get y dimension.
"""
for dim in self.dims:
if dim.getDimType() == DimensionType.Y:
return dim
return None
def zdim(self):
"""
Get z dimension.
"""
for dim in self.dims:
if dim.getDimType() == DimensionType.Z:
return dim
return None
def tdim(self):
"""
Get time dimension.
"""
if self.dims is None:
return None
for dim in self.dims:
if dim.getDimType() == DimensionType.T:
return dim
return None
def islondim(self, idx=0):
"""
Check a dimension is a longitude dimension or not.
:param idx: (*int*) Dimension index.
"""
dim = self.dims[idx]
if dim.getDimType() == DimensionType.X and self.proj.isLonLat():
return True
else:
return False
def islatdim(self, idx=0):
"""
Check a dimension is a latitude dimension or not.
:param idx: (*int*) Dimension index.
"""
dim = self.dims[idx]
if dim.getDimType() == DimensionType.Y and self.proj.isLonLat():
return True
else:
return False
def islonlatdim(self, idx=0):
"""
Check a dimension is a longitude or latitude dimension or not.
:param idx: (*int*) Dimension index.
"""
return self.islondim(idx) or self.islatdim(idx)
def istimedim(self, idx=0):
"""
Check a dimension is a time dimension or not.
:param idx: (*int*) Dimension index.
"""
dim = self.dims[idx]
if dim.getDimType() == DimensionType.T:
return True
else:
return False
def asgriddata(self, xdata=None, ydata=None, fill_value=-9999.0):
if xdata is None or ydata is None:
xdata = self.dimvalue(1)
ydata = self.dimvalue(0)
data = self.array
if xdata[1] < xdata[0]:
xdata = xdata[::-1]
data = self.array[:,::-1]
if ydata[1] < ydata[0]:
ydata = ydata[::-1]
data = self.array[::-1,:]
gdata = GridData(data._array, xdata._array, ydata._array, fill_value, self.proj)
return gdata
def asgridarray(self, xdata=None, ydata=None, fill_value=-9999.0):
if xdata is None or ydata is None:
xdata = self.dimvalue(1)
ydata = self.dimvalue(0)
data = self.array
if xdata[1] < xdata[0]:
xdata = xdata[::-1]
data = self.array[:,::-1]
if ydata[1] < ydata[0]:
ydata = ydata[::-1]
data = self.array[::-1,:]
gdata = GridArray(data._array, xdata._array, ydata._array, fill_value, self.proj)
return gdata
def squeeze(self):
"""
Remove single-dimensional entries from the shape of an array.
:returns: (*array_like*) The self array, but with all or a subset of the dimensions of length 1
removed.
"""
r = self._array.reduce()
dims = []
for dim in a.dims:
if dim.getLength() > 1:
dims.append(dim)
return DimArray(r, dims, self.proj)
def mean(self, axis=None, keepdims=False):
"""
Compute tha arithmetic mean along the specified axis.
:param axis: (*int*) Axis along which the value is computed.
The default is to compute the value of the flattened array.
:param keepdims: (*bool*) If this is set to True, the axes which are reduced are
left in the result as dimensions with size one. Default if `False`.
returns: (*array_like*) Mean result
"""
r = super(DimArray, self).mean(axis, keepdims)
if isinstance(r, numbers.Number):
return r
else:
dims = []
for i in range(0, self.ndim):
if isinstance(axis, (list, tuple)):
if not i in axis:
dims.append(self.dims[i])
else:
if i != axis:
dims.append(self.dims[i])
return DimArray(r, dims, self.proj)
def rot90(self, k=1):
"""
Rotate an array by 90 degrees in the counter-clockwise direction. The first two dimensions
are rotated if the array has more than 2 dimensions.
:param k: (*int*) Number of times the array is rotated by 90 degrees
:returns: (*array_like*) Rotated array.
"""
r = ArrayMath.rot90(self._array, k)
dims = []
if math.fabs(k) == 1 or math.fabs(k) == 3:
dims.append(self.dims[1])
dims.append(self.dims[0])
for i in range(2, len(self.dims)):
dims.append(self.dims[i])
else:
for i in range(0, len(self.dims)):
dims.append(self.dims[i])
return DimArray(r, dims, self.proj)
def maskout(self, mask):
"""
Maskout data by polygons - the elements outside polygons will be set as NaN.
:param mask: (*list*) Polygon list as mask borders.
:returns: (*DimArray*) Maskouted data.
"""
if isinstance(mask, NDArray):
r = ArrayMath.maskout(self.asarray(), mask.asarray())
return DimArray(NDArray(r), self.dims, self.proj)
else:
x = self.dims[1].getDimValue()
y = self.dims[0].getDimValue()
xy = ArrayUtil.meshgrid([x,y])
x = xy[0]
y = xy[1]
if not isinstance(mask, (list, ArrayList)):
mask = [mask]
r = GeometryUtil.maskout(self.asarray(), x, y, mask)
r = DimArray(NDArray(r), self.dims, self.proj)
return r
def maskin(self, mask):
"""
Maskin data by polygons - the elements inside polygons will be set as NaN.
:param mask: (*list*) Polygon list as mask borders.
:returns: (*DimArray*) Maskined data.
"""
if isinstance(mask, NDArray):
r = ArrayMath.maskin(self.asarray(), mask.asarray())
return DimArray(r, self.dims, self.proj)
else:
x = self.dimvalue(1)
y = self.dimvalue(0)
if not isinstance(mask, (list, ArrayList)):
mask = [mask]
r = GeometryUtil.maskin(self._array, x._array, y._array, mask)
r = DimArray(r, self.dims, self.proj)
return r
def transpose(self, axes=None):
"""
Permute the dimensions of an array.
:param axes: (*list of int*) By default, reverse the dimensions, otherwise permute the axes according to the
values given.
:returns: Permuted array.
"""
if axes is None:
axes = [self.ndim-i-1 for i in range(self.ndim)]
r = super(DimArray, self).transpose(axes)
if self.ndim == 1:
dims = self.dims
else:
dims = []
for ax in axes:
dims.append(self.dims[ax])
return DimArray(r, dims, self.proj)
T = property(transpose)
def swapaxes(self, axis1, axis2):
"""
Interchange two axes of an array.
:param axis1: (*int*) First axis.
:param axis2: (*int*) Second axis.
:returns: Axes swapped array.
"""
if self.ndim == 1:
return self
if axis1 < 0:
axis1 = self.ndim + axis1
if axis2 < 0:
axis2 = self.ndim + axis2
if axis1 == axis2:
return self
r = self._array.transpose(axis1, axis2)
dims = []
for i in range(self.ndim):
if i == axis1:
dims.append(self.dims[axis2])
elif i == axis2:
dims.append(self.dims[axis1])
else:
dims.append(self.dims[i])
return DimArray(r, dims, self.proj)
def inv(self):
"""
Calculate inverse matrix array.
:returns: Inverse matrix array.
"""
r = super(DimArray, self).inv()
return DimArray(r, self.dims, self.proj)
I = property(inv)
def lonflip(self):
"""
Reorder global array from 0 - 360 longitude to -180 - 180 longitude or vice versa.
:returns: Reordered array.
"""
lon = self.dimvalue(self.ndim - 1)
if lon.max() > 180:
return self.lonpivot(180)
else:
return self.lonpivot(0)
def lonpivot(self, pivot):
"""
Pivots an array about a user-specified longitude. The rightmost dimension must be the
longitude dimension, which must be global and without a cyclic point.
:param pivot: (*float*) The longitude value around which to pivot.
:returns: Result array after longitude pivot.
"""
lon = self.dimvalue(self.ndim - 1)
minlon = lon.min()
maxlon = lon.max()
dlon = lon[1] - lon[0]
if pivot < minlon:
pivot += 360
elif pivot > maxlon:
pivot -= 360
keys1 = []
keys2 = []
for i in range(self.ndim - 1):
keys1.append(slice(None,None,None))
keys2.append(slice(None,None,None))
keys1.append('%f:%f' % (pivot, maxlon))
keys2.append('%f:%f' % (minlon, pivot - dlon))
r1 = self.__getitem__(tuple(keys1))
r2 = self.__getitem__(tuple(keys2))
lon1 = r1.dimvalue(r1.ndim - 1)
lon2 = r2.dimvalue(r2.ndim - 1)
if maxlon > 180:
lon1 = lon1 - 360
else:
lon2 = lon2 + 360
r = r1.join(r2, self.ndim - 1)
lon1 = lon1.aslist()
lon1.extend(lon2.aslist())
r.setdimvalue(self.ndim - 1, lon1)
return r
def month_to_season(self, season):
"""
Computes a user-specified three-month seasonal mean
(DJF, JFM, FMA, MAM, AMJ, MJJ, JJA, JAS, ASO, SON, OND, NDJ).
The first average (DJF=JF) and the last average (NDJ=ND) are actually
two-month averages.
The time (leftmost) dimension must be divisible by 12. The data are assumed
to be monthly mean data and the first record is assumed to be January.
:param season: (*string*) A string representing the season to
calculate: e.g., "JFM", "JJA".
:returns: Season averaged data array.
"""
nmonth = self.dimlen(0)
nyear = nmonth / 12
seasons = ['DJF','JFM','FMA','MAM','AMJ','MJJ','JJA','JAS','ASO','SON','OND','NDJ']
season = season.upper()
if not season in seasons:
print('Season string is not valid: "' + season + '"!')
raise KeyError()
idx = seasons.index(season) - 1
keys = []
keys.append(slice(0,nyear,1))
for i in range(1, self.ndim):
keys.append(slice(None,None,None))
r = self.__getitem__(tuple(keys))
si = idx
for i in range(nyear):
ei = si + 3
if si < 0:
si = 0
if ei > nmonth:
ei = nmonth
keys[0] = slice(si,ei,1)
sdata = self.__getitem__(tuple(keys))
sdata = ArrayMath.mean(sdata.asarray(), 0)
keys[0] = i
r.__setitem__(tuple(keys), sdata)
si += 12
return r
def interpn(self, xi):
"""
Multidimensional interpolation on regular grids.
:param xi: (*list*) The coordinates to sample the gridded data at.
:returns: (*float*) Interpolated value at input coordinates.
"""
points = []
for i in range(self.ndim):
points.append(ArrayUtil.array(self.dims[i].getDimValue()))
if isinstance(xi, (list, tuple)):
if isinstance(xi[0], NDArray):
nxi = []
for x in xi:
nxi.append(x._array)
else:
nxi = []
for x in xi:
if isinstance(x, datetime.datetime):
x = miutil.date2num(x)
nxi.append(x)
nxi = NDArray(nxi)._array
else:
nxi = nxi._array
r = ArrayUtil.interpn(points, self._array, nxi)
if isinstance(r, Array):
return NDArray(r)
else:
return r
def project(self, x=None, y=None, toproj=None, method='bilinear'):
"""
Project array
:param x: To x coordinates.
:param y: To y coordinates.
:param toproj: To projection.
:param method: Interpolation method: ``bilinear`` or ``nearest``. Default is ``bilinear``.
:returns: (*NDArray*) Projected array
"""
yy = self.dims[self.ndim - 2].getDimValue()
xx = self.dims[self.ndim - 1].getDimValue()
if toproj is None:
toproj = self.proj
if x is None or y is None:
pr = Reproject.reproject(self._array, xx, yy, self.proj, toproj)
r = pr[0]
x = pr[1]
y = pr[2]
dims = self.dims
ydim = Dimension(DimensionType.Y)
ydim.setDimValues(NDArray(y).aslist())
dims[-2] = ydim
xdim = Dimension(DimensionType.X)
xdim.setDimValues(NDArray(x).aslist())
dims[-1] = xdim
rr = DimArray(NDArray(r), dims, toproj)
return rr
if method == 'bilinear':
method = ResampleMethods.Bilinear
else:
method = ResampleMethods.NearestNeighbor
if isinstance(x, (list, tuple)):
x = NDArray(x)
if isinstance(y, (list, tuple)):
y = NDArray(y)
if x.ndim == 1:
x, y = ArrayUtil.meshgrid(x.asarray(), y.asarray())
else:
x = x._array
y = y._array
r = Reproject.reproject(self._array, xx, yy, x, y, self.proj, toproj, method)
return NDArray(r)
def join(self, b, dimidx):
r = ArrayMath.join(self._array, b._array, dimidx)
dima = self.dimvalue(dimidx)
dimb = b.dimvalue(dimidx)
dimr = []
if dima[0] < dimb[0]:
for i in range(0, len(dima)):
dimr.append(dima[i])
for i in range(0, len(dimb)):
dimr.append(dimb[i])
else:
for i in range(0, len(dimb)):
dimr.append(dimb[i])
for i in range(0, len(dima)):
dimr.append(dima[i])
rdims = []
for i in range(0, len(self.dims)):
if i == dimidx:
ndim = Dimension()
ndim.setDimValues(dimr)
rdims.append(ndim)
else:
rdims.append(self.dims[i])
return DimArray(NDArray(r), rdims, self.proj)
def savegrid(self, fname, format='surfer', **kwargs):
"""
Save the array data to an ASCII or binary file. The array must be 2 dimension.
:param fname: (*string*) File name.
:param format: (*string*) File format [surfer | bil | esri_ascii | micaps4].
:param description: (*string*) Data description - only used for ``micaps4`` file.
:param date: (*datetime*) Data datetime - only used for ``micaps4`` file.
:param hours: (*int*) Data forcasting hours - only used for ``micaps4`` file.
:param level: (*float*) Data vertical level - only used for ``micaps4`` file.
:param smooth: (*int*) 1 or 0 - only used for ``micaps4`` file.
:param boldvalue: (*int*) Bold contour value - only used for ``micaps4`` file.
:param proj: (*ProjectionInfo*) Data ProjectionInfo - only used for ``micaps4`` file.
:param float_format: (*string*) Float number format, such as '%.2f'.
"""
if self.ndim != 2:
print('The array must be 2 dimensional!')
return
gdata = self.asgridarray()
float_format = kwargs.pop('float_format', None)
if format == 'surfer':
gdata.saveAsSurferASCIIFile(fname)
elif format == 'bil':
gdata.saveAsBILFile(fname)
elif format == 'esri_ascii':
gdata.saveAsESRIASCIIFile(fname)
elif format == 'micaps4':
desc = kwargs.pop('description', 'var')
date = kwargs.pop('date', datetime.datetime.now())
date = miutil.jdate(date)
hours = kwargs.pop('hours', 0)
level = kwargs.pop('level', 0)
smooth = kwargs.pop('smooth', 1)
boldvalue =kwargs.pop('boldvalue', 0)
proj = kwargs.pop('proj', self.proj)
if proj is None:
gdata.saveAsMICAPS4File(fname, desc, date, hours, level, smooth, boldvalue, float_format)
else:
if proj.isLonLat():
gdata.saveAsMICAPS4File(fname, desc, date, hours, level, smooth, boldvalue, float_format)
else:
gdata.saveAsMICAPS4File(fname, desc, date, hours, level, smooth, boldvalue, float_format, proj)
def dim_array(a, dims=None):
"""
Create a dimension array (DimArray).
:param a: (*array_like*) Array (NDArray) or data list.
:param dims: (*list*) List of dimensions.
:returns: (*DimArray*) Dimension array.
"""
if not isinstance(a, NDArray):
a = array(a)
if dims is None:
dims = []
for i in range(a.ndim):
dim = Dimension()
dim.setDimValues(range(a.shape[i]))
dims.append(dim)
return DimArray(a, dims)