Module sprit.sprit_plot

Functions

def parse_plot_string(plot_string)
Expand source code
def parse_plot_string(plot_string):
    """Function to parse a plot string into a list readable by plotting functions

    Parameters
    ----------
    plot_string : str
        Plot string used by sprit.plot_hvsr to define results plot

    Returns
    -------
    list
        A list readable by various sprit plotting functions to show what to include in the results plot.
    """
    plot_list = plot_string.split()

    hvsrList = ['hvsr', 'hv', 'h']
    compList = ['component', 'comp', 'c']
    compPlus = [item + '+' for item in compList]
    specList = ['spectrogram', 'specgram', 'spec','sg', 's']

    hvInd = np.nan
    compInd = np.nan
    specInd = np.nan

    hvIndFound = False
    compIndFound = False
    specIndFound = False

    for i, item in enumerate(plot_list):
        if item.lower() in hvsrList and not hvIndFound:
            # assign the index
            hvInd = i
            hvIndFound = True
        if (item.lower() in compList or item.lower() in compPlus) and not compIndFound:
            # assign the index
            compInd = i
            compIndFound = True
        if item.lower() in specList and not specIndFound:
            # assign the index
            specInd = i
            specIndFound = True

    # Get individual plot lists (should already be correctly ordered)
    if hvInd is np.nan:
        hvsr_plot_list = ['HVSR']

    if compInd is np.nan:
        comp_plot_list = []
        if specInd is np.nan:
            if hvInd is not np.nan:
                hvsr_plot_list = plot_list
            spec_plot_list = []
        else:
            if hvInd is not np.nan:
                hvsr_plot_list = plot_list[hvInd:specInd]
            spec_plot_list = plot_list[specInd:]
    else:
        if hvInd is not np.nan:
            hvsr_plot_list = plot_list[hvInd:compInd]
        
        if specInd is np.nan:
            comp_plot_list = plot_list[compInd:]
            spec_plot_list = []
        else:
            comp_plot_list = plot_list[compInd:specInd]
            spec_plot_list = plot_list[specInd:]

    # Figure out how many subplots there will be
    plot_list_list = [hvsr_plot_list, comp_plot_list, spec_plot_list]

    return plot_list_list

Function to parse a plot string into a list readable by plotting functions

Parameters

plot_string : str
Plot string used by sprit.plot_hvsr to define results plot

Returns

list
A list readable by various sprit plotting functions to show what to include in the results plot.
def plot_cross_section(hvsr_data,
title=None,
fig=None,
ax=None,
use_elevation=True,
show_feet=False,
primary_unit='m',
show_curves=True,
annotate_curves=False,
curve_alignment='peak',
grid_size='auto',
orientation='WE',
interpolation_type='cloughtocher',
interpolate_log_values=True,
surface_elevations=None,
show_peak_points=True,
smooth_bedrock_surface=False,
depth_limit=150,
minimum_elevation=None,
show_bedrock_surface=True,
return_data_batch=True,
show_cross_section=True,
verbose=False,
**kwargs)
Expand source code
def plot_cross_section(hvsr_data,  title=None, fig=None, ax=None, use_elevation=True, show_feet=False, primary_unit='m', 
                       show_curves=True, annotate_curves=False, curve_alignment='peak',
                       grid_size='auto', orientation='WE', 
                       interpolation_type='cloughtocher', interpolate_log_values=True,
                       surface_elevations=None, show_peak_points=True, smooth_bedrock_surface=False,
                       depth_limit=150, minimum_elevation=None, show_bedrock_surface=True,
                       return_data_batch=True, show_cross_section=True, verbose=False,
                       **kwargs):
    """Function to plot a cross section given an HVSRBatch or similar object

    Parameters
    ----------
    hvsr_data : HVSRBatch, list, or similar
        HVSRBatch (intended usage) object with HVSRData objects to show in profile/cross section view
    title : str, optional
        Title to use for plot, by default None
    fig : matplotlib.Figure, optional
        Figure to use for plot, by default None
    ax : matplotlib.Axis, optional
        Axis to use for plot, by default None
    use_elevation : bool, optional
        Whether to use elevation (if True) or depth (if False), by default True
    show_feet : bool, optional
        Whether to show feet (if True) or meters (if False), by default False
    primary_unit : str, optional
        Primary unit to use ('m' or 'ft'), by default 'm'
    show_curves : bool, optional
        Whether to also show curves on plot, by default True
    annotate_curves : bool, optional
        Whether to annotate curves by plotting site names next to them, by default False
    curve_alignment : str, optional
        How to horizontally align the curve. 
        If "peak" the curve will be aligned so that the peak is at the correct latitude or longitude.
        If "max" will align the maximum point of the curve to the correct location.
        If any other value, will align at the surface (i.e., highest frequency). By default 'peak'.
    grid_size : list, optional
        Two item list with height and width of grid for interpolation.
        If "auto" this will be calculated based on the data, by default 'auto'.
    orientation : str, optional
        The orientation of the cross section. 
        Should be either "WE", "EW", "NS", or "SN", by default 'WE'.
    interpolation_type : str, optional
        Interpolation type to use. Uses scipy.interpolation.
        Options include: 'cloughtocher', 'nearest neighbor', 'linear', 
        or 'radial basis function', by default 'cloughtocher'.
    interpolate_log_values : bool, optional
        Whether to use log values of the H/V curve for interpolation (can be useful for better normalizing data)
    surface_elevations : shapely.LineString, optional
        A shapely.LineString object containing surface elevation coordinates along cross section path.
        If None, uses elevation data in HVSRBatch specified by hvsr_data, by default None.
    show_peak_points : bool, optional
        Whether to plot small triangles where peaks were picked, by default True
    smooth_bedrock_surface : bool, optional
        Whether to smooth the bedrock surface when plotting, by default False
    depth_limit : int, optional
        Depth limit for the plot, by default 150
    minimum_elevation : _type_, optional
        Minimum elevation of the plot, by default None
    show_bedrock_surface : bool, optional
        Whether to show the bedrock surface, by default True
    return_data_batch : bool, optional
        Whether to return the HVSRBatch object, by default True
    show_cross_section : bool, optional
        Whether to show the cross section plot, by default True
    verbose : bool, optional
        Whether to print information about the process to terminal, by default False

    Returns
    -------
    figure
        Currently only matplotlib figure supported
    """
    if verbose:
        print("Getting cross section plot configuration")
        
    if fig is None and ax is None:
        fig, ax = plt.subplots()
    elif ax is None and fig is not None:
        fig = fig
        ax = fig.get_axes()[0]
    elif fig is None and ax is not None:
        ax = ax
        fig = plt.figure()
        fig.axes.append(ax)
    else:
        fig = fig
        ax = ax
    plt.sca(ax)
    
    if verbose:
        print("Getting data batch for cross section plot")
    batchExt = None
    if isinstance(hvsr_data, (str, os.PathLike, pathlib.Path)):
        if pathlib.Path(hvsr_data).exists() and pathlib.Path(hvsr_data).is_dir():
            batchExt = 'hvsr'
    hvDataBatch = sprit_hvsr.HVSRBatch(hvsr_data, batch_ext=batchExt)
    
    if verbose:
        print("Sorting and Orienting data")
    # Get orientation/order of data
    nsList = ['ns', "north-south", 'northsouth', 'south', 's']
    snList = ['sn', "south-north", 'southnorth', 'north', 'n']
    weList = ['we', "west-east", 'westeast', 'east', 'e']
    ewList = ['ew', "east-west", 'eastwest', 'west', 'w']

    if str(orientation).lower() in nsList:
        ordercoord = 'latitude'
        order = 'descending'
        profile_direction = 'north-south'
    elif str(orientation).lower() in snList:
        ordercoord = 'latitude'
        order = 'ascending'
        profile_direction  = 'south-north'
    elif str(orientation).lower() in weList:
        ordercoord = 'longitude'
        order = 'ascending'
        profile_direction = 'west-east'
    elif str(orientation).lower() in ewList:
        ordercoord = 'longitude'
        order = 'descending'
        profile_direction = 'east-west'
    else:
        if verbose:
            print(f"Value for orientation={orientation} is not recognized. Using West-East orientation.")
        order = 'ascending'
        ordercoord='longitude'
        profile_direction = 'west-east (default)'

    # Get data in correct order, as specified by orientation parameter
    reverseit = (order == 'descending')
    sorted_sites = sorted(hvDataBatch, key=lambda site: hvDataBatch[site][ordercoord], reverse=reverseit)
    hvDataSorted = [hvDataBatch[h] for h in sorted_sites]

    if verbose:
        print(f'Plotting {len(hvDataBatch.sites)} sites, {profile_direction}.')
        [print(f"\t{hvdata.site[:12]:<12}: {hvdata.longitude:>8.4f}, {hvdata.latitude:>8.4f}, {hvdata.elevation:<6.1f}") for hvdata in hvDataSorted]

    # Get cross section profile
    shapelyPoints = []
    interpData = []
    interpCoords = {'longitude':[], 'latitude':[], 'elevation':[]}
    for i, hvData in enumerate(hvDataSorted):
        if not hasattr(hvData, 'x_elev_m'):
            calc_depth_kwargs = {k: v for k, v in kwargs.items() if k in tuple(inspect.signature(sprit_calibration.calculate_depth).parameters.keys())}
            hvData = sprit_calibration.calculate_depth(hvData, **calc_depth_kwargs, verbose=verbose)
        
        # Create shapely Point objects at each profile location
        x = hvData['longitude']
        y = hvData['latitude']
        z = hvData['elevation']

        shapelyPoints.append(shapely.Point(x, y, z))

        # Points arranged for interpolation
        if interpolate_log_values:
            interpData.extend(list(np.log10(hvData.hvsr_curve)))
        else:
            interpData.extend(list(hvData.hvsr_curve))
        for i, pt in enumerate(hvData.hvsr_curve):
            interpCoords['longitude'].append(x)
            interpCoords['latitude'].append(y)
            interpCoords['elevation'].append(hvData['x_elev_m']['Z'][i])

        # Since already doing loop, ensure hvData has all depth/elev info it needs
        if not hasattr(hvData, 'x_elev_m'):
            calc_depth_kwargs = {k: v for k, v in kwargs.items()
                                      if k in tuple(inspect.signature(sprit_calibration.calculate_depth).parameters.keys())}
            if 'calculate_depth_in_feet' not in calc_depth_kwargs:
                calc_depth_kwargs['calculate_depth_in_feet'] = True
            hvDataSorted[i] = sprit_calibration.calculate_depth(hvData, **calc_depth_kwargs, verbose=verbose)

    xSectionProfile = shapely.LineString(shapelyPoints)
    profileXs, profileYs = xSectionProfile.xy
    
    orderCoordValues = profileXs
    if ordercoord == 'latitude':
        orderCoordValues = profileYs

    minX = min(profileXs)
    minY = min(profileYs)
    maxX = max(profileXs)
    maxY = max(profileYs)

    # Generate grid
    if verbose:
        print("Generating Grid: ", end='')
    xSectionLength = xSectionProfile.length
    if grid_size == 'auto':
        grid_size=(50, 100)

        cellHNumber = grid_size[0]
        cellWNumber = grid_size[1]
    elif isinstance(grid_size, (list, tuple)):
        cellHNumber = grid_size[0]
        cellWNumber = grid_size[1]
    else:
        grid_size=(50, 100)

        cellHNumber = grid_size[0]
        cellWNumber = xSectionLength/grid_size[1]

        if verbose:
            print(f'grid_size value ({grid_size} not recognized, using grid 100 cells wide and 50 cells high: grid_size=(50, 100))')

    cellWSize = xSectionLength/cellWNumber
    
    max_surf_elev = max([hvd.elevation for hvd in hvDataSorted])
    min_br_elev = min([hvd.Table_Report['Peak'][0] for hvd in hvDataSorted])
    elev_range = max_surf_elev - min_br_elev

    max_grid_elev = math.ceil(max_surf_elev)

    # Minimum grid elevation is determined by depth_limit and minimum_elevation
    if str(minimum_elevation).lower() == 'auto':
        min_grid_elev = min_br_elev - (elev_range) * 0.1
    elif isinstance(minimum_elevation, numbers.Number):
        min_grid_elev = minimum_elevation
    elif minimum_elevation is None:
        min_grid_elev = max_grid_elev - depth_limit
    
    xSectionDepth = max_grid_elev - min_grid_elev
    cellHSize = xSectionDepth/cellHNumber

    # Get grid coordinates (all coords in z direction (depth/elev))
    gridZcoords = np.linspace(min_grid_elev, max_grid_elev, cellHNumber)

    gridXDists = np.linspace(0, xSectionProfile.length, cellWNumber)
    gridXcoords = [] # All coords in the "x" direction (along profile)
    for xdist in gridXDists:
        x, y = xSectionProfile.interpolate(xdist).xy
        if 'east' in profile_direction:
            gridXcoords.append(x[0])
        else:
            gridXcoords.append(y[0])
    gridXcoords = np.array(gridXcoords)
    if verbose:
        print(f'Grid generated ({cellWNumber*cellHNumber} cells)\n\tx-range: {xSectionLength:.5f} ({cellWNumber:d} cells, each {cellWSize:.5f} units in size)\n\tz-range: {xSectionDepth:.2f} ({cellHNumber:d} cells, each {cellHSize:.5f} units in size)')

    #print('x', len(interpCoords['longitude']))
    #print('y', len(interpCoords['latitude']))
    #print('z', len(interpCoords['elevation']))
    #print('interp', np.array(interpData).shape)
    if verbose:
        print(f'Beginning interpolation ({interpolation_type})... ', end='')
    
    ctList = ['cloughtocher2dinterpolator', 'cloughtocher', 'ct', 'clough-tocher', 'clough tocher', 'cubic', 'c']
    nearList = ['nearestnd', 'nearest', 'near', 'n']
    linList = ['linearnd', 'linear', 'lin', 'l']
    rbfList = ['radial basis function', 'rbf', 'rbfinterpolator']
    
    if str(interpolation_type).lower() in ctList:
        interp = interpolate.CloughTocher2DInterpolator(list(zip(interpCoords[ordercoord], interpCoords['elevation'])), interpData)
    elif str(interpolation_type).lower() in rbfList:
        interp = interpolate.RBFInterpolator(list(zip(interpCoords[ordercoord], interpCoords['elevation'])), interpData)        
    elif str(interpolation_type).lower() in linList:
        interp = interpolate.LinearNDInterpolator(list(zip(interpCoords[ordercoord], interpCoords['elevation'])), interpData)
    else: # elif str(interpolation_type).lower() in nearList:
        interp = interpolate.NearestNDInterpolator(list(zip(interpCoords[ordercoord], interpCoords['elevation'])), interpData)
        
    xx, zz = np.meshgrid(gridXcoords, gridZcoords)
    interpData = interp(xx, zz)
    interpDataflat = interpData[:-1, :-1]
    if verbose:
        print('Data interpolated')
        print('Plotting colormesh')
    
    # kwargs-defined pcolormesh kwargs
    pcolormeshKwargs = {k: v for k, v in kwargs.items() if k in tuple(inspect.signature(plt.pcolormesh).parameters.keys())}
    
    # Set defaults for cmap and shading (if not overriden in kwargs)
    if 'cmap' not in pcolormeshKwargs:
        pcolormeshKwargs['cmap'] = 'nipy_spectral'
    if 'shading' not in pcolormeshKwargs:
        pcolormeshKwargs['shading'] = 'flat'
                
    ax.pcolormesh(xx, zz, interpDataflat, zorder=0, **pcolormeshKwargs)
    
    if show_curves:
        if verbose:
            print('Plotting curves')
        norm_div = 1
        normal_factor = np.diff(orderCoordValues)
        normal_factor = np.nanmedian(normal_factor[normal_factor != 0]) / norm_div
        
        zAttr = 'x_depth_m'
        if use_elevation:
            zAttr = 'x_elev_m'

        for hvData in hvDataSorted:
            hvData['Normalized_HVCurve'] = (hvData['hvsr_curve'] / np.nanmax(hvData['hvsr_curve'])) * normal_factor
            locatedCurve = hvData['Normalized_HVCurve'] + hvData[ordercoord]
            if curve_alignment.lower() == 'peak':
                normal_peak_factor = (hvData["BestPeak"]['HV']['A0'] / np.nanmax(hvData['hvsr_curve'])) * normal_factor
                locatedCurve = locatedCurve  - normal_peak_factor
            elif curve_alignment.lower() == 'max':
                locatedCurve = locatedCurve  - normal_factor
            else:
                pass
            
            if max(locatedCurve) > max(gridXcoords):
                locatedCurve = locatedCurve - (max(locatedCurve) - max(gridXcoords))
            if min(locatedCurve) < min(gridXcoords):
                locatedCurve = locatedCurve + (min(gridXcoords) - min(locatedCurve))
                
            ax.plot(locatedCurve, hvData[zAttr]['Z'][:-1], c='k', linewidth=0.5, zorder=3)

    if annotate_curves:
        for hvData in hvDataSorted:
            if len(hvData.site) > 10:
                sitename = hvData.site[:8]+ '...'
            else:
                sitename = hvData.site
            ax.text(hvData[ordercoord], y=min_grid_elev, s=sitename, ha='right', va='bottom', rotation='vertical')
    
    if smooth_bedrock_surface:
        show_bedrock_surface = True

    if show_peak_points or show_bedrock_surface:
        brX = []
        brZ = []
        for hvData in hvDataSorted:
            if 'BedrockElevation' in hvData['Table_Report'].columns:
                brX.append(hvData[ordercoord])
                brZ.append(hvData['Table_Report'].loc[0,'BedrockElevation'][()])
        if show_peak_points:
            ax.scatter(brX, brZ, zorder=5, c='k', marker='v')

        
        if smooth_bedrock_surface:
            #brSurfZ = scipy.signal.savgol(brZ, window_length=len(brZ))
            if brX[0] > brX[-1]:
                brX = np.flip(brX)
                brZ = np.flip(brZ)
                doFlip=True
            else:
                doFlip=False
                
            newX = np.sort(gridXcoords)
            brSurfZ = np.interp(newX, brX, brZ)
            brSurfX = newX
        else:
            brSurfX = brX
            brSurfZ = brZ
        
        zMinPts = list(np.array(brSurfZ) * 0 + min(gridZcoords))
        
        if show_bedrock_surface:
            ax.fill_between(brSurfX, brSurfZ, zMinPts,facecolor='w', alpha=0.5, zorder=1)
            ax.plot(brSurfX, brSurfZ, c='k', zorder=2)

    # Plot surfaces
    if verbose:
        print('Plotting surfaces')
    
    if surface_elevations is None:
        surfPts_shapely = []
        surfPtsX = []
        surfPtsZ = []
        
        surface_elevations = shapely.LineString([shapely.Point(hvData['longitude'], 
                                                               hvData["latitude"], 
                                                               hvData["elevation"]) 
                                                 for hvData in hvDataSorted])
    
    xPts = []
    zPts = []
    for surf_pt in surface_elevations.coords:
        surfPtDict = {'longitude':surf_pt[0], 
                      'latitude': surf_pt[1],
                      'elevation': surf_pt[2]}
        xPts.append(surfPtDict[ordercoord])
        zPts.append(surfPtDict['elevation'])
    
    zMaxPts = list(np.array(zPts) * 0 + max_grid_elev)

    # Fill in above surface so interpolation is cleaner and surface topo shape is clear
    ax.fill_between(xPts, zPts, zMaxPts, facecolor='w', zorder=1000)

    # Plot surface topography
    ax.plot(xPts, zPts, c='g', linewidth=1.5, zorder=1001)

    # Plot configuration
    if verbose:
        print('Configuring plot')
    ax.set_xlim([min(gridXcoords), max(gridXcoords)])
    ax.set_ylim([min_grid_elev, max_grid_elev])
    
    ax.tick_params(top=True, labeltop=True, bottom=False, labelbottom=False)
    ax.set_xlabel(str(ordercoord).title())
    ax.xaxis.set_label_position('top')
    ax.set_ylabel('Elevation [Meters]')
    if title is None:
        title = 'HVSR Cross Section Profile'
    ax.set_title(title)
    
    # Display orientation of cross section profile
    # Calculate angle
    profile_angle = math.degrees(math.atan2(shapelyPoints[-1].y - shapelyPoints[0].y, shapelyPoints[-1].x - shapelyPoints[0].x))
    
    # Convert angle to geographic coordinates
    profile_angle = (profile_angle * -1) + 90
    if profile_angle < 0:
        profile_angle += 360

    if verbose:
        print(f"Calculated profile angle to be {profile_angle:.3f} degrees.")
    # Calculate angle name of cross section profile
    if profile_angle < -11.25 + 22.5 * 1:
        profileStart = 'S'
        profileEnd = 'N'
    elif profile_angle < -11.25 + 22.5 * 2:
        profileEnd = 'NNE'
        profileStart = 'SSW'        
    elif profile_angle < -11.25 + 22.5 * 3:
        profileEnd = 'NE'
        profileStart = 'SW'
    elif profile_angle < -11.25 + 22.5 * 4:
        profileEnd = 'ENE'
        profileStart = 'WSW'        
    elif profile_angle < -11.25 + 22.5 * 5:
        profileEnd = 'E'
        profileStart = 'W'        
    elif profile_angle < -11.25 + 22.5 * 6:
        profileEnd = 'ESE'
        profileStart = 'WNW'
    elif profile_angle < -11.25 + 22.5 * 7:
        profileEnd = 'SE'
        profileStart = 'NW'
    elif profile_angle < -11.25 + 22.5 * 8:
        profileEnd = 'SSE'
        profileStart = 'NNW'
    elif profile_angle < -11.25 + 22.5 * 9:
        profileEnd = 'S'
        profileStart = 'N'
    elif profile_angle < -11.25 + 22.5 * 10:
        profileEnd = 'SSW'
        profileStart = 'NNE'        
    elif profile_angle < -11.25 + 22.5 * 11:
        profileEnd = 'SW'
        profileStart = 'NE'
    elif profile_angle < -11.25 + 22.5 * 12:
        profileEnd = 'WSW'
        profileStart = 'ENE'
    elif profile_angle < -11.25 + 22.5 * 13:
        profileEnd = 'W'
        profileStart = 'E'
    elif profile_angle < -11.25 + 22.5 * 14:
        profileEnd = 'WNW'
        profileStart = 'ESE'
    elif profile_angle < -11.25 + 22.5 * 15:
        profileEnd = 'NW'
        profileStart = 'SE'
    elif profile_angle < -11.25 + 22.5 * 16:
        profileEnd = 'NNW'
        profileStart = 'SSE'
    elif profile_angle <= 360:
        profileEnd = 'N'
        profileStart = 'S'

    # Always orient north and east to the right
    if 'north' in profile_direction[:5] or 'east' in profile_direction[:5]:
        ax.invert_xaxis()

    plt.sca(ax)
    plt.figtext(0.1,0.95, s=profileStart)
    plt.figtext(0.9,0.95, s=profileEnd)

    if show_cross_section:
        if verbose:
            print('Displaying plot')
        plt.sca(ax)
        plt.show()
        
    if return_data_batch:
        hvBatch = sprit_hvsr.HVSRBatch(hvDataSorted)
        hvBatch['Cross_Section_Plot'] = fig
        return hvBatch
            
    return fig

