small improvements to psd analysis, added plot_slopes methods, extrac…

…t profile adapted to return object.

source/pyProfile/psd.py

@@ -5,7 +5,7 @@
 #from pySurf.data2D import projection
 
 def psd2prof(f,p,phase=None,N=None):
-    """build a profile from PSD, if phase (in rad) is not passed use random 
+    """build a signal from PSD, if phase (in rad) is not passed use random 
     values. A real profile of npoints is assumed. Note that profiles with 
     (len(p)-1)*2 and (len(p)-1)*2+1 points give psd with same number of points,
     so both reconstructions are possible. Odd number of points is the default,
@@ -45,10 +45,10 @@ def normPSD(N,L=None,form=1):
     """
 
     if form==0:
-        factor=1. #units: [Y**2] no normalization. "sum squared amplitude" in NR 13.4.1
+        factor=1. #units: [Y**2] no normalization. "sum squared amplitude" in NR 13.4.1. rms is the sum. 
     if form==1:
-        factor=1./N**2*L #[Y**2][X] This has the advantage that measurements on same range with different number of points match
-        #and it is the way it is usually plotted. The rms is the integral of PSD over frequency rather than the sum.
+        factor=1./N**2*L #[Y**2][X] This has the advantage that measurements on same range with different number of points match, since the PSD is actually normalized by the frequency interval.
+        #This is the default and the way it is usually plotted. The rms is the integral of PSD over frequency rather than the sum.
     elif form==2:
         factor=1./N**2 #[Y**2] here rms is the sum of PSD, this is RA normalization (however it differs in RA when a window is used).
         #"" 13.4.5 in NR
@@ -287,7 +287,7 @@ def plot_psd(f,p,units=None,label=None,span=0,psdrange=None,
     if units is None: units=['[X]','[Y]','[Z]']
     if len(units)==2: units = [None,units[0],units[1]]
     #print ("FIX THIS ROUTINE BEFORE USING IT, see newview_plotter")
-    plt.ylabel('axial PSD ('+units[2]+'$^2$ '+units[1]+')')
+    plt.ylabel('PSD ('+units[2]+'$^2$ '+units[1]+')')
     plt.xlabel('Freq. ('+units[1]+'$^{-1}$)')
     #pdb.set_trace()
     if len(p.shape)==2: #experimentally deals with 2D array automatically plotting average

source/pySurf/data2D.py

@@ -240,11 +240,11 @@ def _validate_degree(degree, axis):
             raise ValueError("Degree must be scalar when leveling by line (axis != None): received {}".format(degree))
         return degree
 
-    def _level_by_line(data, degree):
+    def _level_by_line(data, degree, *args, **kwargs):
         leg = fitlegendre(y, data, degree, *args, **kwargs)
         return leg
 
-    def _level_plane(data, degree):
+    def _level_plane(data, degree, *args, **kwargs):
         xo, yo = degree if np.size(degree) == 2 else (degree, degree)
         xl, yl = [f.flatten() for f in np.meshgrid(np.arange(xo + 2), np.arange(yo + 1))]
         leg = legendre2d(data, x, y, xl=xl, yl=yl, *args, **kwargs)[0]
@@ -258,11 +258,11 @@ def _level_plane(data, degree):
     degree = _validate_degree(degree, axis)
 
     if axis == 0:
-        leg = _level_by_line(data, degree)
+        leg = _level_by_line(data, degree, *args, **kwargs)
     elif axis == 1:
-        leg = _level_by_line(data.T, degree).T
+        leg = _level_by_line(data.T, degree, *args, **kwargs).T
     elif axis is None:
-        leg = _level_plane(data, degree)
+        leg = _level_plane(data, degree, *args, **kwargs)
     else:
         raise ValueError("Invalid axis value: {}. Axis must be None, 0, or 1.".format(axis))
 
@@ -1686,6 +1686,7 @@ def plot_slope_slice(wdata,x,y,scale=(1.,1.,1.),vrange=None,srange=None,filter=F
     data outside of the range are also excluded from rms calculation.
     """
 
