peak-period-estimation/observed_data.py at main · sofarocean/peak-period-estimation · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
"""
Contents: Routines to load data from observed data files. Also so convinience routines that save intermediate results
to disk so that plotting can be done faster once results are calculated.

These files serve as the companion to the manuscript:

    "Continuous peak-period estimates from discrete surface-wave spectra", Smit et al., 2023

For details we on methods and results, please refer to the manuscript.

Copyright (C) 2023
Sofar Ocean Technologies

Authors: Pieter Bart Smit
"""
from roguewavespectrum import FrequencySpectrum
from scipy.signal.windows import get_window
import xarray
import os
import numpy
from xarray import DataArray
import pickle

def get_data():
    kinds = ['reference','target','monotone','natural']
    spectra = {}
    tp = {}
    hm0 = {}
    for kind in kinds:
        spectra[kind] = get_spectrum(kind)
        tp[kind] = get_peak_period(kind)
        hm0[kind] = spectra[kind].significant_waveheight

    return spectra, tp, hm0

def get_spectrum(kind) -> FrequencySpectrum:
    """
    Get spectra from disk or calculate it if it does not exist yet. The returned spectrum takes the form
    of a roguewave.FrequencySpectrum object.

    :param kind: one of 'reference', 'target', 'monotone', 'natural'
    :return: spectrum
    """


    if kind == 'reference':
        kwargs = {'segment_length_seconds': 3600, 'use_u':True,'use_v':True}
        name = f'./data/spectrum_reference.nc'


        return FrequencySpectrum.from_netcdf(name) # type: FrequencySpectrum


    elif kind == 'target':
        kwargs = {'window':get_window('hann', 2048),'spectral_window':numpy.ones(9),
                  'segment_length_seconds':3600, 'use_u':True,'use_v':True}
        name = f'./data/spectrum_target.nc'

        return FrequencySpectrum.from_netcdf(name)
    elif kind == 'monotone':
        name = './data/spectrum_monotone.nc'

        if os.path.exists(name):
            with open(name, "rb") as file_handle:
                spec = pickle.load(file_handle)
            return spec

        else:
            native = get_spectrum('reference')
            hr = get_spectrum('target')
            spec = native.interpolate_frequency(hr.frequency, method='spline')
            with open(name, "wb") as file_handle:
                pickle.dump(spec, file_handle)
            return spec

    elif kind == 'natural':
        name = './data/spectrum_natural.nc'

        if os.path.exists(name):
            if os.path.exists(name):
                with open(name, "rb") as file_handle:
                    spec = pickle.load(file_handle)
                return spec

        else:
            native = get_spectrum('reference')
            hr = get_spectrum('target')
            spec = native.interpolate_frequency(hr.frequency, method='spline', monotone_interpolation=False)
            with open(name, "wb") as file_handle:
                pickle.dump(spec, file_handle)
            return spec

    else:
        raise Exception(f'unknown kind {kind}')

def get_peak_period(kind) -> DataArray:
    """
    Get peak period from disk or calculate it if it does not exist yet. The returned peak period takes the form
    of a xarray.DataArray object.

    :param kind: one of 'reference', 'target', 'monotone', 'natural'
    :return: peak periods
    """

    if kind == 'reference':
        return get_spectrum(kind).peak_period()

    elif kind == 'target':
        return get_spectrum(kind).peak_period()

    elif kind == 'monotone':
        name = './data/tp_monotone.nc'

        if os.path.exists(name):
            with open(name, "rb") as file_handle:
                tp = pickle.load(file_handle)
            return tp(name)  # type: DataArray
        else:
            spec = get_spectrum('reference')
            tp = spec.peak_period(use_spline=True)
            with open(name, "wb") as file_handle:
                pickle.dump(tp, file_handle)
            return tp

    elif kind == 'natural':
        name = './data/tp_natural.nc'

        if os.path.exists(name):
            with open(name, "rb") as file_handle:
                tp = pickle.load(file_handle)
            return tp(name)  # type: DataArray
        else:
            spec = get_spectrum('reference')
            tp = spec.peak_period(use_spline=True, monotone_interpolation=False)
            with open(name, "wb") as file_handle:
                pickle.dump(tp, file_handle)
            return tp
    else:
        raise Exception(f'unknown kind {kind}')


def get_displacements():
    """
    Get raw displacement data from netcdf file. Note that the displacement data is not include as part of the repository
    due to its size. Instead, derived spectra are included.

    :return: Pandas dataframe with displacement data.
    """
    dataset = xarray.open_dataset('./data/displacement.nc')
    return dataset.to_dataframe()