Function to plot a cross section given an HVSRBatch or similar object

Parameters

hvsr_data : HVSRBatch, list, or similar
HVSRBatch (intended usage) object with HVSRData objects to show in profile/cross section view
title : str, optional
Title to use for plot, by default None
fig : matplotlib.Figure, optional
Figure to use for plot, by default None
ax : matplotlib.Axis, optional
Axis to use for plot, by default None
use_elevation : bool, optional
Whether to use elevation (if True) or depth (if False), by default True
show_feet : bool, optional
Whether to show feet (if True) or meters (if False), by default False
primary_unit : str, optional
Primary unit to use ('m' or 'ft'), by default 'm'
show_curves : bool, optional
Whether to also show curves on plot, by default True
annotate_curves : bool, optional
Whether to annotate curves by plotting site names next to them, by default False
curve_alignment : str, optional
How to horizontally align the curve. If "peak" the curve will be aligned so that the peak is at the correct latitude or longitude. If "max" will align the maximum point of the curve to the correct location. If any other value, will align at the surface (i.e., highest frequency). By default 'peak'.
grid_size : list, optional
Two item list with height and width of grid for interpolation. If "auto" this will be calculated based on the data, by default 'auto'.
orientation : str, optional
The orientation of the cross section. Should be either "WE", "EW", "NS", or "SN", by default 'WE'.
interpolation_type : str, optional
Interpolation type to use. Uses scipy.interpolation. Options include: 'cloughtocher', 'nearest neighbor', 'linear', or 'radial basis function', by default 'cloughtocher'.
interpolate_log_values : bool, optional
Whether to use log values of the H/V curve for interpolation (can be useful for better normalizing data)
surface_elevations : shapely.LineString, optional
A shapely.LineString object containing surface elevation coordinates along cross section path. If None, uses elevation data in HVSRBatch specified by hvsr_data, by default None.
show_peak_points : bool, optional
Whether to plot small triangles where peaks were picked, by default True
smooth_bedrock_surface : bool, optional
Whether to smooth the bedrock surface when plotting, by default False
depth_limit : int, optional
Depth limit for the plot, by default 150
minimum_elevation : _type_, optional
Minimum elevation of the plot, by default None
show_bedrock_surface : bool, optional
Whether to show the bedrock surface, by default True
return_data_batch : bool, optional
Whether to return the HVSRBatch object, by default True
show_cross_section : bool, optional
Whether to show the cross section plot, by default True
verbose : bool, optional
Whether to print information about the process to terminal, by default False

