Index: /issm/trunk-jpl/src/m/plot/colormaps/demmap.py =================================================================== --- /issm/trunk-jpl/src/m/plot/colormaps/demmap.py (revision 26568) +++ /issm/trunk-jpl/src/m/plot/colormaps/demmap.py (revision 26569) @@ -64,6 +64,4 @@ cmx = nland * interval - clim = [cmn, cmx] - if strcmpi(colorscheme, 'dem'): # concatenate and transpose to match matplotlib's colormap format Index: /issm/trunk-jpl/src/m/plot/plotmodel.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plotmodel.py (revision 26568) +++ /issm/trunk-jpl/src/m/plot/plotmodel.py (revision 26569) @@ -14,5 +14,5 @@ def plotmodel(md, *args): ''' - PLOTMODEL - At command prompt, type 'plotdoc()' for additional + PLOTMODEL - At command prompt, type 'plotdoc()' for additional documentation. @@ -31,5 +31,5 @@ subplotwidth = ceil(sqrt(numberofplots)) - # TODO: Check that commenting this out is correct; we do not need a hold + # TODO: Check that commenting this out is correct; we do not need a hold # under matplotlib, right? # @@ -60,5 +60,5 @@ # Go through plots # - # NOTE: The following is where Python + matplolib differs substantially in + # NOTE: The following is where Python + matplolib differs substantially in # implementation and inteface from MATLAB. # @@ -83,5 +83,5 @@ plotnum = None - # NOTE: The inline comments for each of the following parameters are + # NOTE: The inline comments for each of the following parameters are # taken from https://matplotlib.org/api/_as_gen/mpl_toolkits.axes_grid1.axes_grid.ImageGrid.html # @@ -95,5 +95,5 @@ # # TODO: - # - Add 'edge' option (research if there is a corresponding option in + # - Add 'edge' option (research if there is a corresponding option in # MATLAB)? # @@ -118,5 +118,5 @@ # rect(float, float, float, float) or int # - # The axes position, as a (left, bottom, width, height) tuple or as a + # The axes position, as a (left, bottom, width, height) tuple or as a # three-digit subplot position code (e.g., "121"). # @@ -127,5 +127,4 @@ direction=direction, axes_pad=axes_pad, - add_all=add_all, share_all=share_all, label_mode=label_mode, Index: /issm/trunk-jpl/src/m/plot/processdata.py =================================================================== --- /issm/trunk-jpl/src/m/plot/processdata.py (revision 26568) +++ /issm/trunk-jpl/src/m/plot/processdata.py (revision 26569) @@ -15,5 +15,5 @@ See also: PLOTMODEL, PROCESSMESH """ - # {{{ Initialization and grabbing auxiliaries + # Initialization and grabbing auxiliaries {{{ # check format if (len(data) == 0 or (len(data) == 1 and not isinstance(data, dict) and np.isnan(data).all())): @@ -22,8 +22,15 @@ if 'numberofvertices2d' in dir(md.mesh): numberofvertices2d = md.mesh.numberofvertices2d - numberofelements2d = md.mesh.numberofelements2d else: numberofvertices2d = np.nan - numberofelements2d = np.nan + + if options.exist('amr'): + step = options.getfieldvalue('amr', 0) + numberofvertices = len(md.results.TransientSolution[step].MeshX) + numberofelements = np.shape(md.results.TransientSolution[step].MeshElements)[0] + else: + numberofvertices = md.mesh.numberofvertices + numberofelements = md.mesh.numberofelements + procdata = np.copy(data) #initialize datatype @@ -44,5 +51,93 @@ raise ValueError('data passed to plotmodel has bad dimensions; check that column vectors are rank - 1') # }}} - # {{{ process NaN's if any + + # log {{{ + if options.exist('log'): + cutoff = options.getfieldvalue('log', 1) + procdata[np.where(procdata < cutoff)] = cutoff + # }}} + + # quiver plot {{{ + if datasize[1] > 1 and datasize[0] != numberofvertices + 1: + if datasize[0] == numberofvertices and datasize[1] == 2: + datatype = 3 + else: + raise ValueError('plotmodel error message: data should have two columns of length md.mesh.numberofvertices for a quiver plot') + # }}} + + # element data{{{ + if datasize[0] == numberofelements and datasize[1] == 1: + #initialize datatype if non patch + if datatype != 4 and datatype != 5: + datatype = 1 + # AMR {{{ + if options.exist('amr'): + nonan = np.nonzero(~np.isnan(md.results.TransientSolution[step].MeshElements)) + procdata = procdata[nonan] + # }}} + # mask {{{ + if options.exist('mask'): + flags = options.getfieldvalue('mask') + hide = np.invert(flags) + if np.size(flags) == numberofvertices: + EltMask = np.asarray([np.any(np.in1d(index, np.where(hide))) for index in md.mesh.elements - 1]) + procdata = np.ma.array(procdata, mask=EltMask) + options.addfielddefault('cmap_set_bad', 'w') + elif np.size(flags) == numberofelements: + procdata = np.ma.array(procdata, mask=hide) + options.addfielddefault('cmap_set_bad', 'w') + else: + print('Warning: processdata.py: mask length not supported yet (supported length are md.mesh.numberofvertices and md.mesh.numberofelements') + # }}} + + # }}} + + # node data {{{ + if datasize[0] in [numberofvertices, numberofvertices2d] and datasize[1] == 1: + datatype = 2 + # AMR {{{ + if options.exist('amr'): + nonan = np.nonzero(~np.isnan(md.results.TransientSolution[step].MeshX)) + procdata = procdata[nonan] + # }}} + # Mask {{{ + if options.exist('mask'): + flags = options.getfieldvalue('mask') + hide = np.invert(flags) + if np.size(flags) == numberofvertices: + procdata = np.ma.array(procdata, mask=hide) + options.addfielddefault('cmap_set_bad', 'w') + elif np.size(flags) == numberofelements: + NodeMask = np.zeros(np.shape(md.mesh.x), dtype=bool) + HideElt = md.mesh.elements[np.where(hide)[0]] - 1 + NodeMask[HideElt] = True + procdata = np.ma.array(procdata, mask=NodeMask) + options.addfielddefault('cmap_set_bad', 'w') + else: + print('Warning: processdata.py: mask length not supported yet (supported length are md.mesh.numberofvertices and md.mesh.numberofelements') + # }}} + # }}} + + # spc time series {{{ + if datasize[0] == numberofvertices + 1: + datatype = 2 + spccol = options.getfieldvalue('spccol', 0) + print('multiple-column spc field; specify column to plot using option "spccol"') + print(('column ', spccol, ' plotted for time: ', procdata[-1, spccol])) + procdata = procdata[0:-1, spccol] + + #mask? + + #layer projection? + + #control arrow density if quiver plot + # }}} + + # convert rank - 2 array to rank - 1 {{{ + if np.ndim(procdata) == 2 and np.shape(procdata)[1] == 1: + procdata = procdata.reshape(-1, ) + # }}} + + # process NaN's if any {{{ nanfill = options.getfieldvalue('nan', -9999) if np.any(np.isnan(procdata)): @@ -53,83 +148,13 @@ ub = ub + 0.5 nanfill = lb - 1 - #procdata[np.isnan(procdata)] = nanfill + procdata[np.isnan(procdata)] = nanfill procdata = np.ma.array(procdata, mask=np.isnan(procdata)) - options.addfielddefault('clim', [lb, ub]) + #clim looks to be deprecated and replaced by caxis + #options.addfielddefault('clim', [lb, ub]) options.addfielddefault('cmap_set_under', '1') print(("WARNING: nan's treated as", nanfill, "by default. Change using pairoption 'nan', nan_fill_value in plotmodel call")) # }}} - # {{{ log - if options.exist('log'): - cutoff = options.getfieldvalue('log', 1) - procdata[np.where(procdata < cutoff)] = cutoff - # }}} - # {{{ quiver plot - if datasize[1] > 1 and datasize[0] != md.mesh.numberofvertices + 1: - if datasize[0] == md.mesh.numberofvertices and datasize[1] == 2: - datatype = 3 - else: - raise ValueError('plotmodel error message: data should have two columns of length md.mesh.numberofvertices for a quiver plot') - # }}} - # {{{ element data - if datasize[0] == md.mesh.numberofelements and datasize[1] == 1: - #initialize datatype if non patch - if datatype != 4 and datatype != 5: - datatype = 1 - # {{{mask - if options.exist('mask'): - flags = options.getfieldvalue('mask') - hide = np.invert(flags) - if np.size(flags) == md.mesh.numberofvertices: - EltMask = np.asarray([np.any(np.in1d(index, np.where(hide))) for index in md.mesh.elements - 1]) - procdata = np.ma.array(procdata, mask=EltMask) - options.addfielddefault('cmap_set_bad', 'w') - elif np.size(flags) == md.mesh.numberofelements: - procdata = np.ma.array(procdata, mask=hide) - options.addfielddefault('cmap_set_bad', 'w') - else: - print('Warning: processdata.py: mask length not supported yet (supported length are md.mesh.numberofvertices and md.mesh.numberofelements') - # }}} - - # }}} - # {{{ node data - if datasize[0] == md.mesh.numberofvertices and datasize[1] == 1: - datatype = 2 - # {{{ Mask - if options.exist('mask'): - flags = options.getfieldvalue('mask') - hide = np.invert(flags) - if np.size(flags) == md.mesh.numberofvertices: - procdata = np.ma.array(procdata, mask=hide) - options.addfielddefault('cmap_set_bad', 'w') - elif np.size(flags) == md.mesh.numberofelements: - NodeMask = np.zeros(np.shape(md.mesh.x), dtype=bool) - HideElt = md.mesh.elements[np.where(hide)[0]] - 1 - NodeMask[HideElt] = True - procdata = np.ma.array(procdata, mask=NodeMask) - options.addfielddefault('cmap_set_bad', 'w') - else: - print('Warning: processdata.py: mask length not supported yet (supported length are md.mesh.numberofvertices and md.mesh.numberofelements') - # }}} - # }}} - # {{{ spc time series - if datasize[0] == md.mesh.numberofvertices + 1: - datatype = 2 - spccol = options.getfieldvalue('spccol', 0) - print('multiple-column spc field; specify column to plot using option "spccol"') - print(('column ', spccol, ' plotted for time: ', procdata[-1, spccol])) - procdata = procdata[0:-1, spccol] - - #mask? - - #layer projection? - - #control arrow density if quiver plot - # }}} - # {{{ convert rank - 2 array to rank - 1 - if np.ndim(procdata) == 2 and np.shape(procdata)[1] == 1: - procdata = procdata.reshape(-1, ) - # }}} - # {{{ if datatype is still zero, error out + # if datatype is still zero, error out {{{ if datatype == 0: raise ValueError("processdata error: data provided not recognized or not supported") Index: /issm/trunk-jpl/src/m/plot/processmesh.py =================================================================== --- /issm/trunk-jpl/src/m/plot/processmesh.py (revision 26568) +++ /issm/trunk-jpl/src/m/plot/processmesh.py (revision 26569) @@ -41,12 +41,22 @@ if hasattr(md.mesh, 'elements2d'): elements2d = md.mesh.elements2d + numofvertices2d = md.mesh.numberofvertices2d + numofelements2d = md.mesh.numberofelements2d + else: + numofvertices2d = np.nan + numofelements2d = np.nan if options.exist('amr'): step = options.getfieldvalue('amr') - x = md.results.TransientSolution[step].MeshX - y = md.results.TransientSolution[step].MeshY - elements = md.results.TransientSolution[step].MeshElements + nonan = np.nonzero(~np.isnan(md.results.TransientSolution[step].MeshX)) + x = md.results.TransientSolution[step].MeshX[nonan] + y = md.results.TransientSolution[step].MeshY[nonan] + nonan = np.nonzero(~np.isnan(md.results.TransientSolution[step].MeshElements)) + elements = md.results.TransientSolution[step].MeshElements[nonan] - 1 + eldim = np.shape(md.results.TransientSolution[step].MeshElements)[1] + elements = np.reshape(elements, ((int(len(elements) / eldim), eldim))) + else: - elements = md.mesh.elements + elements = md.mesh.elements - 1 if options.getfieldvalue('coord', 'xy') != 'latlon': x = md.mesh.x @@ -65,14 +75,13 @@ z = getattr(md, z) + force2D = numofelements2d in np.shape(data) or numofvertices2d in np.shape(data) #is it a 2D plot? - if md.mesh.dimension() == 2 or options.getfieldvalue('layer', 0) >= 1: + if md.mesh.dimension() == 2 or options.getfieldvalue('layer', 0) >= 1 or force2D: is2d = 1 else: is2d = 0 - elements = md.mesh.elements - 1 - #layer projection? - if options.getfieldvalue('layer', 0) >= 1: + if options.getfieldvalue('layer', 0) >= 1 or force2D: if options.getfieldvalue('coord', 'xy') == 'latlon': raise Exception('processmesh error message: cannot work with 3D meshes for now')