Filter PSDs for wind speed

[1]:

from pathlib import Path
import numpy as np

from obspy.core import UTCDateTime as UTC

from data_quality_control import analysis, dqclogging, timelist

# Only for display in documentation!
from IPython.core.display import display, HTML

Input parameters

[2]:

starttime = UTC("2020-12-20")
endtime = UTC("2021-01-12")
#starttime = UTC("2020-01-01")
#endtime = UTC("2021-12-31")

[3]:

# NSLC
nslc_code = "GR.BFO..BHZ"
datadir = Path('output/usage_demo/')

winddatafile = Path("../testing/winddata_bfo.txt")

Set logger

[4]:

logfilename = datadir.joinpath("read_winddata.log")
dqclogging.configure_handlers("INFO", "DEBUG", logfilename)

24-02-09 13:26:04 - data_quality_control - INFO - Find log file at output/usage_demo/read_winddata.log

Get infos on data availability

[5]:

lyza = analysis.Analyzer(datadir, nslc_code,
                            fileunit="year")

24-02-09 13:26:04 - data_quality_control.base.BaseProcessedData - INFO - Start/end time set to: None - None

[6]:

print(lyza)

Analyzer for station GR.BFO..BHZ
Datadir: output/usage_demo
HDF5-file covers 1 year
Filename pattern: output/usage_demo/GR.BFO..BHZ_{year:04d}.hdf5
Loglevel: 10
No data attached.

[7]:

files = lyza.get_available_datafiles()
print(files)

24-02-09 13:26:04 - data_quality_control.analysis.Analyzer - INFO - Looking for pattern output/usage_demo/GR.BFO..BHZ_[0-9][0-9][0-9][0-9].hdf5

['output/usage_demo/GR.BFO..BHZ_2021.hdf5', 'output/usage_demo/GR.BFO..BHZ_2020.hdf5']

[8]:

sdate, edate = lyza.get_available_timerange()
print(sdate, edate)

24-02-09 13:26:04 - data_quality_control.analysis.Analyzer - INFO - Looking for pattern output/usage_demo/GR.BFO..BHZ_[0-9][0-9][0-9][0-9].hdf5
24-02-09 13:26:04 - data_quality_control.base.BaseProcessedData - INFO - Start/end time set to: None - None
24-02-09 13:26:04 - data_quality_control.base.BaseProcessedData - INFO - Reading file output/usage_demo/GR.BFO..BHZ_2020.hdf5
24-02-09 13:26:04 - data_quality_control.base.BaseProcessedData - INFO - Start/end time set to: 2020-12-25T01:00:00.000000Z - 2021-01-01T00:00:00.000000Z
24-02-09 13:26:04 - data_quality_control.base.BaseProcessedData - INFO - Start/end time set to: None - None
24-02-09 13:26:04 - data_quality_control.base.BaseProcessedData - INFO - Reading file output/usage_demo/GR.BFO..BHZ_2021.hdf5
24-02-09 13:26:04 - data_quality_control.base.BaseProcessedData - INFO - Start/end time set to: 2021-01-01T00:00:00.000000Z - 2021-01-09T23:00:00.000000Z

2020-12-25T01:00:00.000000Z 2021-01-09T23:00:00.000000Z

Get wind data

Read wind data

First, we read the raw wind speed data from the text file.

[9]:

#reload(timelist)
awind = timelist.read_winddata(winddatafile)

Interpolate wind data

Wind speed is available about every 6 hours. The power spectral densities are available every 1 hour. We want to filter the PSDs for times when wind speed is in a specific range. To do so, we interpolate the wind speed to the finer grid.

[10]:

timedelta_wind = np.ediff1d(awind[:,0])
print("Report interval of wind speed is approx. {:.0f} h".format(np.mean(timedelta_wind) / 3600))

Report interval of wind speed is approx. 6 h

[11]:

x, f = timelist.read_interp_winddata(winddatafile, starttime, endtime, 3600)

[12]:

timedelta_wind = np.ediff1d(x)
print("New interval of wind speed is {:.1f} h".format(np.mean(timedelta_wind) / 3600))

New interval of wind speed is 1.0 h

Filter wind speed

From the regularly sampled wind data, we can extract those times when windspeed is e.g. \(3 \leq c \leq 6\) m/s.

[13]:

# Extract times with certain wind speed
cmin = 3
cmax = 6
xsel = x[np.logical_and(f >= cmin, f <= cmax)]
tlist = [UTC(xi) for xi in xsel]
print("For {:d} samples of {} is wind speed between {:g} and {:g} m/s".format(
    len(tlist), len(x), cmin, cmax))