Returns

figure
Currently only matplotlib figure supported
def plot_depth_curve(hvsr_results,
use_elevation=True,
show_feet=False,
normalize_curve=True,
depth_limit=250,
depth_model=None,
annotate=True,
depth_plot_export_path=None,
fig=None,
ax=None,
show_depth_curve=True)
Expand source code
def plot_depth_curve(hvsr_results, use_elevation=True, show_feet=False, normalize_curve=True, 
                     depth_limit=250, depth_model=None,
                     annotate=True, depth_plot_export_path=None, 
                     fig=None, ax=None, show_depth_curve=True):
    """Function to plot depth curves, given hvsr_results with depth_model specified.

    Parameters
    ----------
    hvsr_results : sprit.HVSRData or sprit.HVSRBatch
        HVSRData object with depth information (or `depth_model` specified).
    use_elevation : bool, optional
        Whether to use elevation (True) or just depth (False), by default True
    show_feet : bool, optional
        Whether to show elevation/depth in feet on Y axis, by default False
    normalize_curve : bool, optional
        Whether to normalize amplitude of H/V curve (x-axis) using maximum/minimum amplitudes, by default True
    depth_limit : int, optional
        Depth limit at which to cut off plot, by default 250 (meters)
    depth_model : None or tuple, optional
        If depth_model not already specified, this can be used to run sprit.calculate_depth() before generating plot.
    annotate : bool, optional
        Whether to annotate plot, by default True
    depth_plot_export_path : filepath-like object, optional
        If specified, will export depth plot to location specifed, by default None
    fig : matplotlib.figure.Figure, optional
        Maplotlib Figure to use for plotting, by default None (new one will be created)
    ax : matplotlib.axis.Axis, optional
        Maplotlib Axis to use for plotting, by default None (new one will be created)
    show_depth_curve : bool, optional
        Whether to diplay the depth curve chart after generating it, by default True

    Returns
    -------
    HVSRData
        HVSRData object with additional property for Depth_Plot
    """


    if fig is None and ax is None:
        fig, ax = plt.subplots(layout='constrained')#, figsize=(5, 15))
        fig.suptitle(hvsr_results['site'])
        ax.set_title('Calibrated Depth to Interface', size='small')
    ax.tick_params(top=True, labeltop=True, bottom=False, labelbottom=False)
    
    if depth_model is not None:
        hvsr_results = sprit_calibration.calculate_depth(hvsr_results, depth_model=depth_model, show_depth_curve=False)

    surfElev = hvsr_results.Table_Report['Elevation'][0]
    bedrockElev = hvsr_results.Table_Report['BedrockElevation'][0]
    bedrockDepth = hvsr_results.Table_Report['BedrockDepth'][0]
    curveRange = max(hvsr_results.hvsr_curve) - min(hvsr_results.hvsr_curve)

    if normalize_curve:
        curvePlot = (hvsr_results.hvsr_curve - min(hvsr_results.hvsr_curve)) / curveRange
        xBase = 0
        xCap = 1
        xLims = [-0.25, 1.25]
        ax.set_xticks([0, 1])
    else:
        curvePlot = hvsr_results.hvsr_curve
        xBase = min(hvsr_results.hvsr_curve)
        xCap = hvsr_results.BestPeak['HV']['A0']#max(hvsr_results.hvsr_curve)
        xLims = [xBase-(0.15*curveRange), xCap+(0.15*curveRange)]
 
    if use_elevation:
        yLims = [hvsr_results.x_elev_m['Z'][0] - depth_limit,
              hvsr_results.x_elev_m['Z'][0]]
        yVals = hvsr_results.x_elev_m['Z'][:-1]
        ax.set_ylabel('Elevation [m]')
        bedrockVal = bedrockElev
        if annotate:
            #Annotate surface elevation
            ax.text(x=xLims[0],
                    y=surfElev,
                    s=str(round(float(surfElev), 2))+'m',
                    ha='right',
                    va='bottom',
                    size='x-small')

            # Annotate bedrock elevation
            ax.text(x=xBase,
                    y=bedrockElev,
                    s=' ' + str(round(float(bedrockElev), 2))+'m\n elevation',
                    ha='left',
                    va='center',
                    size='x-small')
            
            # Annotate bedrock depth
            ax.text(x=xBase,
                    y=max(yLims),
                    s=str(round(float(bedrockDepth), 2))+'m deep ',
                    ha='right',
                    va='top',
                    size='x-small',
                    rotation='vertical')
    else:
        yLims = [depth_limit, hvsr_results.x_depth_m['Z'][0]]
        yVals = hvsr_results.x_depth_m['Z'][:-1]
        ax.set_ylabel('Depth [m]')
        bedrockVal = bedrockDepth
        if annotate:
            # Annotate surface elevation
            ax.text(x=xLims[0],
                    y=0,
                    s=str(round(float(surfElev), 2))+'m',
                    ha='right',
                    va='bottom',
                    size='x-small')
            
            # Annotate Bedrock elevation
            ax.text(x=xBase,
                    y=bedrockVal,
                    s=str(round(float(bedrockElev), 2))+'m\nelevation',
                    ha='center',
                    va='center',
                    size='x-small')

            # Annotate Bedrock depth
            ax.text(x=xBase,
                    y=(min(yLims)+float(bedrockDepth))/2,
                    s=str(round(float(bedrockDepth), 2))+'m deep',
                    ha='right',
                    va='top',
                    size='x-small',
                    rotation='vertical')

    # Plot curve
    ax.fill_betweenx(y=yVals, x1=xBase, x2=curvePlot, alpha=0.2, facecolor='k')
    ax.plot(curvePlot, yVals, c='k', linewidth=0.5)
    if show_feet:
        ax_ft = ax.twinx()
        ax_ft.plot(curvePlot, yVals*3.281, alpha=0)
        ax_ft.set_ylim(yLims[0]*3.281, yLims[1]*3.281)
        ax_ft.set_ylabel('Elevation [ft]')
        if not use_elevation:
            ax_ft.set_ylabel('Depth [ft]')
        
    # Plot peak location
    ax.axhline(y=bedrockVal,
               linestyle='dotted', c='k', linewidth=0.5)
    ax.scatter(xBase, y=bedrockVal, c='k', s=0.5)
    ax.scatter(xCap, y=bedrockVal, c='k', s=0.5)
    
    # Plot "base" line
    ax.axvline(x=xBase, linestyle='dotted', c='k', linewidth=0.5)

    ax.set_ylim(yLims)
    ax.set_xlim(xLims)
    
    xlabel = "H/V Ratio"
    if normalize_curve:
        xlabel += '\n(Normalized)'
        ax.set_xticks([])

    ax.set_xlabel('H/V Ratio')
    ax.xaxis.set_label_position('top')

    ax.set_title(hvsr_results['site'])

    plt.sca(ax)
    fig.set_size_inches(4, 8)
    if show_depth_curve:
        plt.show()
    else:
        plt.close()
        
    if depth_plot_export_path is not None:
        if isinstance(depth_plot_export_path, os.PathLike):
            fig.savefig(depth_plot_export_path)
        else:
            print(f'Please specify a valid path for depth_plot_export_path, not {depth_plot_export_path}')
    
    hvsr_results['Depth_Plot'] = fig
    return hvsr_results

