Changeset 26569

Timestamp:

11/09/21 06:58:37 (3 years ago)

Author:

bdef

Message:

CHG:modification to enable AMR plotting

Location:

issm/trunk-jpl/src/m/plot

Files:

• colormaps/demmap.py (modified) (1 diff)
• plotmodel.py (modified) (7 diffs)
• processdata.py (modified) (4 diffs)
• processmesh.py (modified) (2 diffs)

issm/trunk-jpl/src/m/plot/colormaps/demmap.py

@@ -64 +64 @@
         cmx = nland * interval
 
-    clim = [cmn, cmx]
-
     if strcmpi(colorscheme, 'dem'):
         # concatenate and transpose to match matplotlib's colormap format

issm/trunk-jpl/src/m/plot/plotmodel.py

@@ -14 +14 @@
 def plotmodel(md, *args):
     '''
-    PLOTMODEL - At command prompt, type 'plotdoc()' for additional 
+    PLOTMODEL - At command prompt, type 'plotdoc()' for additional
     documentation.
 
@@ -31 +31 @@
     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 +60 @@
     # 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 +83 @@
             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 +95 @@
         #
         # 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 +118 @@
         #   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 +127 @@
             direction=direction,
             axes_pad=axes_pad,
-            add_all=add_all,
             share_all=share_all,
             label_mode=label_mode,

issm/trunk-jpl/src/m/plot/processdata.py

@@ -15 +15 @@
     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 +22 @@
     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 +51 @@
         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 +148 @@
             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")

issm/trunk-jpl/src/m/plot/processmesh.py

@@ -41 +41 @@
     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 +75 @@
         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')