#----------------------------------------------------- # 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 ._ndarray import NDArray import math import datetime import numbers import mipylib.miutil as miutil from java.lang import Double from java.util import ArrayList nan = Double.NaN def dimension(value, name='null', type=None): """ Create a new Dimension. :param value: (*array_like*) Dimension value. :param name: (*string*) Dimension name. :param type: (*string*) Dimension type ['X' | 'Y' | 'Z' | 'T']. """ if isinstance(value, NDArray): value = value.aslist() dtype = DimensionType.OTHER if not type is None: if type.upper() == 'X': dtype = DimensionType.X elif type.upper() == 'Y': dtype = DimensionType.Y elif type.upper() == 'Z': dtype = DimensionType.Z elif type.upper() == 'T': dtype = DimensionType.T dim = Dimension(dtype) dim.setDimValues(value) dim.setShortName(name) return dim # 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, squeeze) 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) 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 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)