# Automatically adapted for numpy.oldnumeric Aug 01, 2007 by
# Further modified to be pure new numpy June 24th 2008
"""CDMS HorizontalGrid objects"""
from __future__ import print_function
import numpy
import cdms2
import os
import os.path
# import PropertiedClasses
from .error import CDMSError
from .grid import AbstractGrid, LongitudeType, LatitudeType, CoordTypeToLoc
from .axis import TransientVirtualAxis
from .axis import getAutoBounds, allclose
from cdms2 import bindex
from cdms2 import _bindex
from functools import reduce
import copy
MethodNotImplemented = "Method not yet implemented"
[docs]def _flatten(boundsar):
boundsshape = boundsar.shape
if len(boundsshape) > 2:
newshape = (reduce((lambda x, y: x * y),
boundsshape[:-1], 1), boundsshape[-1])
boundsar.shape = newshape
return boundsar
class AbstractHorizontalGrid(AbstractGrid):
""" Create an horizontal grid
Parameters
----------
latAxis
lonAxis
id - Default None
maskvar - Default None
tempmask - Default None
node - Default None
"""
[docs] def __init__(self, latAxis, lonAxis, id=None,
maskvar=None, tempmask=None, node=None):
"""Create a horizontal grid.
"""
AbstractGrid.__init__(self, node)
if id is None:
self.id = "<None>"
else:
self.id = id
self._lataxis_ = latAxis
self._lonaxis_ = lonAxis
self._maskVar_ = maskvar
self._tempMask_ = tempmask
# Generate default bounds
[docs] def genBounds(self):
"""
Not documented
"""
raise CDMSError(MethodNotImplemented)
# Get the n-th axis. naxis is 0 or 1.
[docs] def getAxis(self, naxis):
"""
Not documented
"""
raise CDMSError(MethodNotImplemented)
[docs] def getBounds(self):
"""Get the grid cell boundaries, as a tuple (latitudeBounds, longitudeBounds)
"""
latbnds, lonbnds = (self._lataxis_.getExplicitBounds(),
self._lonaxis_.getExplicitBounds())
if (latbnds is None or lonbnds is None) and getAutoBounds() in [1, 2]:
nlatbnds, nlonbnds = self.genBounds()
if latbnds is None:
latbnds = nlatbnds
if lonbnds is None:
lonbnds = nlonbnds
return (latbnds, lonbnds)
[docs] def getLatitude(self):
"""Get the latitude coordinates."""
return self._lataxis_
[docs] def getLongitude(self):
"""Get the longitude coordinates."""
return self._lonaxis_
[docs] def getMask(self):
"""Get the mask array, if any, otherwise None is returned."""
if self._maskVar_ is not None:
return self._maskVar_
else:
return self._tempMask_
[docs] def getMesh(self):
"""Get the mesh array used by the meshfill plot."""
raise CDMSError(MethodNotImplemented)
[docs] def getWeightsArray(self):
"""
Get Weights Array
Returns
-------
normalized area weights, as an array of the same
shape as the grid.
"""
raise CDMSError(MethodNotImplemented)
[docs] def listall(self, all=None):
"""
Not documented
"""
result = []
result.append('Grid has Python id %s.' % hex(id(self)))
return result
[docs] def setMask(self, mask, permanent=0):
"""
Not documented
"""
self._maskVar_ = mask
[docs] def subGridRegion(self, latRegion, lonRegion):
"""
Not documented
"""
raise CDMSError(MethodNotImplemented)
[docs] def hasCoordType(self, coordType):
"""
Not documented
"""
return ((coordType == LatitudeType) or (coordType == LongitudeType))
[docs] def checkConvex(self):
"""Check that each cell of the grid is convex in lon-lat space, with nodes defined counter-clockwise.
Returns
-------
a 1D numpy array of cells that fail the cross-product test.
"""
from numpy import zeros, where, less, logical_or, compress
latb, lonb = self.getBounds()
saveshape = lonb.shape
lonb = _flatten(lonb)
latb = _flatten(latb)
ncell, nnode = lonb.shape
badmask = zeros((ncell,))
for n0 in range(nnode):
n1 = (n0 + 1) % nnode
n2 = (n1 + 1) % nnode
vec0lon = lonb[:, n1] - lonb[:, n0]
vec0lat = latb[:, n1] - latb[:, n0]
vec1lon = lonb[:, n2] - lonb[:, n1]
vec1lat = latb[:, n2] - latb[:, n1]
cross = vec0lon * vec1lat - vec0lat * vec1lon
mask = where(less(cross, 0.0), 1, 0)
badmask = logical_or(mask, badmask)
badcells = compress(badmask, list(range(len(badmask))))
lonb.shape = saveshape
latb.shape = saveshape
return badcells
[docs] def fixCutCells(self, nonConvexCells, threshold=270.0):
"""
For any mapping from a spherical to a planar surface, there is a linear cut.
Grid cells that span the cut may appear to be nonconvex, which causes problems
with meshfill graphics. This routine attempts to 'repair' the cut cell boundaries
so that meshfill recognizes they are convex.
Parameters
----------
nonConvexCells : 1D numpy array of indices of nonconvex cells, as returned from checkConvex.
threshold : positive floating-point value in degrees.
If the difference in longitude values of consecutive boundaries nodes
exceeds the threshold, the cell is considered a cut cell. On return,
the grid boundaries are modified.
Returns
-------
value is a 1D array of indices of cells that cannot be repaired.
"""
from numpy import take, array
latb, lonb = self.getBounds()
saveshape = lonb.shape
lonb = _flatten(lonb)
latb = _flatten(latb)
ncell, nnode = lonb.shape
lonb2 = take(lonb, nonConvexCells, axis=0)
latb2 = take(latb, nonConvexCells, axis=0)
newbadcells = []
for k in range(len(nonConvexCells)):
savelons = lonb2[k]
# Loop twice
for node in range(2 * nnode):
n0 = node % nnode
n1 = (n0 + 1) % nnode
vec0lon = lonb2[k, n1] - lonb2[k, n0]
if vec0lon > threshold:
lonb2[k, n1] -= 360.0
elif vec0lon < -threshold:
lonb2[k, n1] += 360.0
# If the cross-product test still fails, restore
# the original values and add to the nonConvexCells list
for n0 in range(nnode):
n1 = (n0 + 1) % nnode
n2 = (n1 + 1) % nnode
vec0lon = lonb2[k, n1] - lonb2[k, n0]
vec0lat = latb2[k, n1] - latb2[k, n0]
vec1lon = lonb2[k, n2] - lonb2[k, n1]
vec1lat = latb2[k, n2] - latb2[k, n1]
cross = vec0lon * vec1lat - vec0lat * vec1lon
if cross < 0:
lonb2[k] = savelons
newbadcells.append(nonConvexCells[k])
break
# Scatter the repaired cell bounds back to the original bounds
# and reset the grid bounds.
for k in range(len(nonConvexCells)):
lonb[nonConvexCells[k]] = lonb2[k]
lonb.shape = saveshape
self.getLongitude().setBounds(lonb)
return array(newbadcells)
class AbstractCurveGrid(AbstractHorizontalGrid):
"""Create a curvilinear grid.
"""
[docs] def __init__(self, latAxis, lonAxis, id=None,
maskvar=None, tempmask=None, node=None):
"""Create a curvilinear grid.
"""
if latAxis.shape != lonAxis.shape:
raise CDMSError(
'Latitude and longitude axes must have the same shape.')
AbstractHorizontalGrid.__init__(
self, latAxis, lonAxis, id, maskvar, tempmask, node)
self._index_ = None
[docs] def clone(self, copyData=1):
"""
Not documented
"""
newlat = self._lataxis_.clone(copyData)
newlon = self._lonaxis_.clone(copyData)
return TransientCurveGrid(newlat, newlon, id=self.id)
def __repr__(self):
return "<CurveGrid, id: %s, shape: %s>" % (self.id, repr(self.shape))
__str__ = __repr__
[docs] def getMesh(self, transpose=None):
"""Generate a mesh array for the meshfill graphics method.
If transpose is defined to a tuple, say (1,0), first transpose
latbounds and lonbounds according to the tuple, (1,0,2) in this case.
"""
if self._mesh_ is None:
LAT = 0
LON = 1
latbounds, lonbounds = self.getBounds()
# following work aronud a numpy.ma bug
# latbounds=latbounds.filled()
# lonbounds=lonbounds.filled()
if latbounds is None or lonbounds is None:
raise CDMSError(
'No boundary data is available for grid %s' %
self.id)
if (transpose is not None) and (transpose[1] == 0):
latbounds = numpy.transpose(latbounds, (1, 0, 2))
lonbounds = numpy.transpose(lonbounds, (1, 0, 2))
mesh = numpy.zeros(
(self.size(), 2, latbounds.shape[-1]), latbounds.dtype.char)
mesh[:, LAT, :] = numpy.reshape(
latbounds, (self.size(), latbounds.shape[-1]))
mesh[:, LON, :] = numpy.reshape(
lonbounds, (self.size(), latbounds.shape[-1]))
self._mesh_ = mesh
return self._mesh_
[docs] def _getShape(self):
"""
Not documented
"""
return self._lataxis_.shape
# Don't try to generate bounds for curvilinear grids
[docs] def genBounds(self):
"""
Not documented
"""
return (None, None)
# Get the n-th index axis. naxis is 0 or 1.
[docs] def getAxis(self, naxis):
"""
Not documented
"""
return self._lataxis_.getAxis(naxis)
[docs] def getMask(self):
"""Get the mask array, if any, otherwise None is returned."""
if self._maskVar_ is None:
return self._tempMask_
else:
return self._maskVar_[:]
[docs] def size(self):
"""
Not documented
"""
return self._lataxis_.size()
[docs] def writeScrip(self, cufile, gridTitle=None):
"""Write a grid to a SCRIP file.
Parameter
---------
cufile : is a Cdunif file, NOT a CDMS file.
gridtitle : is a string identifying the grid.
"""
lat = numpy.ma.filled(self._lataxis_[:])
lon = numpy.ma.filled(self._lonaxis_[:])
blat, blon = self.getBounds()
mask = self.getMask()
ni, nj = self.shape
if mask is None:
mask = numpy.ones((ni, nj), numpy.int32)
else:
tmp = 1 - mask
mask[:] = tmp.astype(mask.dtype.char)
mask = mask.astype(numpy.int32)
ngrid = ni * nj
centerLat = copy.copy(lat)
centerLat.shape = (ngrid,)
centerLon = copy.copy(lon)
centerLon.shape = (ngrid,)
mask.shape = (ngrid,)
clat = numpy.ma.filled(copy.copy(blat))
clat.shape = (ngrid, 4)
clon = numpy.ma.filled(copy.copy(blon))
clon.shape = (ngrid, 4)
# Write the file
if gridTitle is None:
gridTitle = self.id
cufile.title = gridTitle
cufile.createDimension("grid_size", ngrid)
cufile.createDimension("grid_corners", 4)
cufile.createDimension("grid_rank", 2)
griddims = cufile.createVariable("grid_dims", 'i', ("grid_rank",))
gridcenterlat = cufile.createVariable(
"grid_center_lat", 'd', ("grid_size",))
gridcenterlat.units = "degrees"
gridcenterlon = cufile.createVariable(
"grid_center_lon", 'd', ("grid_size",))
gridcenterlon.units = "degrees"
gridimask = cufile.createVariable("grid_imask", 'i', ("grid_size",))
gridimask.units = "unitless"
gridcornerlat = cufile.createVariable(
"grid_corner_lat", 'd', ("grid_size", "grid_corners"))
gridcornerlat.units = "degrees"
gridcornerlon = cufile.createVariable(
"grid_corner_lon", 'd', ("grid_size", "grid_corners"))
gridcornerlon.units = "degrees"
griddims[:] = numpy.array([nj, ni], numpy.int32)
gridcenterlat[:] = centerLat
gridcenterlon[:] = centerLon
gridimask[:] = mask
gridcornerlat[:] = clat
gridcornerlon[:] = clon
[docs] def toGenericGrid(self, gridid=None):
"""
Not documented
"""
from .auxcoord import TransientAuxAxis1D
from .coord import TransientVirtualAxis
from .gengrid import TransientGenericGrid
lat = numpy.ma.filled(self._lataxis_[:])
latunits = self._lataxis_.units
lon = numpy.ma.filled(self._lonaxis_[:])
lonunits = self._lonaxis_.units
blat, blon = self.getBounds()
mask = self.getMask()
ni, nj = self.shape
ngrid = ni * nj
centerLat = copy.copy(lat)
centerLat.shape = (ngrid,)
centerLon = copy.copy(lon)
centerLon.shape = (ngrid,)
if mask is not None:
mask.shape = (ngrid,)
cornerLat = numpy.ma.filled(copy.copy(blat))
cornerLat.shape = (ngrid, 4)
cornerLon = numpy.ma.filled(copy.copy(blon))
cornerLon.shape = (ngrid, 4)
iaxis = TransientVirtualAxis("cell", ngrid)
lataxis = TransientAuxAxis1D(centerLat, axes=(iaxis,), bounds=cornerLat,
attributes={'units': latunits}, id="latitude")
lonaxis = TransientAuxAxis1D(centerLon, axes=(iaxis,), bounds=cornerLon,
attributes={'units': lonunits}, id="longitude")
grid = TransientGenericGrid(lataxis, lonaxis, id=gridid, tempmask=mask)
return grid
[docs] def toCurveGrid(self, gridid=None):
"""
Not documented
"""
if gridid is None:
gridid = self.id
result = self.clone()
result.id = gridid
return result
[docs] def writeToFile(self, file):
"""
Not documented
"""
latvar = self._lataxis_.writeToFile(file)
lonvar = self._lonaxis_.writeToFile(file)
if self._maskVar_ is not None:
maskid = "mask_" + self.id
file.write(self._maskVar_, id=maskid)
latvar.maskid = maskid
lonvar.maskid = maskid
return (latvar, lonvar)
[docs] def writeg(self, file):
"""Write self as a Gridspec file representing a curvilinear grid.
The file, normally a CdmsFile, should already be open for writing
and will be closed."""
import time
# Set attributes
if (hasattr(file, 'Conventions')):
if (file.Conventions.find('Gridspec') < 0):
file.Conventions = file.Conventions + ' Gridspec-0.0'
else:
file.Conventions = 'Gridspec-0.0'
if (hasattr(file, 'gs_filetypes')):
if (file.gs_filetypes.find('Curvilinear_Tile') < 0):
file.gs_filetypes = file.gs_filetypes + ' Curvilinear_Tile'
else:
file.gs_filetypes = 'Curvilinear_Tile'
t = time.time()
id = int((t - int(t)) * 1.0e9)
file.gs_id = id
file.gs_originalfilename = os.path.basename(file.id)
newhistory = '\n' + time.ctime() + ' CDAT/CDMS AbstractCurveGrid'
# ... The \n gives a newline in the CDMS Python and in the Cdunif C code
# which gets called to write to a file. Someplace in the NetCDF libraries,
# or possibly the ncdump utility, the newline is converted back to a "\n"
# string, so you don't see a newline when you view the file. If we want
# a real newline as the CF specification says, the libraries or ncdump
# will have to be changed.
# ... In the future we may want to add more history information.
file.history = getattr(self, 'history', '') + \
getattr(file, 'history', '') + newhistory
# former tile variable and attributes
if (hasattr(self, 'long_name') and self.long_name is not None):
file.long_name = self.long_name
else:
file.long_name = 'gridspec_tile'
# gs_geometryType is no longer required of Gridspec files, but as yet
# there is no other proposal for describing the geometry (July 2010)
if (hasattr(self, 'gs_geometryType') and
self.gs_geometryType is not None):
file.gs_geometryType = self.gs_geometryType
else:
file.gs_geometryType = 'spherical'
# gs_discretizationType is no longer required of Gridspec files, but it's
# harmless and may come in useful
if (hasattr(self, 'gs_discretizationType') and
self.gs_discretizationType is not None):
file.gs_discretizationType = self.gs_discretizationType
else:
file.gs_discretizationType = 'logically_rectangular'
file.gs_lonv = 'gs_x'
file.gs_latv = 'gs_y'
# Set up and write variables. When written, cdms writes not only the arrays
# but also their coordinates, e.g. gs_nip.
x = self._lonaxis_
if (not hasattr(x, 'units')):
print("Warning, no units found for longitude")
x.units = 'degree_east'
if (not hasattr(x, 'standard_name')):
print("Warning, no standard_name found for longitude axis")
x.standard_name = 'longitude'
if (x.standard_name == 'geographic_longitude'):
# temporary for updating test files
x.standard_name = 'longitude'
x.id = file.gs_lonv
# _lonaxis_ is a TransientAxis2D, hence a TransientVariable
# But I don't know where the attribute _TransientVariable__domain comes
# from
y = self._lataxis_
if (not hasattr(y, 'units')):
print("Warning, no units found for latitude")
y.units = 'degree_north'
if (not hasattr(y, 'standard_name')):
print("Warning, no standard_name found for latitude axis")
y.standard_name = 'latitude'
if (y.standard_name == 'geographic_latitude'):
# temporary for updating test files
y.standard_name = 'latitude'
y.id = file.gs_latv
if(not hasattr(x, '_TransientVariable__domain')):
# There probably doesn't exist enough information to write a correct
# grid, but this will help.
x._TransientVariable__domain = [(x,), (y,)]
x._TransientVariable__domain[0][0].id = 'gs_njp'
x._TransientVariable__domain[1][0].id = 'gs_nip'
if (not hasattr(y, '_TransientVariable__domain')):
# There probably doesn't exist enough information to write a correct
# grid, but this will help.
y._TransientVariable__domain = [(x,), (y,)]
y._TransientVariable__domain[0][0].id = 'gs_njp'
y._TransientVariable__domain[1][0].id = 'gs_nip'
file.write(x)
file.write(y)
file.close()
[docs] def write_gridspec(self, filename):
"""Writes this grid to a Gridspec-compliant file, or does nothing if there is
already a known file corresponding to this grid. The filename should be a
complete path."""
# This method was never completed because the libCF functionality I had planned to
# use never appeared (yet).
# The functionality (other than checking gsfile) is now done by the writeg
# method above.
if (not hasattr(self, "gsfile")):
self.gsfile = None
self.gspath = None
if (self.gsfile is not None):
return (self.gsfile, self.gspath)
else:
raise RuntimeError(
'The libCF/Gridspec API does not provide for writing CurveGrids<<<')
[docs] def init_from_gridspec(self, filename):
"""Reads to grid from a Gridspec-compliant file. The filename should be a
complete path. The contents of the file may overwrite data in the existing
grid object."""
# - This is really a kind of init function. The __init__ function should
# determine what kind of initialization is being done (from a file, from
# another object, from arguments specifying the contents e.g. axes) and branch
# to call a method such as this one.
# - Another way to read a file is with the standard CDMS
# pattern file=cdms2.open(...); g=file('grid') or g=file('')
try:
f = cdms2.open(filename)
except IOError:
print("Cannot open grid file for reading: ", filename)
return
self.init_from_gridspec_file(self, f)
f.close()
[docs] def init_from_gridspec_file(self, f):
"""Reads to grid from a Gridspec-compliant file, f. This f should be a
CdmsFile object, already open for reading. The contents of the file may
overwrite data in the existing grid object."""
# As for the above init_from_gridspec method, this is really a kind of
# init function and should be called from __init__ .
ax, ay, gs_attr = f.gridspec_file_contents()
# ... gridspec_file_contents is defined in dataset.py
self.__init__(ay, ax)
self.gsfile = gs_attr['filebase']
self.gspath = gs_attr['filepath']
self.long_name = gs_attr['long_name']
# gs_geometryType is no longer required of Gridspec files, but as yet
# there is no other proposal for describing the geometry (July 2010)
self.gs_geometryType = gs_attr['gs_geometryType']
# gs_discretizationType is no longer required of Gridspec files, but it's
# harmless and may come in useful
self.gs_discretizationType = gs_attr['gs_discretizationType']
return self
[docs] def subSlice(self, *specs, **keys):
"""Get a transient subgrid based on an argument list <specs> of slices."""
newlat = self._lataxis_.subSlice(*specs, **keys)
newlon = self._lonaxis_.subSlice(*specs, **keys)
if self._maskVar_ is None:
newmask = None
else:
newmask = self._maskVar_.subSlice(*specs, **keys)
result = TransientCurveGrid(newlat, newlon, maskvar=newmask)
return result
[docs] def getGridSlices(self, domainlist, newaxislist, slicelist):
"""Determine which slices in slicelist correspond to the lat/lon elements
of the grid.
Parameters
----------
domainlist : is a list of axes of a variable.
newaxislist : is a list of result axes after the slicelist is applied to domainlist.
slicelist : is a list of slices.
All lists are of equal length.
Returns
-------
value : is (newslicelist, gridaxislist) where
newslicelist : is the elements of slicelist that correspond to the grid, in the
preferred order of the grid.
gridaxislist : is the elements of newaxislist that correspond to the grid, in the
preferred order of the grid.
"""
iaxis = self._lataxis_.getAxis(0)
jaxis = self._lataxis_.getAxis(1)
k = 0
i = j = -1
for d in domainlist:
if d.shape == iaxis.shape:
if numpy.allclose(d[:], iaxis[:]) is True:
inewaxis = newaxislist[k]
islice = slicelist[k]
i = k
if d.shape == jaxis.shape:
if numpy.allclose(d[:], jaxis[:]) is True:
jnewaxis = newaxislist[k]
jslice = slicelist[k]
j = k
k += 1
if i == -1 or j == -1:
raise RuntimeError(
'Grid lat/lon domains do not match variable domain')
return ((islice, jslice), (inewaxis, jnewaxis))
[docs] def getIndex(self):
"""Get the grid index"""
if self._index_ is None:
# Trying to stick in Stephane Raynaud's patch for autodetection
nj, ni = self._lataxis_.shape
dlon = numpy.max(self._lonaxis_[:]) - numpy.min(self._lonaxis_[:])
dx = max(dlon / ni, dlon / nj)
dlat = numpy.max(self._lataxis_[:]) - numpy.min(self._lataxis_[:])
dy = max(dlat / ni, dlat / nj)
latlin = numpy.ravel(numpy.ma.filled(self._lataxis_[:]))
lonlin = numpy.ravel(numpy.ma.filled(self._lonaxis_[:]))
_bindex.setDeltas(dx, dy)
self._index_ = bindex.bindexHorizontalGrid(latlin, lonlin)
return self._index_
[docs] def intersect(self, spec):
"""Intersect with the region specification.
Parameters
----------
'spec' : is a region specification of the form defined in the grid module.
Returns
-------
(mask, indexspecs) where
'mask' : is the mask of the result grid AFTER self and region spec are interested.
'indexspecs' : is a list of index specifications suitable for slicing a
variable with the given grid.
"""
ni, nj = self.shape
index = self.getIndex()
latspec = spec[CoordTypeToLoc[LatitudeType]]
lonspec = spec[CoordTypeToLoc[LongitudeType]]
latlin = numpy.ravel(numpy.ma.filled(self._lataxis_[:]))
lonlin = numpy.ravel(numpy.ma.filled(self._lonaxis_[:]))
points = bindex.intersectHorizontalGrid(
latspec, lonspec, latlin, lonlin, index)
if len(points) == 0:
raise CDMSError(
'No data in the specified region, longitude=%s, latitude=%s' %
(repr(lonspec), repr(latspec)))
fullmask = numpy.ones(ni * nj)
numpy.put(fullmask, points, 0)
fullmask = numpy.reshape(fullmask, (ni, nj))
iind = points // nj
jind = points - iind * nj
imin, imax, jmin, jmax = (
min(iind), max(iind) + 1, min(jind), max(jind) + 1)
submask = fullmask[imin:imax, jmin:jmax]
yid = self.getAxis(0).id
xid = self.getAxis(1).id
indexspecs = {yid: slice(imin, imax), xid: slice(jmin, jmax)}
return submask, indexspecs
def getAxisList(self):
"""
Not documented
"""
return (self._lataxis_.getAxis(0), self._lataxis_.getAxis(1))
[docs] def isClose(self, g):
"""
Is Close
Returns
-------
1 if g : is a grid of the same type and shape. A real element-by-element
comparison would be too expensive here."""
if g is None:
return 0
elif self.shape != g.shape:
return 0
elif not isinstance(g, AbstractCurveGrid):
return 0
else:
return 1
[docs] def checkAxes(self, axes):
"""
Check Axes
Returns
-------
1 iff every element of self.getAxisList() is in the list 'axes'."""
for item in self.getAxisList():
# if all [False, False, ....] result=0
if not any([allclose(item[:], axis[:]) for axis in axes]):
result = 0
break
else:
result = 1
return result
[docs] def reconcile(self, axes):
"""
Reconcile
Returns
-------
a grid that is consistent with the axes, or None.
For curvilinear grids this means that the grid-related axes are
contained in the 'axes' list.
"""
result = self
selfaxes = self.getAxisList()
missing = []
for i in range(2):
if selfaxes[i] not in axes:
missing.append(i)
result = None
# Some of the grid axes are not in the 'axes' list
if result is None:
result = self.clone()
used = [] # axes already matched
for i in missing:
for item in axes:
if (item not in used) and len(selfaxes[i]) == len(
item) and allclose(selfaxes[i], item):
result._lataxis_.setAxis(i, item)
result._lonaxis_.setAxis(i, item)
used.append(item)
break
else:
result = None
break
return result
[docs] def flatAxes(self):
"""
Flat Axes
Returns
-------
(flatlat, flatlon) where flatlat is a 1D NumPy array
having the same length as the number of cells in the grid, similarly for flatlon.
"""
if self._flataxes_ is None:
from . import MV2 as MV
alat = MV.filled(self.getLatitude())
alon = MV.filled(self.getLongitude())
alatflat = numpy.ravel(alat)
alonflat = numpy.ravel(alon)
self._flataxes_ = (alatflat, alonflat)
return self._flataxes_
shape = property(_getShape, None)
# PropertiedClasses.set_property (AbstractCurveGrid, 'shape',
# AbstractCurveGrid._getShape, nowrite=1,
# nodelete=1)
class DatasetCurveGrid(AbstractCurveGrid):
def __init__(self, latAxis, lonAxis, id, parent=None,
maskvar=None, tempmask=None, node=None):
"""Create a file curvilinear grid.
"""
AbstractCurveGrid.__init__(
self, latAxis, lonAxis, id, maskvar, tempmask, node)
self.parent = parent
def __repr__(self):
return "<DatasetCurveGrid, id: %s, shape: %s>" % (
self.id, repr(self.shape))
class FileCurveGrid(AbstractCurveGrid):
"""
Not documented
"""
def __init__(self, latAxis, lonAxis, id, parent=None,
maskvar=None, tempmask=None, node=None):
"""Create a file curvilinear grid.
"""
AbstractCurveGrid.__init__(
self, latAxis, lonAxis, id, maskvar, tempmask, node)
self.parent = parent
def __repr__(self):
return "<FileCurveGrid, id: %s, shape: %s>" % (
self.id, repr(self.shape))
class TransientCurveGrid(AbstractCurveGrid):
"""
Not documented
"""
grid_count = 0
[docs] def __init__(self, latAxis, lonAxis, id=None, maskvar=None, tempmask=None):
"""Create a file curvilinear grid.
"""
if id is None:
TransientCurveGrid.grid_count += 1
id = 'grid_' + str(TransientCurveGrid.grid_count)
AbstractCurveGrid.__init__(
self, latAxis, lonAxis, id, maskvar, tempmask)
def __repr__(self):
return "<TransientCurveGrid, id: %s, shape: %s>" % (
self.id, repr(self.shape))
[docs] def toCurveGrid(self, gridid=None):
"""
Not documented
"""
if gridid is None:
result = self
else:
result = self.clone()
result.id = gridid
return result
[docs]def readScripCurveGrid(fileobj, dims, whichType, whichGrid):
"""Read a 'native' SCRIP grid file, returning a transient curvilinear grid.
Parameters
----------
fileobj : is an open CDMS dataset or file object.
dims : is the grid shape.
whichType : is the type of file, either "grid" or "mapping"
f whichType is "mapping", whichGrid is the choice of grid, either "source" or "destination"
Returns
-------
"""
from .coord import TransientAxis2D
if 'S' in list(fileobj.variables.keys()):
if whichType == "grid":
gridCornerLatName = 'grid_corner_lat'
gridCornerLonName = 'grid_corner_lon'
gridMaskName = 'grid_imask'
gridCenterLatName = 'grid_center_lat'
gridCenterLonName = 'grid_center_lon'
titleName = 'title'
elif whichGrid == "destination":
gridCornerLatName = 'yv_b'
gridCornerLonName = 'xv_b'
gridMaskName = 'mask_b'
gridCenterLatName = 'yc_b'
gridCenterLonName = 'xc_b'
titleName = 'dest_grid'
else:
gridCornerLatName = 'yv_a'
gridCornerLonName = 'xv_a'
gridMaskName = 'mask_a'
gridCenterLatName = 'yc_a'
gridCenterLonName = 'xc_a'
titleName = 'source_grid'
else:
if whichType == "grid":
gridCornerLatName = 'grid_corner_lat'
gridCornerLonName = 'grid_corner_lon'
gridMaskName = 'grid_imask'
gridCenterLatName = 'grid_center_lat'
gridCenterLonName = 'grid_center_lon'
titleName = 'title'
elif whichGrid == "destination":
gridCornerLatName = 'dst_grid_corner_lat'
gridCornerLonName = 'dst_grid_corner_lon'
gridMaskName = 'dst_grid_imask'
gridCenterLatName = 'dst_grid_center_lat'
gridCenterLonName = 'dst_grid_center_lon'
titleName = 'dest_grid'
else:
gridCornerLatName = 'src_grid_corner_lat'
gridCornerLonName = 'src_grid_corner_lon'
gridMaskName = 'src_grid_imask'
gridCenterLatName = 'src_grid_center_lat'
gridCenterLonName = 'src_grid_center_lon'
titleName = 'source_grid'
vardict = fileobj.variables
cornerLat = fileobj(gridCornerLatName)
cornerLon = fileobj(gridCornerLonName)
ncorners = cornerLat.shape[-1]
ni = dims[1]
nj = dims[0]
gridshape = (ni, nj)
boundsshape = (ni, nj, ncorners)
if hasattr(cornerLat, 'units') and cornerLat.units[0:6].lower() == 'radian':
cornerLat = (cornerLat * (180.0 / numpy.pi)).reshape(boundsshape)
cornerLon = (cornerLon * (180.0 / numpy.pi)).reshape(boundsshape)
else:
cornerLat = cornerLat.reshape(boundsshape)
cornerLon = cornerLon.reshape(boundsshape)
iaxis = TransientVirtualAxis("i", ni)
jaxis = TransientVirtualAxis("j", nj)
if gridMaskName in vardict:
# SCRIP convention: 0 for invalid data
# numpy.ma convention: 1 for invalid data
mask = 1 - fileobj(gridMaskName)
mask = mask.reshape(gridshape)
else:
mask = None
if gridCenterLatName in vardict:
centerLat = fileobj(gridCenterLatName).reshape(gridshape)
gclat = fileobj[gridCenterLatName]
if hasattr(gclat, "units") and gclat.units.lower() == 'radians':
centerLat *= (180.0 / numpy.pi)
else:
centerLat = cornerLat[:, :, 0]
if gridCenterLonName in vardict:
centerLon = fileobj(gridCenterLonName).reshape(gridshape)
gclon = fileobj[gridCenterLonName]
if hasattr(gclon, "units") and gclon.units.lower() == 'radians':
centerLon *= (180.0 / numpy.pi)
else:
centerLon = cornerLon[:, :, 0]
if hasattr(fileobj, titleName):
gridid = getattr(fileobj, titleName)
gridid = gridid.strip().replace(' ', '_')
else:
gridid = "<None>"
lataxis = TransientAxis2D(centerLat, axes=(iaxis, jaxis), bounds=cornerLat,
attributes={'units': 'degrees_north'}, id="latitude")
lonaxis = TransientAxis2D(centerLon, axes=(iaxis, jaxis), bounds=cornerLon,
attributes={'units': 'degrees_east'}, id="longitude")
grid = TransientCurveGrid(lataxis, lonaxis, id=gridid, tempmask=mask)
return grid
