Merge branch 'main' of github.com:SCECcode/ucvm_plotting

pycvm/__init__.py

@@ -10,10 +10,10 @@
 from .vs30_slice import Vs30Slice
 from .elevation_slice import ElevationSlice
 from .vs30_etree_slice import Vs30EtreeSlice
-from .vs30_etree_difference_slice import Vs30EtreeDifferenceSlice
+from .horizontal_difference_slice import HorizontalDifferenceSlice
+from .cross_difference_section import CrossDifferenceSection
 from .map_grid_horizontal_slice import MapGridHorizontalSlice
 from .basin_slice import BasinSlice, Z10Slice, Z25Slice
 from .depth_profile import DepthProfile
 from .elevation_profile import ElevationProfile
 from .difference import Difference
-from .cybershake import CyberShake

pycvm/common.py

@@ -962,7 +962,7 @@ def import_binary(self, fname, num_x, num_y):
         if( k != -1) : 
             rawfile = rawfile[:k] + "_data.bin"
         try :
-            fh = open(rawfile, 'r') 
+            fh = open(rawfile, 'rb') 
         except:
             print("ERROR: binary data does not exist.")
             exit(1)
@@ -973,16 +973,22 @@ def import_binary(self, fname, num_x, num_y):
         if len(floats) == 2 * (num_x * num_y) :
           fh.seek(0,0)
           floats = np.fromfile(fh, dtype=np.float)
+        fh.close()
+
+## maybe it is a np array
+        if len(floats) != (num_x * num_y) :
+          ffh = open(rawfile, 'rb') 
+          ffloats = np.load(ffh)
+          ffh.close()
+          floats=ffloats.reshape([num_x * num_y]);
 
         print("TOTAL number of binary data read:"+str(len(floats))+"\n")
 
         # sanity check,  
         if len(floats) != (num_x * num_y) :
-            print("import_binary(), wrong size !!!"+ str(len(floats))+ " expecting "+ str(num_x * num_y))
+            print("import_binary(), wrong size !!!"+ str(len(floats)) + " expecting "+ str(num_x * num_y))
             exit(1)
 
-        fh.close()
-
         if len(floats) == 1:
             return floats[0]
 
@@ -991,7 +997,7 @@ def import_binary(self, fname, num_x, num_y):
 #  export np float array to an exernal file
 #  
     def export_np_float_array(self, floats, fname):
-#        print("calling export_np_float_array")
+#       print("calling export_np_float_array -",len(floats))
         rawfile = fname
         if rawfile is None :
             rawfile="data.bin"
@@ -1005,11 +1011,13 @@ def export_np_float_array(self, floats, fname):
             exit(1)
 
         np.save(fh, floats)
-        fh.close
+        fh.close()
+
+
 
 #  export raw floats nxy ndarray  to an external file 
     def export_binary(self, floats, fname):
-#        print("calling export_binary")
+        print("calling export_binary -",len(floats))
         rawfile = fname
         if rawfile is None :
             rawfile="data.bin"
@@ -1044,7 +1052,7 @@ def import_metadata(self, fname):
 
 #  export ascii meta data to an external file 
     def export_metadata(self,meta,fname):
-#        print("calling export_metadata")
+        print("calling export_metadata")
         metafile=fname
         if metafile is None:
           metafile = "meta.json"
@@ -1075,7 +1083,7 @@ def import_matprops(self, fname):
 
 #  export material properties in JSON form to an external file 
     def export_matprops(self,blob,fname):
-#        print("calling export_matprops")
+        print("calling export_matprops")
         matpropsfile=fname
         if matpropsfile is None :
             matpropsfile="matprops.json"