+    # TODO: plot of y axis slope (ax3) is not aligned with surface plot.
     #Note very important that filtering is made pointwise on the 2D slope figure:
     #  all values out of range are masked and all lines containing at least one value are
     #  excluded.

source/pySurf/data2D_class.py

@@ -16,6 +16,8 @@
                            read_data, register_data, resample_data,
                            rotate_data, save_data, slope_2D, subtract_data,
                            sum_data, transpose_data)
+from .data2D import plot_slope_slice, plot_slope_2D
+from .psd2d import avgpsd2d
 from pySurf.points import matrix_to_points2, points_autoresample
 from pySurf.psd2d import (plot_psd2d, plot_rms_power, psd2d, psd2d_analysis,
                           rms_power)
@@ -109,24 +111,35 @@ def __array_finalize__(self, obj):
 
 class Data2D(object):  # np.ndarrays
     """A class containing 2D data with x and y coordinates and methods for analysis.
-    """ 
     
-    # Args:
-    #     object ([type]): [description]
+        Can be initialized with data | data, x, y | file | file, x, y.
+        if x and y are coordinates if match number of elements, 
+            or used as range if two element. 
+        If provided together with file, they override x and y 
+            read from file.
+        
+        Function methods return a copy with new values and don't alter original
+           object. Reassign to variable to use as modifier:
+           e.g. a=a.level( ... )
 
-    # Raises:
-    #     ValueError: [description]
-    #     ValueError: [description]
-    #     ValueError: [description]
-    #     ValueError: [description]
-    #     ValueError: [description]
-    #     ValueError: [description]
+        Args:
+            data (2D array or string): 2D data or file name (suitable reader must 
+            provided).
+            
+            x, y (array): coordinates or ranges.
+            
+            file (str): alternative way to provide a data file.
+            
+            units (str array): 3-element array with units symbols for `x`, `y`, `data`  
+            
+            reader (function): reader function (see `pySurf.readers.instrumentReader`).
+            
+            name (str): sets the name of created object.  
 
-    # Returns:
-    #     [type]: [description]
-
+            *args, **kwargs: optional arguments for `pySurf.data2D.register_data`. 
+                        
+        """
 
-    """
 
     def __new__(subtype, shape, dtype=float, buffer=None, offset=0,
                 strides=None, order=None, info=None):
@@ -188,33 +201,7 @@ def __init__(
         *args,
         **kwargs
     ):
-        """
-        Can be initialized with data | data, x, y | file | file, x, y.
-        if x and y are coordinates if match number of elements, 
-            or used as range if two element. 
-        If provided together with file, they override x and y 
-            read from file.
-        
-        Function methods return a copy with new values and don't alter original
-           object. Reassign to variable to use as modifier:
-           e.g. a=a.level( ... )
-
-        Args:
-            data (2D array or string): 2D data or file name (suitable reader must 
-            provided).
-            
-            x, y (array): coordinates or ranges.
-            
-            file (str): alternative way to provide a data file.
-            
-            units (str array): 3-element array with units symbols for `x`, `y`, `data`  
-            
-            reader (function): reader function (see `pySurf.readers.instrumentReader`).
-            
-            name (str): sets the name of created object.  
 
-            *args, **kwargs: optional arguments for `pySurf.data2D.register_data`. 
-        """
 
         # from pySurf.readers.instrumentReader import reader_dic
 
@@ -696,10 +683,15 @@ def remove_outliers(self, fill_value=np.nan, mask=False, *args, **kwargs):
 
     remove_outliers = update_docstring(remove_outliers, outliers.remove_outliers)
 
-    def extract_profile(self, *args, **kwargs):
+    def extract_profile(self, raw=False, *args, **kwargs):
+        """ Extract one or more profiles from start to end points. Return a `profile_class.Profile` object unless `raw` is True."""
+        
         p = self.topoints()
         prof = points.extract_profile(p,  *args, **kwargs)
-        return prof
+        if raw:
+            return prof
+        else:
+            return Profile(*prof,units=[self.units[0],self.units[2]])
 
     extract_profile = update_docstring(extract_profile, points.extract_profile)
     
@@ -753,7 +745,30 @@ def slope(self, *args, **kwargs):
 
     slope = update_docstring(slope, slope_2D)
 