For 206 samples of 553 is wind speed between 3 and 6 m/s

[14]:

DATA = lyza.get_data(tlist)
print(lyza.startdate, lyza.enddate)
print(lyza.amplitudes.shape)
print("PSD shape:", lyza.psds.shape)
print("len(tlist):", len(tlist))

fig = lyza.plot_spectrogram()#cmap=plt.cm.gist_heat, vmin=-0.7, vmax=2.4)
#fig_cont.savefig(datadir.joinpath("spectrogram_wind_{:g}-{:g}mps.png"))

24-02-09 13:26:07 - data_quality_control.analysis.Analyzer - INFO - Renewing data
24-02-09 13:26:07 - data_quality_control.analysis.Analyzer - INFO - Extending from file output/usage_demo/GR.BFO..BHZ_2020.hdf5
24-02-09 13:26:07 - data_quality_control.analysis.Analyzer - INFO - Reading file output/usage_demo/GR.BFO..BHZ_2020.hdf5
24-02-09 13:26:07 - data_quality_control.analysis.Analyzer - INFO - Extending from file output/usage_demo/GR.BFO..BHZ_2021.hdf5
24-02-09 13:26:07 - data_quality_control.base.BaseProcessedData - INFO - Start/end time set to: None - None
24-02-09 13:26:07 - data_quality_control.base.BaseProcessedData - INFO - Reading file output/usage_demo/GR.BFO..BHZ_2021.hdf5
24-02-09 13:26:08 - data_quality_control.analysis.Analyzer - INFO - Start/end time set to: 2020-12-25T01:00:00.000000Z - 2021-01-09T23:00:00.000000Z
24-02-09 13:26:08 - data_quality_control.analysis.Analyzer - INFO - Adjusting starttime to available: 2020-12-25T01:00:00.000000Z
24-02-09 13:26:08 - data_quality_control.analysis.Analyzer - INFO - Adjusting endtime to available: 2021-01-09T23:00:00.000000Z

2020-12-25T01:00:00.000000Z 2021-01-09T23:00:00.000000Z
(382,)
PSD shape: (133, 1025)
len(tlist): 206
_images/02_read_winddata_21_2.png 
[15]:

fig_psd = lyza.plot3d_spectrogram()
display(HTML(fig_psd.to_html(include_mathjax="cdn")))

24-02-09 13:26:10 - data_quality_control.analysis.Analyzer - WARNING - Check size of HTML-figure! Your browser might crash!

Note, that the plotting currently uses the time list as axis. Therefore pixels are not equally wide but stretched/squeezed until the next timestamp.

Continuous time range

For comparison, let’s see how the spectrogram looks like when we use the full available time range.

[16]:

DATA = lyza.get_data(starttime, endtime)

24-02-09 13:26:11 - data_quality_control.analysis.Analyzer - INFO - Renewing data
24-02-09 13:26:11 - data_quality_control.analysis.Analyzer - INFO - Extending from file output/usage_demo/GR.BFO..BHZ_2020.hdf5
24-02-09 13:26:11 - data_quality_control.analysis.Analyzer - INFO - Reading file output/usage_demo/GR.BFO..BHZ_2020.hdf5
24-02-09 13:26:11 - data_quality_control.analysis.Analyzer - INFO - Extending from file output/usage_demo/GR.BFO..BHZ_2021.hdf5
24-02-09 13:26:11 - data_quality_control.base.BaseProcessedData - INFO - Start/end time set to: None - None
24-02-09 13:26:11 - data_quality_control.base.BaseProcessedData - INFO - Reading file output/usage_demo/GR.BFO..BHZ_2021.hdf5
24-02-09 13:26:11 - data_quality_control.analysis.Analyzer - INFO - Start/end time set to: 2020-12-25T01:00:00.000000Z - 2021-01-09T23:00:00.000000Z
24-02-09 13:26:11 - data_quality_control.analysis.Analyzer - INFO - Adjusting starttime to available: 2020-12-25T01:00:00.000000Z
24-02-09 13:26:11 - data_quality_control.analysis.Analyzer - INFO - Adjusting endtime to available: 2021-01-09T23:00:00.000000Z

[17]:

print(DATA)

None

[18]:

fig = lyza.plot_spectrogram()
_images/02_read_winddata_27_0.png 
[19]:

fig_psd = lyza.plot3d_spectrogram()
display(HTML(fig_psd.to_html(include_mathjax="cdn")))

24-02-09 13:26:12 - data_quality_control.analysis.Analyzer - WARNING - Check size of HTML-figure! Your browser might crash!

[ ]: