.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/pqpf.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_pqpf.py: PQPF ======================================================== This example demonstrates how to perform operational PQPF up to three hours from raingauge and radar data. .. GENERATED FROM PYTHON SOURCE LINES 10-13 Definitions -------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 13-40 .. code-block:: default import os import numpy as np import pandas as pd from pyresample import utils import xarray import cartopy.feature as cfeature import cartopy.crs as ccrs import matplotlib import matplotlib.pyplot as plt from matplotlib.colors import BoundaryNorm, LinearSegmentedColormap from swirlspy.rad.irisref import read_irisref from swirlspy.qpe.utils import dbz2rr, rr2rf, locate_file, timestamps_ending from swirlspy.qpe.utils import temporal_interp, multiple_acc from swirlspy.obs.rain import Rain from swirlspy.qpf.rover import rover from swirlspy.qpf.sla import sla from swirlspy.qpf.rover import rover from swirlspy.core.resample import grid_resample plt.switch_backend('agg') THIS_DIR = os.getcwd() os.chdir(THIS_DIR) start_time = pd.Timestamp.now() .. rst-class:: sphx-glr-script-out .. code-block:: pytb Traceback (most recent call last): File "/tmp/build/docs/swirlspy/swirlspy/examples/pqpf.py", line 29, in from swirlspy.qpf.rover import rover File "/tmp/build/docs/swirlspy/swirlspy/qpf/rover.py", line 6, in from rover.rover import rover as rover_api ImportError: libopencv_core.so.3.4: cannot open shared object file: No such file or directory .. GENERATED FROM PYTHON SOURCE LINES 41-49 Initialising --------------------------------------------------- This section demonstrates extracting radar reflectivity data. Step 1: Define a basetime. .. GENERATED FROM PYTHON SOURCE LINES 49-53 .. code-block:: default # Supply basetime basetime = pd.Timestamp('201902190800') .. GENERATED FROM PYTHON SOURCE LINES 54-56 Step 2: Using basetime, generate timestamps of desired radar files timestamps_ending() and locate files using locate_file(). .. GENERATED FROM PYTHON SOURCE LINES 56-68 .. code-block:: default # Obtain radar files dir = THIS_DIR + '/../tests/samples/iris/ppi' located_files = [] radar_ts = timestamps_ending( basetime+pd.Timedelta(minutes=6), duration=pd.Timedelta(minutes=60) ) for timestamp in radar_ts: located_files.append(locate_file(dir, timestamp)) .. GENERATED FROM PYTHON SOURCE LINES 69-71 Step 3: Read data from radar files into xarray.DataArray using read_irisref(). .. GENERATED FROM PYTHON SOURCE LINES 71-79 .. code-block:: default reflectivity_list = [] # stores reflec from read_irisref() for filename in located_files: reflec = read_irisref( filename ) reflectivity_list.append(reflec) .. GENERATED FROM PYTHON SOURCE LINES 80-84 Step 4: Define the target grid as a pyresample AreaDefinition. Since the data is in Centered Azimuthal Projection, the source grid also must also be defined as a pyresample AreaDefinition. .. GENERATED FROM PYTHON SOURCE LINES 84-127 .. code-block:: default # Defining target grid area_id = "hk1980_250km" description = ("A 250 m resolution rectangular grid " "centred at HKO and extending to 250 km " "in each direction in HK1980 easting/northing coordinates") proj_id = 'hk1980' projection = ('+proj=tmerc +lat_0=22.31213333333334 ' '+lon_0=114.1785555555556 +k=1 +x_0=836694.05 ' '+y_0=819069.8 +ellps=intl +towgs84=-162.619,-276.959,' '-161.764,0.067753,-2.24365,-1.15883,-1.09425 +units=m ' '+no_defs') x_size = 500 y_size = 500 area_extent = (587000, 569000, 1087000, 1069000) area_def_tgt = utils.get_area_def( area_id, description, proj_id, projection, x_size, y_size, area_extent ) # Defining source grid radius = 256000 sitecoords = reflectivity_list[0].attrs['site'] res = 480 area_id = 'aeqd' description = ("Azimuthal Equidistant Projection " "centered at the radar site " "extending up to {radius:f}m " "in each direction " "with a {res:f}x{res:f} grid resolution ").format( radius=radius, res=res ) proj_id = 'aeqd' projection = ('+proj=aeqd +lon_0={lon:f} ' + '+lat_0={lat:f} +ellps=WGS84 +datum=WGS84 ' + '+units=m +no_defs').format( lon=sitecoords[0], lat=sitecoords[1]) x_size = res y_size = res area_extent = (-radius, -radius, radius, radius) area_def_src = utils.get_area_def( area_id, description, proj_id, projection, x_size, y_size, area_extent ) .. GENERATED FROM PYTHON SOURCE LINES 128-130 Step 5: Reproject the radar data from read_irisref() from Centered Azimuthal (source) projection to HK 1980 (target) projection. .. GENERATED FROM PYTHON SOURCE LINES 130-139 .. code-block:: default reproj_reflectivity_list = [] for reflec in reflectivity_list: reproj_reflec = grid_resample( reflec, area_def_src, area_def_tgt, coord_label=['easting', 'northing'] ) reproj_reflectivity_list.append(reproj_reflec) .. GENERATED FROM PYTHON SOURCE LINES 140-142 Step 6: Assigning reflectivity xarrays at the last three timestamps to variables for use during ROVER QPF. .. GENERATED FROM PYTHON SOURCE LINES 142-149 .. code-block:: default xarray1 = reproj_reflectivity_list[-3] xarray2 = reproj_reflectivity_list[-2] xarray3 = reproj_reflectivity_list[-1] initialising_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 150-159 Running ROVER and Semi-Lagrangian Advection ------------------------------------------- 1. Concatenate required reflectivity xarrays along time dimension. 2. Run ROVER on all members, with the concatenated xarray as the input. 3. Store motion field xarrays as a list. 4. Perform Semi-Lagrangian Advection on all members using the motion fields from rover. 5. Store forecasted reflectivities as list. .. GENERATED FROM PYTHON SOURCE LINES 159-220 .. code-block:: default # Combining two reflectivity DataArrays # the order of the coordinate keys is now ['y', 'x', 'time'] # as opposed to ['time', 'x', 'y'] reflec_concat_6min = xarray.concat([xarray2, xarray3], dim='time') reflec_concat_12min = xarray.concat([xarray1, xarray3], dim='time') # Running rover on 4 members # Mem-1 u1, v1 = rover( reflec_concat_6min, start_level=1, max_level=7, sigma=2.5 ) # Mem-2 u2, v2 = rover( reflec_concat_12min, start_level=2, max_level=7, sigma=2.5 ) # Rover-A u3, v3 = rover( reflec_concat_6min ) # Mem-4 u4, v4 = rover( reflec_concat_12min ) # Storing motion fields for quiver plot motion_list = [[u1, v1], [u2, v2], [u3, v3], [u4, v4]] rover_time = pd.Timestamp.now() # Running SLA on all members z1 = sla( reflec_concat_6min, u1, v1, steps=30 ) z2 = sla( reflec_concat_12min, u2, v2, steps=15 ) z3 = sla( reflec_concat_6min, u3, v3, steps=30 ) z4 = sla( reflec_concat_12min, u4, v4, steps=15 ) # appending all reflectivities to list z_sla_list = [z1, z2, z3, z4] sla_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 221-229 Concatenating observed and forecasted reflectivities --------------------------------------------------- 1. Add forecasted reflectivity to reproj_reflectivity_list. 2. Concatenate observed and forecasted reflectivity xarray.DataArrays along the time dimension. 3. Concatenate reflectivities of different members along a fourth dimension. .. GENERATED FROM PYTHON SOURCE LINES 229-251 .. code-block:: default z_cat_list = [] for reflectivity in z_sla_list: z_all = reproj_reflectivity_list + [reflectivity[1:, ...]] z_cat = xarray.concat(z_all, dim='time') z_cat.attrs['long_name'] = 'reflectivity' z_cat_list.append(z_cat) # Concatenating different members z_ens_6min = xarray.concat( [z_cat_list[0], z_cat_list[2]], xarray.IndexVariable('member', ['Mem-1', 'Mem-3']) ) z_ens_12min = xarray.concat( [z_cat_list[1], z_cat_list[3]], xarray.IndexVariable('member', ['Mem-2', 'Mem-4']) ) concat_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 252-269 Generating radar reflectivity maps ----------------------------------- 1. Define the colour scale and the format of the reflectivity plots. 2. Defining the projection. 3. Defining the zoom or area extent. Tuple order is (x0, x1, y0, y1) as opposed to pyresample (x0, y0, x1, y1). 4. Define figure. 5. Define axes. 6. Set area extent. 7. Plot gridlines and GSHHS coastlines using axes methods. 8. Plot data using xarray.plot(). 9. Plot quiver using axes method. In this example, only hourly images will be plotted. .. GENERATED FROM PYTHON SOURCE LINES 269-354 .. code-block:: default # Defining colour levels levels = [ -32768, 10, 15, 20, 24, 28, 32, 34, 38, 41, 44, 47, 50, 53, 56, 58, 60, 62 ] cmap = matplotlib.colors.ListedColormap([ '#FFFFFF', '#08C5F5', '#0091F3', '#3898FF', '#008243', '#00A433', '#00D100', '#01F508', '#77FF00', '#E0D100', '#FFDC01', '#EEB200', '#F08100', '#F00101', '#E20200', '#B40466', '#ED02F0' ]) norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True) # Plotting maps for forecasts with 6 minute tracking intervals # Defining the crs and coastline type crs = area_def_tgt.to_cartopy_crs() hires = cfeature.GSHHSFeature( scale='h', levels=[1], edgecolor='black', facecolor='none' ) # Define area extent zoom = (712000, 962000, 695000, 945000) # Generating a timelist for every hour timelist = [ (basetime + pd.Timedelta(hours=i)).to_datetime64() for i in range(4) ] for idx, member in enumerate(z_ens_6min.coords['member'].values): # Defining quiver for each member qx = motion_list[idx][0].coords['easting'].values[::5] qy = motion_list[idx][0].coords['northing'].values[::-5] qu = motion_list[idx][0].values[::5, ::5] qv = motion_list[idx][1].values[::5, ::5] for time in timelist: fig = plt.figure(figsize=(20, 15)) ax = plt.axes(projection=crs) ax.set_extent(zoom, crs=crs) ax.add_feature(hires) ax.gridlines() z_ens_6min.sel( member=member, time=time ).plot.pcolormesh(cmap=cmap, norm=norm) ax.quiver(qx, qy, qu, qv, pivot='mid', regrid_shape=20) plt.savefig( THIS_DIR + "/../tests/outputs/rover-output-map-" f"{member}-{pd.Timestamp(time).strftime('%Y%m%d%H%M')}.png" ) # Plotting for members with 12 min tracking interval for idx, member in enumerate(z_ens_12min.coords['member'].values): # Defining quiver for each member qx = motion_list[idx][0].coords['easting'].values[::5] qy = motion_list[idx][0].coords['northing'].values[::-5] qu = motion_list[idx][0].values[::5, ::5] qv = motion_list[idx][1].values[::5, ::5] for time in timelist: fig = plt.figure(figsize=(20, 15)) plt.axis('equal') ax = plt.axes(projection=crs) ax.set_extent(zoom, crs=crs) ax.add_feature(hires) ax.gridlines() z_ens_12min.sel( member=member, time=time ).plot.pcolormesh(cmap=cmap, norm=norm) ax.quiver(qx, qy, qu, qv, pivot='mid', regrid_shape=20) plt.savefig( THIS_DIR + "/../tests/outputs/rover-output-map-" f"{member}-{pd.Timestamp(time).strftime('%Y%m%d%H%M')}.png" ) radar_image_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 355-363 Accumulating hourly rainfall for 3-hour forecast ------------------------------------------------ Hourly accumulated rainfall is calculated every 30 minutes, the first endtime is the basetime i.e. T+0min. | Step 1: Convert reflectivity in dBZ to rainrates in mm/h with dbz2rr(). | Step 2: Convert rainrates to rainfalls in 6 and 12 minutes with rr2rf(). .. GENERATED FROM PYTHON SOURCE LINES 363-379 .. code-block:: default rainrates_ens_6min = dbz2rr(z_ens_6min, a=58.53, b=1.56) rainfalls_ens_6min = rr2rf(rainrates_ens_6min) rainrates_ens_12min = dbz2rr(z_ens_12min, a=58.53, b=1.56) rainfalls_ens_12min = xarray.concat( [rr2rf(rainrates_ens_12min.sel( time=slice(pd.Timestamp('201902190706'), pd.Timestamp('201902190800')), ), scan_duration=6), rr2rf(rainrates_ens_12min.sel( time=slice(pd.Timestamp('201902190806'), pd.Timestamp('201902191100')), ), scan_duration=12)], dim='time' ) .. GENERATED FROM PYTHON SOURCE LINES 380-385 Step 3: Compute hourly accumulated rainfall every 30 minutes. In this example, since the intervals between the time coordinates in rainfalls_ens_12min are uneven, an interpolation is required. .. GENERATED FROM PYTHON SOURCE LINES 385-400 .. code-block:: default acc_rf_6min = multiple_acc( rainfalls_ens_6min, basetime, basetime+pd.Timedelta(hours=3) ) rainfalls_ens_12min_interp = temporal_interp( rainfalls_ens_12min, pd.Timestamp('201902190700'), pd.Timestamp('201902191100') ) acc_rf_12min = multiple_acc( rainfalls_ens_12min_interp, basetime, basetime+pd.Timedelta(hours=3) ) acc_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 401-421 Calculating Probability Exceeding Threshold ---------------------------------------------- In this example, thresholds are 0.5mm, 5mm, 10mm, 30mm, 50mm and 70 mm. Probabilities are 0%, 25%, 50%, 75% and 100%, as there are four members. Steps: 1. Define threshold. 2. Concatenate rainfalls from 6 minute and 12 minute ensembles along the member dimension. 3. Use a loop to calculate the probability of exceeding threshold for each gridcell and store results in list. 4. Concatenate the contents of the list. Result is an xarray with dimensions (threshold, time, y, x). 5. Plot the results using xarray.plot(). In this example, only the 0.5mm, 10mm and 50mm rainfall every hour will be plotted. .. GENERATED FROM PYTHON SOURCE LINES 421-489 .. code-block:: default # Define threshold threshold = [0.5, 5., 10., 30., 50., 70.] # Concatenating rainfalls_ens_6min and rainfalls_ens_12min acc_rf = xarray.concat( [acc_rf_6min, acc_rf_12min], dim='member' ) # Define list to store probabilities of exceeding rainfall thresholds # List corresponds to threshold prob_list = [] # Calculate probability for th in threshold: bool_forecast = acc_rf >= th count = bool_forecast.sum(dim='member') prob = count/len(bool_forecast.coords['member'].values) * 100 prob_list.append(prob) # Generate coordinate xarray for threshold th_xarray = xarray.IndexVariable( 'threshold', threshold ) # concatenate prob_rainfall = xarray.concat( prob_list, dim=th_xarray ) prob_rainfall.attrs['long_name'] = "Probability of Exceeding Threshold" prob_rainfall.attrs['units'] = "%" # Plot the results cmap = LinearSegmentedColormap.from_list( 'custom blue', ['#FFFFFF', '#000099'] ) for time in prob_rainfall.coords['time'].values[2::2]: for th in prob_rainfall.coords['threshold'].values[::2]: plt.figure(figsize=(20, 15)) ax = plt.axes(projection=crs) ax.set_extent(zoom, crs=crs) ax.add_feature(hires) ax.gridlines() lol = prob_rainfall.sel(threshold=th, time=time) prob_rainfall.sel(threshold=th, time=time).plot(cmap=cmap) t_minus = pd.Timestamp(time) - \ basetime-pd.Timedelta(minutes=60) t_minus = t_minus.to_timedelta64().astype('timedelta64[m]') t_plus = pd.Timestamp(time) - basetime t_plus = t_plus.to_timedelta64().astype('timedelta64[m]') plt.title( "Probability of Exceeding Threshold\n" f"Threshold = {th}mm " f"Basetime: {str(basetime)} " f"Start time: t {str(t_minus)} End time: t {str(t_plus)}" ) plt.savefig( THIS_DIR + f"/../tests/outputs/p_{pd.Timestamp(time).strftime('%Y%m%d%H%M')}" f"_threshold_{th}mm.png" ) prob_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 490-506 Rainfall percentiles ---------------------------- 1. Using xarray.DataArray.mean(), calculate the mean rainfall of all gridpoints. 2. Using xarray.DataArray.min(), find the minimum rainfall of all gridpoints. 3. Using xarray.DataArray.max(), find the maximum rainfall of all gridpoints. 4. Using xarray.DataArray.quantile() find the 25th, 50th and 75th percentile rainfall of all gridpoints. 5. Concatenate rainfall along percentile dimension. 6. Plot results using xarray.plot(). In this example only the minimum, 50th and 75th percentile rainfall every hour will be plotted. .. GENERATED FROM PYTHON SOURCE LINES 506-593 .. code-block:: default # mean mean_rainfall = acc_rf.mean(dim='member', keep_attrs=True) # max max_rainfall = acc_rf.max(dim='member', keep_attrs=True) # min min_rainfall = acc_rf.min(dim='member', keep_attrs=True) # quartiles q_rainfall = acc_rf.quantile( [.25, .5, .75], dim='member', interpolation='nearest', keep_attrs=True) # generate index percentile = xarray.IndexVariable( 'percentile', ['Minimum', '25th Percentile', '50th Percentile', 'Mean', '75th Percentile', 'Maximum'] ) # concatenating p_rainfall = xarray.concat( [min_rainfall, q_rainfall.sel(quantile=.25).squeeze().drop('quantile'), q_rainfall.sel(quantile=.50).squeeze().drop('quantile'), mean_rainfall, q_rainfall.sel(quantile=.75).squeeze().drop('quantile'), max_rainfall], dim=percentile ) p_rainfall.attrs['long_name'] = "Hourly Accumulated Rainfall" p_rainfall.attrs['units'] = "mm" # Defining levels # Defining the colour scheme levels = [ 0, 0.5, 2, 5, 10, 20, 30, 40, 50, 70, 100, 150, 200, 300, 400, 500, 600, 700 ] cmap = matplotlib.colors.ListedColormap([ '#ffffff', '#9bf7f7', '#00ffff', '#00d5cc', '#00bd3d', '#2fd646', '#9de843', '#ffdd41', '#ffac33', '#ff621e', '#d23211', '#9d0063', '#e300ae', '#ff00ce', '#ff57da', '#ff8de6', '#ffe4fd' ]) norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True) # Plotting for time in p_rainfall.coords['time'].values[2::2]: for pos in p_rainfall.coords['percentile'].values[::2]: fig = plt.figure(figsize=(20, 15)) ax = plt.axes(projection=crs) ax.add_feature(hires) ax.set_extent(zoom, crs=crs) ax.gridlines() p_rainfall.sel(percentile=pos, time=time).plot( cmap=cmap, norm=norm ) t_minus = pd.Timestamp(time) - \ basetime-pd.Timedelta(minutes=60) t_minus = t_minus.to_timedelta64().astype('timedelta64[m]') t_plus = pd.Timestamp(time) - basetime t_plus = t_plus.to_timedelta64().astype('timedelta64[m]') plt.title( f"{pos} Rainfall Intensity\n" f"Basetime: {str(basetime)}\n" f"Start time: {str(t_minus)} " f"End time: {str(t_plus)}" ) position = pos.split(" ")[0] plt.savefig( THIS_DIR + "/../tests/outputs/rainfall_" f"{pd.Timestamp(time).strftime('%Y%m%d%H%M')}" f"_{position}.png" ) ptime = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 594-606 Extract the rainfall values at a specified location ------------------------------------------------------------------ In this example, the rainfall values at the location is assumed to be the same as the nearest gridpoint. 1. Read information regarding the radar stations into a pandas.DataFrame. 2. Extract the rainfall values at the nearest gridpoint to location for given timesteps (in this example, 30 minute intervals). 3. Store rainfall values over time in an xarray.DataArray. 4. Plot the time series of rainfall with boxplots at desired station. In this case, the 15th percentile member is plotted. .. GENERATED FROM PYTHON SOURCE LINES 606-673 .. code-block:: default # Getting radar station coordinates df = pd.read_csv( os.path.join(THIS_DIR, "../tests/samples/hk_raingauge.csv"), usecols=[0, 1, 2, 3, 4] ) # Extract rainfall values at gridpoint closest to the # location specified for given timesteps and storing it # in xarray.DataArray. station_rf_list = [] station_name = [] for index, row in df.iterrows(): station_rf_list.append(p_rainfall.sel( northing=row[1], easting=row[2], method='nearest' ).drop('northing').drop('easting')) station_name.append(row[0]) station_name_index = xarray.IndexVariable( 'ID', station_name ) station_rf = xarray.concat( station_rf_list, dim=station_name_index ) # Extracting the 15th ranked station xr_15_percentile = station_rf.quantile( .15, dim='ID', interpolation='nearest').drop('quantile') # Plotting _, tax = plt.subplots(figsize=(20, 15)) plt.plot( np.arange(1, len(xr_15_percentile.coords['time'].values) + 1), xr_15_percentile.loc['Mean'].values, 'ko-' ) # plot line # Storing percentiles as dictionary to call ax.bxp for boxplot stats_list = [] for i in range(len(xr_15_percentile.coords['time'].values)): stats = { 'med': xr_15_percentile.loc['50th Percentile'].values[i], 'q1': xr_15_percentile.loc['25th Percentile'].values[i], 'q3': xr_15_percentile.loc['75th Percentile'].values[i], 'whislo': xr_15_percentile.loc['Maximum'].values[i], 'whishi': xr_15_percentile.loc['Minimum'].values[i] } stats_list.append(stats) # Plot boxplot tax.bxp( stats_list, showfliers=False ) # Labels xcoords = xr_15_percentile.coords['time'].values xticklabels = [pd.to_datetime(str(t)).strftime("%-H:%M") for t in xcoords] tax.set_xticklabels(xticklabels) plt.title('Time Series of Hourly Accumulated Rainfall') plt.ylabel("Hourly Accumulated Rainfall [mm]") plt.xlabel("Time") plt.savefig(THIS_DIR+"/../tests/outputs/pqpf_time_series.png") extract_time = pd.Timestamp.now() .. GENERATED FROM PYTHON SOURCE LINES 674-677 Checking run time of each component -------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 677-709 .. code-block:: default print(f"Start time: {start_time}") print(f"Initialising time: {initialising_time}") print(f"Rover time: {rover_time}") print(f"SLA time: {sla_time}") print(f"Concatenating time: {concat_time}") print(f"Plotting radar image time: {radar_image_time}") print(f"Accumulating rainfall time: {acc_time}") print( "Calculate and plot probability exceeding threshold: " f"{prob_time}" ) print( f"Plotting rainfall maps: {ptime}" ) print(f"Extracting and plotting time series time: {extract_time}") print(f"Time to initialise: {initialising_time-start_time}") print(f"Time to run rover: {rover_time-initialising_time}") print(f"Time to perform SLA: {sla_time-rover_time}") print(f"Time to concatenate xarrays: {concat_time - sla_time}") print(f"Time to plot radar image: {radar_image_time - concat_time}") print(f"Time to accumulate rainfall: {acc_time - radar_image_time}") print( "Time to calculate and plot probability exceeding threshold: " f"{prob_time-acc_time}" ) print(f"Time to plot rainfall maps: {ptime-prob_time}") print( f"Time to extract station data and plot time series: " f"{extract_time-ptime}" ) .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.003 seconds) .. _sphx_glr_download_auto_examples_pqpf.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: pqpf.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: pqpf.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_