+    def plot_slope(self, slice = False, scale = (1.,1.,1.) , *args, **kwargs):
+        """Use `data2D.plot_slope_2D` and `.plot_slope_slice` to ploto 4-panel x and y slopes.
+
+        Plot surface, x and y slope maps (rms profile if `slice` is set) and slope (or rms) distribution.
+        Accept all keywords for `data2D.plot_slope_2D` and `.plot_slope_slice`."""
+        
+        # check for scale
+        if self.units is not None:
+            if len(set(self.units)) != 1 and len(set(scale)) == 1:
+                # more than one unit, but same scaling
+                print ("WARNING: units are different for axis, but scaling is uniform")
+                print ("units: ", self.units)
+                print ("scale: ", scale)
+
+        if slice:
+            plot_slope_slice(self.data, self.x, self.y, *args, **kwargs)
+        else:
+            plot_slope_2D(self.data, self.x, self.y, *args, **kwargs)
+
+    plot_slope = update_docstring(plot_slope, plot_slope_2D)
+    plot_slope = update_docstring(plot_slope, plot_slope_slice)
 
+    
+    
 class PSD2D(Data2D):
     """It is a type of data 2D with customized behavoiur and additional properties
     and methods.
@@ -780,7 +795,7 @@ def avgpsd(self, *args, **kwargs):
         from pyProfile.profile_class import PSD
         
         return PSD(res.x, res.y, units = res.units, name = res.name, *args, **kwargs)
-        
+    avgpsd = update_docstring(avgpsd, avgpsd2d)
     
     def rms_power(self, plot=False, rmsrange=None, *args, **kwargs):
         """Calculate rms slice power by integrating .

source/pySurf/psd2d.py

@@ -148,12 +148,12 @@ def plot_rms_power(f,p,x=None,rmsrange=None,ax2f=None,units=None,*args,**kwargs)
     plt.plot(x,rms,*args,**kwargs)
 
     ax3=plt.gca()
-    tit1,tit2=(['Left y axis','Right y axis'] if np.any(ax2f) else [None,None]) #legend title headers
+    tit1,tit2=(['Left y axis','Right y axis'] if np.any(ax2f) else [None,None]) #legend title headers for multiple y axes if needed.
     l1=ax3.legend(loc=loc1,title=tit1)
 
     #plt.title('Total rms power=%6.3g'%(np.sqrt((rms**2).sum()))+((" "+units[2]) if units[2] is not None else ""))  #wrong math, forgets average?
     plt.title('Total rms power=%6.3g'%(np.sqrt(np.nansum(rms**2)/(np.sum(~np.isnan(rms)))))+((" "+units[2]) if units[2] is not None else ""))
-    c = plt.rcParams['axes.prop_cycle'] # c=plt.gca()._get_lines.prop_cycler
+    c=plt.gca()._get_lines.prop_cycler  # c = plt.rcParams['axes.prop_cycle'] # 2024/01/27 not working any more
 
     rms_v=rms
     if rmsrange is not None:

source/pySurf/scripts/MIRO_analysis_functions.py

@@ -488,7 +488,7 @@ def trim_datadic(ranges,datadic_val,xrange=None,yrange=None,
 '''
 
 def plot_all_psds(pathlabels,outfolder=None,subfix='_psd',xrange=None,yrange=None):
-    """passing a datadic plots all psds for each key.
+    """passing a datadic plots all psds for each key. WARNING: units are hardcoded. 
     
     Generate a formatted plot (N.B.: units are fixed) of all PSDs for zipping pathlabel in a single graph 
     and prints their rms.
@@ -575,7 +575,7 @@ def psds_table(tdic,outfolder=None,subfix='_psdtable',ranges = None):
             print('attenzione')
             fs=fs[1:]
             ps=ps[1:]
-        rms_full = np.sqrt(np.trapz(ps,fs))
+        rms_full = np.sqrt(np.trapz(ps,fs)) 
         print (os.path.basename(name),' : ','full freq. range [%6.3f:%6.3f], rms %6.3f'%(*(span(fs)),rms_full))
         rmss.append([span(fs),rms_full])
         #pdb.set_trace()