Function to plot depth curves, given hvsr_results with depth_model specified.

Parameters

hvsr_results : HVSRData or HVSRBatch
HVSRData object with depth information (or depth_model specified).
use_elevation : bool, optional
Whether to use elevation (True) or just depth (False), by default True
show_feet : bool, optional
Whether to show elevation/depth in feet on Y axis, by default False
normalize_curve : bool, optional
Whether to normalize amplitude of H/V curve (x-axis) using maximum/minimum amplitudes, by default True
depth_limit : int, optional
Depth limit at which to cut off plot, by default 250 (meters)
depth_model : None or tuple, optional
If depth_model not already specified, this can be used to run sprit.calculate_depth() before generating plot.
annotate : bool, optional
Whether to annotate plot, by default True
depth_plot_export_path : filepath-like object, optional
If specified, will export depth plot to location specifed, by default None
fig : matplotlib.figure.Figure, optional
Maplotlib Figure to use for plotting, by default None (new one will be created)
ax : matplotlib.axis.Axis, optional
Maplotlib Axis to use for plotting, by default None (new one will be created)
show_depth_curve : bool, optional
Whether to diplay the depth curve chart after generating it, by default True

Returns

HVSRData
HVSRData object with additional property for Depth_Plot
def plot_input_stream(hv_data,
stream=None,
input_fig=None,
plot_engine='plotly',
spectrogram_component='Z',
decimate=True,
show_plot=True,
return_fig=False,
**kwargs)
Expand source code
def plot_input_stream(hv_data, stream=None, input_fig=None, plot_engine='plotly',
                        spectrogram_component='Z', decimate=True,
                        show_plot=True, return_fig=False, **kwargs):

    """Function to plot input stream using plotly.
    
    Parameters
    ----------
    hv_data : HVSRData
    stream : obspy.stream.Stream
        Can explictly specify stream instead of using hv_data object.
    input_fig : plotly.Figure
        Plotly figure to plot input stream on. If None, creates new one. Default is None.
    spectrogram_component : str, default='Z'
        Which component to use for the spectrogram
    show_plot : bool, default=True
        Whether to show plot or just generate it.
    return_fig : bool, default=False
        Whether to return figure

    Returns
    -------
    plotly figure
        Only returned if return_fig is True
    """

    plotlyList = ['plotly', 'pltly', 'plty', 'p']
    mplList = ['matplotlib', 'pyplot', 'plt', 'mpl', 'm']

    if str(plot_engine).lower() in plotlyList:
        return _plot_input_stream_plotly(hv_data=hv_data, stream=stream, input_fig=input_fig,
                                         spectrogram_component=spectrogram_component,
                                         show_plot=show_plot, return_fig=return_fig)
    elif str(plot_engine).lower() in mplList:
        return _plot_input_stream_mpl(hv_data, stream, input_fig, decimate=True, spectrogram_component='Z', show_plot=True, return_fig=False, **kwargs)
    else:
        try:
            return _plot_simple_stream_mpl(hv_data, stream, fig=None, axes=None, show_plot=False, ylim_std=0.75, return_fig=True)
        except:
            return _plot_simple_stream_obspy()

Function to plot input stream using plotly.

Parameters

hv_data : HVSRData
 
stream : obspy.stream.Stream
Can explictly specify stream instead of using hv_data object.
input_fig : plotly.Figure
Plotly figure to plot input stream on. If None, creates new one. Default is None.
spectrogram_component : str, default='Z'
Which component to use for the spectrogram
show_plot : bool, default=True
Whether to show plot or just generate it.
return_fig : bool, default=False
Whether to return figure

Returns