@@ -1145,7 +1153,7 @@ def makebounds(self,minval=0.0,maxval=5.0,nstep=0,meanval=None, substep=5,all=Tr
         l=0
         nsubstep=substep
         nnstep=float(step)/nsubstep
-        if(meanval != None) :
+        if(meanval != None and step!= 0) :
           l =(meanval-minval) //step
 
         for i in range(0,nstep) :

pycvm/vs30_etree_difference_slice.py → pycvm/cross_difference_section.py

@@ -1,20 +1,20 @@
 ##
-#  @file vs30_etree_difference_slice.py
-#  @brief Take 2 slices of vs30 values, plot a difference horizontal slice
+#  @file cross_difference_section.py
+#  @brief Take 2 slices of cross sections, plot a difference plot
 #  @author Mei-Hui Su - SCEC
 #  @version 
 #
 #  Imports
-from horizontal_slice import HorizontalSlice
+from cross_section import CrossSection
 from common import Point, MaterialProperties, UCVM, UCVM_CVMS, \
                    math, pycvm_cmapDiscretize, cm, mcolors, basemap, np, plt
 
 ##
-#  @class Vs30EtreeSlice
-#  @brief Gets a horizontal slice of Vs30 data.
+#  @class CrossDifferencSection
+#  @brief Gets 2 cross section and make a difference plot
 #
-#  Retrieves 2 horizontal slices of Vs30 values and make a difference plot 
-class Vs30EtreeDifferenceSlice(HorizontalSlice):
+#  Retrieves 2 cross sections and make a difference plot 
+class CrossDifferenceSection(CrossSection):
 
     ##
     #  Initializes the super class and copies the parameters over.
@@ -26,8 +26,8 @@ class Vs30EtreeDifferenceSlice(HorizontalSlice):
     #  
     def __init__(self, upperleftpoint, bottomrightpoint, meta={}):
 
-        #  Initializes the base class which is a horizontal slice.
-        HorizontalSlice.__init__(self, upperleftpoint, bottomrightpoint, meta)
+        #  Initializes the base class which is a cross section.
+        CrossSection.__init__(self, upperleftpoint, bottomrightpoint, meta)
 
         if 'datafile1' in self.meta :
             self.datafile1 = self.meta['datafile1']
@@ -41,7 +41,7 @@ def __init__(self, upperleftpoint, bottomrightpoint, meta={}):
 
 
     ##
-    #  Retrieves the values for this Vs30 slice and stores them in the class.
+    #  Retrieves the values for this cross section and stores them in the class.
     def getplotvals(self, property="vs") :
 
         #  How many y and x values will we need?
@@ -61,7 +61,7 @@ def getplotvals(self, property="vs") :
         else :
            self.num_y = int(math.ceil(self.plot_height / self.spacing)) + 1
 
-        ## The 2D array of retrieved Vs30 values.
+        ## The 2D array of retrieved values.
         self.materialproperties = [[MaterialProperties(-1, -1, -1) for x in range(self.num_x)] for x in range(self.num_y)] 
 
         u = UCVM(install_dir=self.installdir, config_file=self.configfile)
@@ -102,7 +102,7 @@ def getplotvals(self, property="vs") :
 
         i = 0
         j = 0
-        
+
         for idx in range(len(dataA)) :
             self.materialproperties[i][j].vs = dataA[idx]-dataB[idx]
             j = j + 1
@@ -111,8 +111,8 @@ def getplotvals(self, property="vs") :
                 i = i + 1
 
     ##
-    #  Plots the Vs30 data as a horizontal slice. This code is very similar to the
-    #  HorizontalSlice routine.
+    #  Plots the Difference data as a cross section. This code is very similar to the
+    #  CrossSection routine.
     #
     #  @param filename The location to which the plot should be saved. Optional.
     #  @param title The title of the plot to use. Optional.
@@ -131,10 +131,10 @@ def plot(self) :
             cvmdesc = self.cvm
 
         if 'title' not in self.meta:
-            title = "%sVs30 Etree Difference Plot For %s" % (location_text, cvmdesc)
+            title = "%sCross Section Difference Plot For %s" % (location_text, cvmdesc)
             self.meta['title'] = title
 
         self.meta['mproperty']="vs"
-        self.meta['difference']="vs30"
+        self.meta['difference']="vs"
 
-        HorizontalSlice.plot(self)
+        CrossSection.plot(self)

pycvm/cross_section.py

@@ -29,7 +29,7 @@
 #  @brief Plots a cross section between two @link common.Point Points @endlink.
 #
 #  Generates a cross section that can either be saved as a file or displayed
-#  to the user. 
+#  to the user or differenced with another plot. 
 class CrossSection:
 
     ##
@@ -86,8 +86,15 @@ def __init__(self, startingpoint, endingpoint, meta={}):
         self.floors = None
         if 'vsfloor' in self.meta and 'vpfloor' in self.meta and 'densityfloor' in self.meta :
             self.floors=self.meta['vsfloor']+","+self.meta['vpfloor']+","+self.meta['densityfloor']
-
-
+
+        if 'scalemin' in self.meta and 'scalemax' in self.meta :
+            ## user supplied a fixed scale bounds
+            self.scalemin=float(self.meta['scalemin'])
+            self.scalemax=float(self.meta['scalemax'])
+        else:
+            self.scalemin=None
+            self.scalemax=None
+
         ## The CVM to use (must be installed with UCVM).
         if 'cvm' in self.meta :
             self.cvm = self.meta['cvm']
@@ -147,7 +154,6 @@ def getplotvals(self, mproperty='vs'):
 #        print("total lat.."+ str(len(depth_list)))
 
         u = UCVM(install_dir=self.installdir, config_file=self.configfile, z_range=self.z_range, floors=self.floors)
-
 ### MEI -- TODO, need to have separate routine that generates cross section datafile
         if (self.datafile != None) :
             ## Private number of x points.
@@ -157,10 +163,14 @@ def getplotvals(self, mproperty='vs'):
             print("\nUsing -->"+self.datafile)
 ##            print("expecting x "+str(self.num_x)+" y "+str(self.num_y))
 
-            if self.datafile.rfind(".raw") != -1:
+
+            if self.datafile.rfind(".binary") != -1 :
                 data = u.import_binary(self.datafile, self.num_x, self.num_y)
-            else:  ## ends in .bin
-                data = u.import_np_float_array(self.datafile, self.num_x, self.num_y)
+            else :
+                if self.datafile.rfind(".raw") != -1 :
+                    data = u.import_raw_data(self.datafile, self.num_x, self.num_y)
+                else:  ## with .bin file
+                    data = u.import_np_float_array(self.datafile, self.num_x, self.num_y)
 
 ## this set of data is only for --datatype: either 'vs', 'vp', 'rho', or 'poisson'
         ## The 2D array of retrieved material properties.
@@ -178,6 +188,8 @@ def getplotvals(self, mproperty='vs'):
                       self.materialproperties[y][x].setProperty('Poisson',tmp)
                     if(mproperty == 'vs'):
                       self.materialproperties[y][x].setProperty('Vs',tmp)
+
+            print("\nUsing --> "+self.datafile) 
         else:
             data = u.query(point_list, self.cvm)
 
@@ -323,33 +335,43 @@ def plot(self) :
         self.min_val=np.nanmin(newdatapoints)
         self.mean_val=np.mean(newdatapoints)
 
-        BOUNDS = u.makebounds()
-        TICKS = u.maketicks()
+        # Set colormap and range
+        colormap = basemap.cm.GMT_seis
+
+        if self.scalemin != None and self.scalemax != None:
+            BOUNDS= u.makebounds(float(self.scalemin), float(self.scalemax), 5)
+            TICKS = u.maketicks(float(self.scalemin), float(self.scalemax), 5)
+            umax=round(self.scalemax)
+            umin=round(self.scalemin)
+            umean=round((umax+umin)/2)
+        else:
+            ## default BOUNDS are from 0 to 5
+            BOUNDS = u.makebounds()
+            TICKS = u.maketicks()
+            umax=round(self.max_val)
+            umin=round(self.min_val)
+            umean=round(self.mean_val)
 
         if mproperty == "vp":
             BOUNDS = [bound * 1.7 for bound in BOUNDS]
             TICKS = [tick * 1.7 for tick in TICKS]
 
-        # Set default colormap and range
-        colormap = basemap.cm.GMT_seis
-        norm = mcolors.BoundaryNorm(BOUNDS, colormap.N)
-
-        umax=round(self.max_val)
-        if( umax < 5 ) :
-            umax=5 
-        umin=round(self.min_val)
-
+##   s, s_r   0,5 / scalemin,scalemax  
+##   sd       umin,umax
+##   b        0,5 / scalemin,scalemax
+##   d, d_r   0,5 / scalemin,scalemax
+##   dd       umin,umax
         if color_scale == "s":
             colormap = basemap.cm.GMT_seis
-            norm = mcolors.Normalize(vmin=0,vmax=umax)
+            norm = mcolors.Normalize(vmin=BOUNDS[0],vmax=BOUNDS[len(BOUNDS) - 1])
         elif color_scale == "s_r":
             colormap = basemap.cm.GMT_seis_r
-            norm = mcolors.Normalize(vmin=0,vmax=umax)
+            norm = mcolors.Normalize(vmin=BOUNDS[0],vmax=BOUNDS[len(BOUNDS) - 1])
         elif color_scale == "sd":
-            BOUNDS= u.makebounds(self.min_val, self.max_val, 5, self.mean_val, substep=5)
             colormap = basemap.cm.GMT_seis
-            TICKS = u.maketicks(self.min_val, self.max_val, 5)
-            norm = mcolors.Normalize(vmin=self.min_val,vmax=self.max_val)
+            BOUNDS= u.makebounds(umin, umax, 5, umean, substep=5)
+            TICKS = u.maketicks(umin, umax, 5)
+            norm = mcolors.Normalize(vmin=BOUNDS[0],vmax=BOUNDS[len(BOUNDS) - 1])
         elif color_scale == "b":
             C = []
             for bound in BOUNDS :
@@ -366,13 +388,17 @@ def plot(self) :
             colormap = pycvm_cmapDiscretize(basemap.cm.GMT_seis_r, len(BOUNDS) - 1)
             norm = mcolors.BoundaryNorm(BOUNDS, colormap.N)  
         elif color_scale == 'dd':
-            BOUNDS= u.makebounds(self.min_val, self.max_val, 5, self.mean_val, substep=5)
-            TICKS = u.maketicks(self.min_val, self.max_val, 5)
+            BOUNDS= u.makebounds(umin, umax, 5, umean, substep=5)
+            TICKS = u.maketicks(umin, umax, 5)
             colormap = pycvm_cmapDiscretize(basemap.cm.GMT_seis, len(BOUNDS) - 1)
             norm = mcolors.BoundaryNorm(BOUNDS, colormap.N)
         else: 
             print("ERROR: unknown option for colorscale.")
-
+
+        if 'difference' in self.meta :
+            bwr = cm.get_cmap('bwr')
+            colormap = pycvm_cmapDiscretize(bwr, len(BOUNDS) - 1)
+
 
 ## MEI, TODO this is a temporary way to generate an output of a cross_section input file
         if( self.datafile == None ):
@@ -412,7 +438,10 @@ def plot(self) :
             if(mproperty.title() == "Density") :
               cbar.set_label(mproperty.title() + " (g/cm^3)")
             else:
-              cbar.set_label(mproperty.title() + " (km/s)")
+              if 'difference' in self.meta :
+                    cbar.set_label(mproperty.title() + " (km)")
+              else:
+                   cbar.set_label(mproperty.title() + " (km/s)")
         else:
             cbar.set_label("Poisson(Vs,Vp)")