I am using the ObsPy Python library and I want to plot the PPSD but without using the PPSD.plot() function of ObsPy, because I look for using some matplotlib utilities which are not implemented in the obspy.signal.spectral_estimation.PPSD class (see ObsPy doc).

For simplicity, first, I try to create this kind of figure:

ppsd.plot()

The best thing I managed to do is by using the ppsd.current_histogram.T*100/(ppsd.current_histogram_count) to get the probability of ppsd.

Here’s the function I implemented:

from obspy.signal.spectral_estimation import get_nlnm, get_nhnm
from obspy.signal import PPSD
import matplotlib.pyplot as plt
from obspy.imaging.cm import pqlx
import numpy as np
plt.close("all")

file = "CCN_301_DNZ_ppsd.npz" # <= .npz file with the stored ppsd
ppsd = PPSD.load_npz(file)

freq = False

def plot_ppsd_data(ppsd, xaxis_frequency = False, show_percentiles = True, show_mean = True, show_model = True):
    
    """
    
    Manual plot of the ppsd data from obspy library between 0.01 (s or Hz) and 100 (s or Hz)
    The function offers the option to plot the noise models from (Peterson, 1993) 
    and basic statistics.
    
    Parameters
    -----------
    ppsd: obspy class => can be obtained using:
       from obspy.signal import PPSD 
       PPSD.load_npz(file)
    xaxis_frequency: bol, default = False
        if True, shows the plot in frequency domain
    show_percentiles: bol, default = True
        if True, shows the percentiles of the ppsd curves
    show_mean: bol, default = True
        if True, shows the mean of the ppsd curves
    show_model: bol, default = True
        if True, shows the high and low new ambient noise models from (Peterson, 1993) 
    
    Returns
    ----------
    None
    
    """
    
    # Read the start and end date of the PPSD
    start_date = ppsd.times_processed[0].date
    end_date = ppsd.times_processed[-1].date
    start_date_string = start_date.strftime("%Y-%m-%d")
    end_date_string = start_date.strftime("%Y-%m-%d")
    
    # PLOT-----------------------
    plt.figure(figsize = (14,8))
    plt.title(ppsd.station + "." + ppsd.channel + "." + start_date_string + "--" + end_date_string)
    ax = plt.subplot()
    
    # Set logarithmic x axis
    x = np.logspace(np.log10(ppsd.psd_periods[0]), np.log10(ppsd.psd_periods[-1]), ppsd.current_histogram.shape[0])
    y = np.linspace(-200, - 50, 150)
    
    ppsd_2d = ppsd.current_histogram.T*100/(ppsd.current_histogram_count)
    # Statistics
    curves = np.array(ppsd.psd_values)
    percentiles = np.percentile(curves, q=[1, 95], axis=0)
    # Import the NLNM and NHNM
    nlnm_periods, nlnm = get_nlnm()
    nhnm_periods, nhnm = get_nhnm()
    
    X, Y = np.meshgrid(x, y)
    plt.pcolor(X,Y, ppsd_2d,  cmap = pqlx, vmax = 30)
    
    plt.xscale("log")
    plt.grid(True, which='both')
    # Labels
    plt.xlabel("Period [s]")
    plt.ylabel("Amplitude [m2/s4/Hz][dB]")
    cbar = plt.colorbar(ax=ax)
    cbar.set_label("[%]")
    
    
    labels = [None, 0.01, 0.1, 1, 10, 100]
    ax.set_xticklabels(labels)
    plt.xlim(0.01, 500)
    plt.ylim(-200, -50)
    
    # Options------------
    if show_model:
        plt.plot(nlnm_periods, nlnm, c = "grey", linewidth = 3)
        plt.plot(nhnm_periods, nhnm, c = "grey", linewidth = 3)
    
    if show_percentiles:
        plt.plot(x, percentiles[0], c = "k", linewidth = 2)
        plt.plot(x, percentiles[1], c = "k", linewidth = 2)
    if show_mean:
        plt.plot(x, np.mean(curves, axis = 0), c = "k", linewidth = 2)
    
    if xaxis_frequency:
        plt.xlabel("Frequency [Hz]")
        ax.set_xlim(right=100)
        # Flip the x-axis
        ax.invert_xaxis()

        # Maintain the same x tick labels and format
        # Get the current x-axis tick labels and format
        xticklabels = ax.get_xticklabels()[::-1]
        # Set the flipped x-axis tick labels and format
        #ax.set_xticks(xticks)
        ax.set_xticklabels(xticklabels)
        
        # Maintain the same x limits
        xlim = ax.get_xlim()
        ax.set_xlim(xlim)
        
        plt.show()
    

if __name__ == "__main__":
    plot_ppsd_data(ppsd, xaxis_frequency= freq) 

This gives me nearly the same result as ppsd.plot():

my_code_ppsd

But I have the impression that this is a pretty ugly method. Can anyone tell me if there are any better options, please?