plotly figure
Only returned if return_fig is True
def plot_outlier_curves(hvsr_data,
plot_engine='plotly',
plotly_module='go',
remove_outliers_during_plot=False,
outlier_threshold=0.98,
use_percentile=True,
use_hv_curves=False,
from_roc=False,
show_plot=True,
verbose=False,
discarded_curves=None)
Expand source code
def plot_outlier_curves(hvsr_data, plot_engine='plotly', plotly_module='go', remove_outliers_during_plot=False,
                        outlier_threshold=0.98, use_percentile=True, use_hv_curves=False, 
                        from_roc=False, show_plot=True, verbose=False, discarded_curves=None):
    """Function to plot outlier curves, including which have been excluded

    Parameters
    ----------
    hvsr_data : HVSRData
        Input data object
    plot_engine : str = {'plotly', 'matplotlib'}
        Which plotting library to use, by default 'plotly'
    plotly_module : str = {'go', 'px'}
        Which plotly module to use if applicable, by default 'go'
    remove_outliers_during_plot : bool, optional
        Whether curves should also be removed when plotted. During sprit.run(), removal happens separately, so this is False.
    outlier_threshold : float, optional
        RMSE threshold (for removing outliers), by default 0.98
    use_percentile : bool, optional
        Whether to use percentile or raw value, by default True
    use_hv_curves : bool, optional
        Whether to perform analysis on HV curves (if True) or PSD curves (if False), by default False
    from_roc : bool, optional
        Helper parameter to determine if this is being called from sprit.remove_outlier_curves function, by default False
    show_plot : bool, optional
        Whether to show plot, by default True
    verbose : bool, optional
        Whether to print information to terminal, by default False

    Returns
    -------
    plotly figure
        Figure type depends on plotly_module
    """
    orig_args = locals().copy()    
    
    hv_data = hvsr_data
    
    plotlyList = ['plotly', 'plty', 'p']
    mplList = ['matplotlib', 'mpl', 'pyplot', 'mtpltlb', 'm']
    
    if outlier_threshold < 1:
        outlier_threshold = outlier_threshold * 100
    
    roc_kwargs = {'outlier_threshold':outlier_threshold,
                    'use_percentile':True,
                    'use_hv_curves':use_hv_curves,
                    'show_outlier_plot':False,
                    'plot_engine':'None',
                    'verbose':verbose
                    }
    
    if str(plot_engine).lower() in plotlyList:
        #outlier_fig = go.FigureWidget()
        pxList = ['px', 'express', 'exp', 'plotly express', 'plotlyexpress']
        if str(plotly_module).lower() in pxList:
            return __plotly_outlier_curves_px(**orig_args)
        
        outlier_fig = go.Figure()
        
        titleText = 'Outlier Curve Plot'
        if use_hv_curves:
            titleText += ' (H/V Curves)'
        else:
            titleText += ' PSD Curves'
        outlier_fig = go.Figure()
            
        if 'generate_psds_status' in hvsr_data.processing_status.keys() and hvsr_data.processing_status['generate_psds_status']:
            #log_textArea.value += f"\n\n{datetime.datetime.now()}\nremove_outlier_curves():\n'{roc_kwargs}"    
            #hvsr_data = sprit_hvsr.remove_outlier_curves(hvsr_data, **roc_kwargs)
            pass
        else:
            #log_textArea.value += f"\n\n{datetime.datetime.now()}\nremove_outlier_curves() attempted, but not completed. hvsr_data.processing_status['generate_psds_status']=False\n'{roc_kwargs}"
            return outlier_fig

        if roc_kwargs['use_hv_curves']:
            no_subplots = 1
            if hasattr(hvsr_data, 'hvsr_windows_df') and 'HV_Curves' in hvsr_data.hvsr_windows_df.columns:
                outlier_fig.data = []
                outlier_fig.update_layout(grid=None)  # Clear the existing grid layout
                outlier_subp = subplots.make_subplots(rows=no_subplots, cols=1, horizontal_spacing=0.01, vertical_spacing=0.1)
                outlier_fig.update_layout(grid={'rows': 1})
                #outlier_fig = go.FigureWidget(outlier_subp)
                outlier_fig = go.Figure(outlier_subp)

                x_data = hvsr_data['x_freqs']['Z']
                curve_traces = []
                for ind, (i, hv) in enumerate(hvsr_data.hvsr_windows_df.iterrows()):
                    nameLabel = f"Window starting at {i.strftime('%H:%M:%S')}<br>Window #{ind}"
                    curve_traces.append(go.Scatter(x=x_data, y=hv['HV_Curves'], 
                                hovertemplate=nameLabel, line=dict(color='rgba(0,0,0,0.1)', width=0.75),
                                showlegend=False))
                outlier_fig.add_traces(curve_traces)
                
                # Calculate a median curve, and reshape so same size as original
                medCurve = np.nanmedian(np.stack(hvsr_data.hvsr_windows_df['HV_Curves']), axis=0)
                outlier_fig.add_trace(go.Scatter(x=x_data, y=medCurve, line=dict(color='rgba(0,0,0,1)', width=1.5),showlegend=False))
                
                minY = np.nanmin(np.stack(hvsr_data.hvsr_windows_df['HV_Curves']))
                maxY = np.nanmax(np.stack(hvsr_data.hvsr_windows_df['HV_Curves']))
                totalWindows = hvsr_data.hvsr_windows_df.shape[0]
                #medCurveArr = np.tile(medCurve, (curr_data.shape[0], 1))

        else:
            no_subplots = 3
            outlier_fig.data = []
            outlier_fig.update_layout(grid=None)  # Clear the existing grid layout
            outlier_subp = subplots.make_subplots(rows=no_subplots, cols=1, horizontal_spacing=0.01, vertical_spacing=0.02,
                                                    row_heights=[1, 1, 1])
            outlier_fig.update_layout(grid={'rows': 3})
            #outlier_fig = go.FigureWidget(outlier_subp)
            outlier_fig = go.Figure(outlier_subp)

            if hasattr(hvsr_data, 'hvsr_windows_df'):
                rowDict = {'Z':1, 'E':2, 'N':3}
                showTLabelsDict={'Z':False, 'E':False, 'N':True}
                def comp_rgba(comp, a):
                    compstr = ''
                    if comp=='Z':
                        compstr = f'rgba(0, 0, 0, {a})'
                    if comp=='E':
                        compstr = f'rgba(50, 50, 250, {a})'
                    if comp=='N':
                        compstr = f'rgba(250, 50, 50, {a})'
                    return compstr                         
                compNames = ['Z', 'E', 'N']
                rmse_to_plot=[]
                med_traces=[]

                noRemoved = 0
                indRemoved = []
                for i, comp in enumerate(compNames):
                    if hasattr(hvsr_data, 'x_freqs'):
                        x_data = hvsr_data['x_freqs'][comp]
                    else:
                        x_data = [1/p for p in hvsr_data['ppsds'][comp]['period_xedges'][1:]]                    
                    column = 'psd_values_'+comp
                    # Retrieve data from dataframe (use all windows, just in case)
                    curr_data = np.stack(hvsr_data['hvsr_windows_df'][column])
                    
                    # Calculate a median curve, and reshape so same size as original
                    medCurve = np.nanmedian(curr_data, axis=0)
                    medCurveArr = np.tile(medCurve, (curr_data.shape[0], 1))
                    medTrace = go.Scatter(x=x_data, y=medCurve, line=dict(color=comp_rgba(comp, 1), width=1.5), 
                                                    name=f'{comp} Component', showlegend=True)
                    # Calculate RMSE
                    rmse = np.sqrt(((np.subtract(curr_data, medCurveArr)**2).sum(axis=1))/curr_data.shape[1])

                    rmse_threshold = np.percentile(rmse, roc_kwargs['outlier_threshold'])
                    
                    # Retrieve index of those RMSE values that lie outside the threshold
                    timeIndex = hvsr_data['hvsr_windows_df'].index
                    for j, curve in enumerate(curr_data):
                        if rmse[j] > rmse_threshold:
                            badTrace = go.Scatter(x=x_data, y=curve,
                                                line=dict(color=comp_rgba(comp, 1), width=1.5, dash='dash'),
                                                #marker=dict(color=comp_rgba(comp, 1), size=3),
                                                name=str(hvsr_data.hvsr_windows_df.index[j]), showlegend=False)
                            outlier_fig.add_trace(badTrace, row=rowDict[comp], col=1)
                            if j not in indRemoved:
                                indRemoved.append(j)
                            noRemoved += 1
                        else:
                            goodTrace = go.Scatter(x=x_data, y=curve,
                                                    line=dict(color=comp_rgba(comp, 0.01)), name=str(hvsr_data.hvsr_windows_df.index[j]), showlegend=False)
                            outlier_fig.add_trace(goodTrace, row=rowDict[comp], col=1)

                    #timeIndRemoved = pd.DatetimeIndex([timeIndex[ind] for ind in indRemoved])
                    #hvsr_data['hvsr_windows_df'].loc[timeIndRemoved, 'Use'] = False

                    outlier_fig.add_trace(medTrace, row=rowDict[comp], col=1)
                    
                    outlier_fig.update_xaxes(showticklabels=False, row=1, col=1)
                    outlier_fig.update_yaxes(title={'text':'Z'}, row=1, col=1)
                    outlier_fig.update_xaxes(showticklabels=False, row=2, col=1)
                    outlier_fig.update_yaxes(title={'text':'E'}, row=2, col=1)
                    outlier_fig.update_xaxes(showticklabels=True, row=3, col=1)
                    outlier_fig.update_yaxes(title={'text':'N'}, row=3, col=1)

                    outlier_fig.update_layout(margin={"l":10, "r":10, "t":30, 'b':0}, showlegend=True,
                                    title=f"{hvsr_data['site']} Outliers")
                    if comp == 'N':
                        minY = np.nanmin(curr_data)
                        maxY = np.nanmax(curr_data)
                    totalWindows = curr_data.shape[0]
                
                outlier_fig.add_annotation(
                    text=f"{len(indRemoved)}/{totalWindows} outlier windows removed",
                    x=np.log10(max(x_data)) - (np.log10(max(x_data))-np.log10(min(x_data))) * 0.01,
                    y=minY+(maxY-minY)*0.01,
                    xanchor="right", yanchor="bottom",#bgcolor='rgba(256,256,256,0.7)',
                    showarrow=False,row=no_subplots, col=1)


        outlier_fig.update_xaxes(type='log')
        outlier_fig.update_layout(paper_bgcolor='white', plot_bgcolor='white',
                                font_color='black', 
                                title=dict(font_color='black',
                                text=titleText))
        #with outlier_graph_widget:
        #    clear_output(wait=True)
        #    display(outlier_fig)
        if show_plot:
            outlier_fig.show()        
    else: # Matplotlib outlier plot
        #if discarded_curves is not None:
        #    for i, b in enumerate(hvsr_data['hvsr_windows_df'].index[pd.Series(discarded_curves)]):
        #        print('DSCURVES', discarded_curves)
        #        print(b)

        # Determine names of hvsr_windows_df columns to use
        if not use_hv_curves:
            compNames = ['Z', 'E', 'N']
            for col_name in hvsr_data['hvsr_windows_df'].columns:
                if 'psd_values' in col_name and 'RMSE' not in col_name:
                    cName = col_name.split('_')[2]
                    if cName not in compNames:
                        compNames.append(cName)
            col_prefix = 'psd_values_'
            colNames = [col_prefix+cn for cn in compNames]
        else:
            compNames = []
            for col_name in hvsr_data['hvsr_windows_df'].columns:
                if col_name.startswith('HV_Curves') and "Log10" not in col_name:
                    compNames.append(col_name)
            colNames = compNames
            col_prefix = 'HV_Curves'
    
        spMosaic = []
        if use_hv_curves:
            spMosaic.append(['HV Curve'])
            fSize = (8.5, 6)
        else:
            for c in compNames:
                spMosaic.append([c])
            fSize = (8.5, len(compNames) * 2)

            # Intialize to only get unique labels
            rem_label_got = False
            keep_label_got = False

        outlier_fig, ax = plt.subplot_mosaic(spMosaic, sharex=True, figsize=fSize)
            
        # Loop through each component, and determine which curves are outliers
        bad_rmse = []
        for i, column in enumerate(colNames):
            if column in compNames:
                if use_hv_curves == False:
                    column = col_prefix + column
                else:
                    column = column

            # Retrieve data from dataframe (use all windows, just in case)
            curr_data = np.stack(hvsr_data['hvsr_windows_df'][column])
            
            # Calculate a median curve, and reshape so same size as original
            medCurve = np.nanmedian(curr_data, axis=0)
            medCurveArr = np.tile(medCurve, (curr_data.shape[0], 1))

            # Calculate RMSE
            rmse = np.sqrt(((np.subtract(curr_data, medCurveArr)**2).sum(axis=1))/curr_data.shape[1])
            hvsr_data['hvsr_windows_df']['RMSE_'+column] = rmse
            if use_percentile is True:
                rmse_threshold = np.percentile(rmse[~np.isnan(rmse)], outlier_threshold)
                if verbose:
                    print(f'\tRMSE at {outlier_threshold}th percentile for {column} calculated at: {rmse_threshold:.2f}')
            else:
                rmse_threshold = outlier_threshold

            # Retrieve index of those RMSE values that lie outside the threshold
            for j, curve in enumerate(curr_data):
                if rmse[j] > rmse_threshold:
                    bad_rmse.append(j)

            # Iterate through each curve to determine if it's rmse is outside threshold, for plot
            keep_label_got = False
            rem_label_got = False
            for j, curve in enumerate(curr_data):
                label = None
                if rmse[j] > rmse_threshold:
                    linestyle = 'dashed'
                    linecolor='darkred'
                    alpha = 1
                    linewidth = 1
                    if not rem_label_got:
                        label='Removed Curve'
                        rem_label_got=True
                else:
                    linestyle='solid'
                    linecolor = 'rosybrown'
                    alpha = 0.25
                    linewidth=0.5
                    if not keep_label_got:
                        keep_label_got=True
                        label='Retained Curve'

                # Plot each individual curve
                if not use_hv_curves:
                    if 'x_freqs' in hvsr_data.keys():
                        ax[compNames[i]].plot(hvsr_data.x_freqs[compNames[i]], curve, linewidth=linewidth, c=linecolor, linestyle=linestyle, alpha=alpha, label=label)
                    else:
                        ax[compNames[i]].plot(1/hvsr_data.ppsds[compNames[i]]['period_bin_centers'], curve, linewidth=linewidth, c=linecolor, linestyle=linestyle, alpha=alpha, label=label)
                else:
                    if 'x_freqs' in hvsr_data.keys():
                        ax["HV Curve"].plot(hvsr_data.x_freqs['Z'][:-1], curve, linewidth=linewidth, c=linecolor, linestyle=linestyle, alpha=alpha, label=label)
                    else:
                        ax["HV Curve"].plot(1/(hvsr_data.ppsds['Z']['period_bin_centers'][:-1]), curve, linewidth=linewidth, c=linecolor, linestyle=linestyle, alpha=alpha, label=label)                    
            
            # Plot the median curve
            if 'HV_Curves' in compNames[i]:
                axName = 'HV Curve'
                keyName = 'Z'
            else:
                axName = keyName = compNames[i]
            
            if not use_hv_curves:
                if 'x_freqs' in hvsr_data.keys():
                    ax[compNames[i]].plot(hvsr_data.x_freqs[compNames[i]], medCurve, linewidth=1, color='k', label='Median Curve')
                else:
                    ax[compNames[i]].plot(1/hvsr_data.ppsds[compNames[i]]['period_bin_centers'],medCurve, linewidth=1, color='k', label='Median Curve')
            else:
                if 'x_freqs' in hvsr_data.keys():
                    ax['HV Curve'].plot(hvsr_data.x_freqs['Z'][:-1], medCurve, linewidth=1, color='k', label='Median Curve')
                else:
                    ax['HV Curve'].plot(1/hvsr_data.ppsds['Z']['period_bin_centers'][:-1],medCurve, linewidth=1, color='k', label='Median Curve')


            # Format axis
            ax[axName].set_ylabel(f"{compNames[i]}")
            ax[axName].legend(fontsize=10, labelspacing=0.1)
            ax[axName].semilogx()

        outlier_fig.suptitle(f"{hvsr_data['site']}\n Curves Removed from Analysis")
        outlier_fig.set_layout_engine('constrained')
                
        hvsr_data['Outlier_Plot'] = outlier_fig 
    
        if show_plot:
            plt.show()
        else:
            plt.close()
    
    if remove_outliers_during_plot:
        bad_rmse = np.unique(bad_rmse)
        if len(bad_rmse) > 0:
            hvsr_data['hvsr_windows_df']['Use'] = hvsr_data['hvsr_windows_df']['Use'] * (rmse_threshold > hvsr_data['hvsr_windows_df']['RMSE_'+column])

    
    hvsr_data['Outlier_Plot'] = outlier_fig 
    return outlier_fig

Function to plot outlier curves, including which have been excluded

Parameters

hvsr_data : HVSRData
Input data object
plot_engine : str = {'plotly', 'matplotlib'}
Which plotting library to use, by default 'plotly'
plotly_module : str = {'go', 'px'}
Which plotly module to use if applicable, by default 'go'
remove_outliers_during_plot : bool, optional
Whether curves should also be removed when plotted. During sprit.run(), removal happens separately, so this is False.
outlier_threshold : float, optional
RMSE threshold (for removing outliers), by default 0.98
use_percentile : bool, optional
Whether to use percentile or raw value, by default True
use_hv_curves : bool, optional
Whether to perform analysis on HV curves (if True) or PSD curves (if False), by default False
from_roc : bool, optional
Helper parameter to determine if this is being called from sprit.remove_outlier_curves function, by default False
show_plot : bool, optional
Whether to show plot, by default True
verbose : bool, optional
Whether to print information to terminal, by default False

Returns

plotly figure
Figure type depends on plotly_module
def plot_results_plotly(hv_data,
plot_string='HVSR p ann C+ p SPEC ann',
azimuth='HV',
results_fig=None,
results_graph_widget=None,
use_figure_widget=False,
return_fig=False,
show_results_plot=False,
html_plot=False,
verbose=False)
Expand source code
def plot_results_plotly(hv_data, plot_string='HVSR p ann C+ p SPEC ann', azimuth='HV',
                results_fig=None, results_graph_widget=None, use_figure_widget=False,
                return_fig=False, show_results_plot=False, html_plot=False,
                verbose=False,):
    
    """Function to plot results using plotly

    Parameters
    ----------
    hv_data : sprit.HVSRData
        Data object to use for plotting
    plot_string : str, optional
        String for designating what to include in plot(s), by default 'HVSR p ann C+ p SPEC ann'
    results_fig : pyplot figure, optional
        Which pyplot figure to plot data onto. If None, makes new figure, by default None.
    results_graph_widget : plotly graph object widget, optional
        Which pyplot figure to plot data onto. If None, makes new widget, if applicable, by default None.
    use_figure_widget : bool, optional
        Whether to use figure widget, by default False
    return_fig : bool, optional
        Whether to return figure, by default False
    show_results_plot : bool, optional
        Wheather to show plot, by default False
    html_plot : bool, optional
        Whether to create an HTML version of the plot, by default False
    verbose : bool, optional
        Whether to print information to terminal, by default False

    Returns
    -------
    plotly figure
        Only if return_fig is True.
    """

    if results_fig is None:
        results_fig = go.FigureWidget()

    hvsr_data = hv_data

    plotymax = max(hvsr_data.hvsrp2['HV']) + (max(hvsr_data.hvsrp2['HV']) - max(hvsr_data.hvsr_curve))
    if plotymax > hvsr_data.BestPeak['HV']['A0'] * 1.5:
        plotymax = hvsr_data.BestPeak['HV']['A0'] * 1.5
    ylim = [0, plotymax]
    if isinstance(hvsr_data, sprit_hvsr.HVSRBatch):
        hvsr_data = hvsr_data[0]

    hvsrDF = hvsr_data.hvsr_windows_df

    plot_list = parse_plot_string(plot_string)

    combinedComp = False
    # By default there 3 subplots
    noSubplots = 3
    # Remove any subplots that are not indicated by plot_type parameter
    noSubplots = noSubplots - plot_list.count([])
    
    # Now, check if comp plot is combined with HV
    if plot_list[1] != [] and '+' not in plot_list[1][0]:
        combinedComp = True
        noSubplots -= 1
    
    # Get all data for each plotted item
    # Get subplot numbers based on plot_list
    spec = []
    if plot_list[0]==[]:
        # if for some reason, HVSR plot was not indicated, add it
        hv_plot_row = 1 # Default first row to hv (may change later)
        noSubplots += 1
        if plot_list[1] == []:
            comp_plot_row = None
            if plot_list[2] == []:
                spec_plot_row = None
                hv_plot_row = 1 #If nothing specified, do basic h/v plot
            else:
                spec_plot_row = 1 # If only spec specified
        else:
            comp_plot_row = 1 # If no HV specified by comp is, comp is subplot 1

            if plot_list[2] == []:
                spec_plot_row = None
            else:
                spec_plot_row = 2 # If only comp and spec specified comp first then spec
    else:
        hv_plot_row = 1 # HV specified explicitly
        if plot_list[1] == []:
            comp_plot_row = None
            if plot_list[2] == []:
                spec_plot_row = None
            else:
                spec_plot_row = 2 # if no comp specified, spec is 2nd subplot
        else:
            if combinedComp:
                comp_plot_row = 1
                if plot_list[2] == []:
                    spec_plot_row = None
                else:
                    spec_plot_row = 2
            else:
                comp_plot_row = 2
                if plot_list[2] == []:
                    spec_plot_row = None
                else:
                    spec_plot_row = 3       

    specList = []
    rHeights = [1]
    if hv_plot_row == 1:
        if comp_plot_row == 1:
            specList.append([{'secondary_y': True}])
            if spec_plot_row == 2:
                specList.append([{'secondary_y': False}])
    else:
        specList.append([{'secondary_y': False}])

        if noSubplots >= 2:
            specList.append([{'secondary_y': False}])
            rHeights = [1.5,1]
        if noSubplots == 3:
            specList.append([{'secondary_y': False}])
            rHeights = [2,1.5,1]
    
    # Failsafes
    while len(specList)<noSubplots:
        specList.append([{}])

    while len(rHeights)<noSubplots:
        rHeights.append(1)

    # Re-initialize results_fig
    results_fig.data = []
    results_fig.update_layout(grid=None)  # Clear the existing grid layout, in case applicable

    results_fig = make_subplots(rows=noSubplots, cols=1, horizontal_spacing=0.01, vertical_spacing=0.07,
                                specs=specList,
                                row_heights=rHeights)
    results_fig.update_layout(grid={'rows': noSubplots})

    if use_figure_widget:
        results_fig = go.FigureWidget(results_fig)

    if plot_list[1] != []:
        results_fig = _parse_comp_plot_list(hvsr_data, results_fig=results_fig, 
                                           comp_plot_list=plot_list[1])
        results_fig.update_xaxes(title_text='Frequency [Hz]',
                                 row=comp_plot_row, col=1)

    # HVSR Plot (plot this after COMP so it is on top COMP and to prevent deletion with no C+)
    results_fig = _parse_hv_plot_list(hvsr_data, hvsr_plot_list=plot_list, results_fig=results_fig)

    # Will always plot the HV Curve
    results_fig.add_trace(go.Scatter(x=hvsr_data.x_freqs['Z'], y=hvsr_data.hvsr_curve,
                        line={'color':'black', 'width':1.5}, marker=None, name='HVSR Curve'),
                        row=1, col='all')

    # SPEC plot
    if plot_list[2] != []:
        results_fig = _parse_spec_plot_list(hvsr_data, spec_plot_list=plot_list[2], subplot_num=spec_plot_row, azimuth=azimuth, results_fig=results_fig)

    # Final figure updating
    resultsFigWidth = 650

    components_HV_on_same_plot = (plot_list[1]==[] or '+' not in plot_list[1][0])
    if components_HV_on_same_plot:
        compxside = 'bottom'
        secondaryY = True
        showHVTickLabels = True
        showComptickLabels = True
    else:
        compxside = 'bottom'
        secondaryY = False
        showHVTickLabels = True
        showComptickLabels = True
    
    # Update H/V Plot
    results_fig.update_xaxes(type='log', title_text='Frequency [Hz]', title_standoff=0,
                    range=[np.log10(hvsr_data['hvsr_band'][0]), np.log10(hvsr_data['hvsr_band'][1])],
                    side='bottom', showticklabels=showHVTickLabels,
                    row=1, col=1)
    results_fig.update_yaxes(title_text='H/V Ratio', row=1, col=1, 
                             secondary_y=False, range=ylim)

    # Update Component plot
    results_fig.update_xaxes(type='log', overlaying='x', showticklabels=showComptickLabels, title_standoff=0,
                             range=[np.log10(hvsr_data['hvsr_band'][0]), np.log10(hvsr_data['hvsr_band'][1])],
                             side=compxside, row=comp_plot_row, col=1)    
    results_fig.update_yaxes(title_text="PPSD Amp\n[m2/s4/Hz][dB]", secondary_y=secondaryY, row=comp_plot_row, col=1)

    # Update Spec plot
    results_fig.update_yaxes(title_text='H/V Over Time', row=noSubplots, col=1)

    # Update entire figure
    titleString = f"{hvsr_data['site']} Results"
    results_fig.update_layout(margin={"l":40, "r":10, "t":35, 'b':0},
                            showlegend=False, autosize=False, width=resultsFigWidth, height=resultsFigWidth*0.7,
                            title=titleString)
    
    # Reset results_graph_widget and display 
    #if results_graph_widget is not None:
    #    with results_graph_widget:
    #        clear_output(wait=True)
    #        display(results_fig)

    if show_results_plot:
        if html_plot:
            results_fig.write_html(titleString.replace(' ', '_') + 'plot.html', auto_open=True)
        else:
            results_fig.show()
    
    if return_fig:
        return results_fig

Function to plot results using plotly

Parameters

hv_data : HVSRData
Data object to use for plotting
plot_string : str, optional
String for designating what to include in plot(s), by default 'HVSR p ann C+ p SPEC ann'
results_fig : pyplot figure, optional
Which pyplot figure to plot data onto. If None, makes new figure, by default None.
results_graph_widget : plotly graph object widget, optional
Which pyplot figure to plot data onto. If None, makes new widget, if applicable, by default None.
use_figure_widget : bool, optional
Whether to use figure widget, by default False
return_fig : bool, optional
Whether to return figure, by default False
show_results_plot : bool, optional
Wheather to show plot, by default False
html_plot : bool, optional
Whether to create an HTML version of the plot, by default False
verbose : bool, optional
Whether to print information to terminal, by default False

Returns

plotly figure
Only if return_fig is True.