Skip to content

Commit 59aa25b

Browse files
pieterbartsmitpieterbartsmit
authored
This pull requests ports the wind estimation from roguewave into roguewavespectrum. Windspeed and direction may now be estimated from methods on the spectrum object. (#9)
- The port does not preserve all original functionality. Specifically - it only implements the now default method by Shimura. - As simpe unit test to confirm results are equivalent is included in the unit tests.
1 parent e529b58 commit 59aa25b

10 files changed

Lines changed: 471 additions & 1 deletion

File tree

‎setup.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@
66

77
setuptools.setup(
88
name="roguewavespectrum",
9-
version="0.2.19",
9+
version="0.2.20",
1010
license="Apache 2 License",
1111
install_requires=[
1212
"numpy",

‎src/roguewavespectrum/_physical_constants.py‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@
1212
VONKARMAN_CONSTANT = 0.4
1313
DYNAMIC_SURFACE_TENSION_WATER = 0.073
1414
PHILLIPS_CONSTANT = 0.0081
15+
CHARNOCK_CONSTANT = 0.012
1516

1617

1718
class PhysicsOptions:
@@ -40,6 +41,7 @@ def __init__(
4041
gravity: float = GRAVITATIONAL_ACCELERATION,
4142
wave_type: str = "gravity",
4243
wave_regime: str = "intermediate",
44+
charnock_constant: float = CHARNOCK_CONSTANT,
4345
):
4446
"""
4547
:param density: Density of the fluid. Default is 1024 kg/m^3
@@ -59,6 +61,7 @@ def __init__(
5961
self.gravity = gravity
6062
self.wave_type = wave_type
6163
self.wave_regime = wave_regime
64+
self.charnock_constant = charnock_constant
6265

6366
@property
6467
def kinematic_surface_tension(self):
@@ -78,6 +81,7 @@ def kinematic_viscosity(self):
7881
gravity=GRAVITATIONAL_ACCELERATION,
7982
wave_type="gravity",
8083
wave_regime="intermediate",
84+
charnock_constant=CHARNOCK_CONSTANT,
8185
)
8286

8387

‎src/roguewavespectrum/_physics.py‎

Lines changed: 170 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,170 @@
1+
"""
2+
Contents: Wind Estimator
3+
4+
Copyright (C) 2022
5+
Sofar Ocean Technologies
6+
7+
Authors: Pieter Bart Smit
8+
"""
9+
import typing
10+
import numpy
11+
from typing import Literal
12+
from xarray import DataArray
13+
from ._physical_constants import PhysicsOptions
14+
15+
_direction_convention = Literal[
16+
"coming_from_clockwise_north", "going_to_counter_clockwise_east"
17+
]
18+
19+
20+
def friction_velocity(
21+
frequency: DataArray,
22+
variance_density: DataArray,
23+
a1: DataArray,
24+
b1: DataArray,
25+
physics_options: PhysicsOptions,
26+
power=4,
27+
directional_spreading_constant=2.5,
28+
beta=0.012,
29+
) -> DataArray:
30+
31+
e_peak, _, _ = equilibrium_range_values(
32+
frequency,
33+
variance_density,
34+
a1,
35+
b1,
36+
power=power,
37+
)
38+
emean = 8.0 * numpy.pi**3 * e_peak
39+
40+
# Get friction velocity from spectrum
41+
return emean / physics_options.gravity / directional_spreading_constant / beta / 4
42+
43+
44+
def estimate_wind_speed_from_wave_spectrum(
45+
frequency: DataArray,
46+
variance_density: DataArray,
47+
a1: DataArray,
48+
b1: DataArray,
49+
height_meter=10,
50+
power=4,
51+
directional_spreading_constant=2.5,
52+
phillips_constant_beta=0.0133,
53+
physics_options: PhysicsOptions = None,
54+
**kwargs,
55+
) -> DataArray:
56+
#
57+
# =========================================================================
58+
# Required Input
59+
# =========================================================================
60+
#
61+
# f :: frequencies (in Hz)
62+
# E :: Variance densities (in m^2 / Hz )
63+
#
64+
# =========================================================================
65+
# Output
66+
# =========================================================================
67+
#
68+
# U10 :: in m/s
69+
# Direction :: in degrees clockwise from North (where wind is *coming from)
70+
#
71+
# =========================================================================
72+
# Named Keywords (parameters to inversion algorithm)
73+
# =========================================================================
74+
# Npower = 4 :: exponent for the fitted f^-N spectral tail
75+
# I = 2.5 :: Philips Directional Constant
76+
# beta = 0.012 :: Equilibrium Constant
77+
# Kapppa = 0.4 :: Von Karman constant
78+
# Alpha = 0.012 :: Constant in estimating z0 from u* in Charnock relation
79+
# grav = 9.81 :: Gravitational acceleration
80+
#
81+
# =========================================================================
82+
# Algorithm
83+
# =========================================================================
84+
#
85+
# 1) Find the part of the spectrum that best fits a f^-4 shape
86+
# 2) Estimate the Phillips equilibrium level "Emean" over that range
87+
# 3) Use Emean to estimate Wind speed (using Charnock and LogLaw)
88+
# 4) Calculate mean direction over equilibrium range
89+
90+
# Get friction velocity from spectrum
91+
ustar = friction_velocity(
92+
frequency,
93+
variance_density,
94+
a1,
95+
b1,
96+
physics_options,
97+
power,
98+
directional_spreading_constant,
99+
phillips_constant_beta,
100+
)
101+
102+
# Find z0 from Charnock Relation
103+
z0 = charnock_roughness_length(ustar, physics_options, **kwargs)
104+
windspeed_at_height = (
105+
ustar / physics_options.vonkarman_constant * numpy.log(height_meter / z0)
106+
)
107+
return windspeed_at_height
108+
109+
110+
def estimate_wind_direction_from_wave_spectrum(
111+
frequency: DataArray,
112+
variance_density: DataArray,
113+
a1: DataArray,
114+
b1: DataArray,
115+
power=4,
116+
**kwargs,
117+
) -> DataArray:
118+
#
119+
120+
e_peak, a1_peak, b1_peak = equilibrium_range_values(
121+
frequency,
122+
variance_density,
123+
a1,
124+
b1,
125+
power=power,
126+
)
127+
128+
# Estimate direction from tail
129+
direction = (180.0 / numpy.pi * numpy.arctan2(b1_peak, a1_peak)) % 360
130+
return direction
131+
132+
133+
def equilibrium_range_values(
134+
frequency: DataArray,
135+
variance_density: DataArray,
136+
a1: DataArray,
137+
b1: DataArray,
138+
power=4,
139+
) -> typing.Tuple[DataArray, DataArray, DataArray]:
140+
"""
141+
:param spectrum:
142+
:param frequency:
143+
:param variance_density:
144+
:param a1:
145+
:param b1:
146+
:param power:
147+
:return:
148+
"""
149+
150+
scaled_spec = variance_density * frequency**power
151+
scaled_spec = scaled_spec.fillna(0)
152+
indices = scaled_spec.argmax(dim="frequency")
153+
a1_peak = a1.isel({"frequency": indices})
154+
b1_peak = b1.isel({"frequency": indices})
155+
e_peak = scaled_spec.isel({"frequency": indices})
156+
return e_peak, a1_peak, b1_peak
157+
158+
159+
def charnock_roughness_length(
160+
friction_velocity: DataArray, physics_options: PhysicsOptions, **kwargs
161+
) -> DataArray:
162+
163+
if not isinstance(friction_velocity, DataArray):
164+
friction_velocity = DataArray(data=friction_velocity)
165+
166+
return (
167+
physics_options.charnock_constant
168+
* friction_velocity**2
169+
/ physics_options.gravity
170+
)

‎src/roguewavespectrum/_spectrum.py‎

Lines changed: 123 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,11 @@
3030
get_angle_convention_and_unit,
3131
)
3232

33+
from ._physics import (
34+
estimate_wind_speed_from_wave_spectrum,
35+
estimate_wind_direction_from_wave_spectrum,
36+
)
37+
3338
from roguewavespectrum._estimators.estimate import (
3439
estimate_directional_spectrum_from_moments,
3540
Estimators,
@@ -2078,6 +2083,124 @@ def waveage(self, windspeed: typing.Union[Number, DataArray]) -> DataArray:
20782083
self.peak_wavespeed() / windspeed, "wave_age", overwrite=True
20792084
)
20802085

2086+
def estimate_wind_speed_at_10_meter(self, **kwargs) -> DataArray:
2087+
"""
2088+
Estimate the wind speed at 10 meter height from the wave spectrum using the method by Shimura et al. (2022)
2089+
as calibrated by Dorsay et al. (2023) for Spotter wave buoys. The method uses the high-frequency tail of the
2090+
wave spectrum to estimate the wind speed and assumes fully developed steady sea conditions. Nearby landmasses
2091+
might influence the accuracy of the estimate. It is reliably for wind speeds between approximately 5 m/s and 20
2092+
m/s. Outside this range caution is advised.
2093+
2094+
Shimura, T., Mori, N., Baba, Y., & Miyashita, T. (2022). Ocean surface wind estimation from waves based on small
2095+
GPS buoy observations in a bay and the open ocean. Journal of Geophysical Research: Oceans,
2096+
127(9), e2022JC018786.
2097+
2098+
Dorsay, C., Egan, G., Houghton, I., Hegermiller, C., & Smit, P. B. (2023). Proxy observations of surface wind
2099+
from a globally distributed network of wave buoys. Journal of Atmospheric and Oceanic Technology,
2100+
40(12), 1591-1603.
2101+
2102+
:param kwargs:
2103+
:return: Estimated wind speed at 10 meter height (m/s)
2104+
"""
2105+
2106+
# these are the defaults used for spotter buoys
2107+
maximum_tail_frequency = kwargs.get("maximum_tail_frequency", 0.5)
2108+
height_meter = kwargs.get("height_meter", 10.0)
2109+
power = kwargs.get("power", 4)
2110+
directional_spreading_constant = kwargs.get(
2111+
"directional_spreading_constant", 2.5
2112+
)
2113+
phillips_constant_beta = kwargs.get("phillips_constant_beta", 0.0133)
2114+
physics_options = kwargs.get("physics_options", None)
2115+
2116+
if physics_options is None:
2117+
# NOTE: Using default physics options with charnock constant overridden
2118+
physics_options = PhysicsOptions(
2119+
density=PHYSICSOPTIONS.density,
2120+
temperature_degrees=PHYSICSOPTIONS.temperature,
2121+
dynamic_viscosity=PHYSICSOPTIONS.dynamic_viscosity,
2122+
vonkarman_constant=PHYSICSOPTIONS.vonkarman_constant,
2123+
dynamic_surface_tension=PHYSICSOPTIONS.surface_tension,
2124+
gravity=PHYSICSOPTIONS.gravity,
2125+
wave_type=PHYSICSOPTIONS.wave_type,
2126+
wave_regime=PHYSICSOPTIONS.wave_regime,
2127+
charnock_constant=0.02,
2128+
)
2129+
2130+
spec1d = self.as_frequency_spectrum()
2131+
spec1d = spec1d.bandpass(fmax=maximum_tail_frequency + 1e-6)
2132+
2133+
wind_speed = estimate_wind_speed_from_wave_spectrum(
2134+
spec1d.frequency,
2135+
spec1d.variance_density,
2136+
spec1d.a1,
2137+
spec1d.b1,
2138+
height_meter,
2139+
power=power,
2140+
directional_spreading_constant=directional_spreading_constant,
2141+
phillips_constant_beta=phillips_constant_beta,
2142+
physics_options=physics_options,
2143+
)
2144+
wind_speed = wind_speed.drop_vars(names=NAME_F, errors="ignore")
2145+
wind_speed = set_conventions(wind_speed, "wind_direction", overwrite=True)
2146+
return wind_speed
2147+
2148+
def estimate_wind_direction(
2149+
self,
2150+
directional_unit: DirectionalUnit = "degree",
2151+
directional_convention: DirectionalConvention = "mathematical",
2152+
**kwargs,
2153+
) -> DataArray:
2154+
"""
2155+
Estimate the wind direction from the wave spectrum using the method by Shimura et al. (2022)
2156+
as calibrated by Dorsay et al. (2023) for Spotter wave buoys. The method uses the high-frequency tail of the
2157+
wave spectrum to estimate the wind speed and assumes fully developed steady sea conditions. Nearby landmasses
2158+
might influence the accuracy of the estimate. It is reliably for wind speeds between approximately 5 m/s and 20
2159+
m/s. Outside this range caution is advised.
2160+
2161+
Shimura, T., Mori, N., Baba, Y., & Miyashita, T. (2022). Ocean surface wind estimation from waves based on small
2162+
GPS buoy observations in a bay and the open ocean. Journal of Geophysical Research: Oceans,
2163+
127(9), e2022JC018786.
2164+
2165+
Dorsay, C., Egan, G., Houghton, I., Hegermiller, C., & Smit, P. B. (2023). Proxy observations of surface wind
2166+
from a globally distributed network of wave buoys. Journal of Atmospheric and Oceanic Technology,
2167+
40(12), 1591-1603.
2168+
2169+
:param directional_unit: units of the output angle, one of 'degree' or 'radians'
2170+
:param directional_convention: convention of the output angle, one of 'mathematical' or 'oceanographical' or
2171+
'meteorological' where:
2172+
- 'mathematical': 0 degrees is East, angles increase anticlockwise, direction is the direction the wind
2173+
points to.
2174+
- 'oceanographical': 0 degrees is North, angles increase clockwise, direction is the direction the wind
2175+
points to.
2176+
- 'meteorological': 0 degrees is North, angles increase clockwise, direction is the direction the wind
2177+
comes from.
2178+
2179+
:param kwargs:
2180+
:return: Estimated wind direction (degrees). By default in degrees anticlockwise from East.
2181+
"""
2182+
2183+
maximum_tail_frequency = kwargs.get("maximum_tail_frequency", 0.5)
2184+
spec1d = self.as_frequency_spectrum()
2185+
spec1d = spec1d.bandpass(fmax=maximum_tail_frequency + 1e-6)
2186+
power = kwargs.get("power", 4)
2187+
2188+
direction = estimate_wind_direction_from_wave_spectrum(
2189+
spec1d.frequency,
2190+
spec1d.variance_density,
2191+
spec1d.a1,
2192+
spec1d.b1,
2193+
power=power,
2194+
**kwargs,
2195+
)
2196+
direction = convert_angle_convention(
2197+
direction, directional_convention, "mathematical", directional_unit
2198+
)
2199+
direction = direction.drop_vars(names=NAME_F, errors="ignore")
2200+
2201+
direction = set_conventions(direction, "wind_direction", overwrite=True)
2202+
return direction
2203+
20812204
# ===================================================================================================================
20822205
# IO
20832206
# ===================================================================================================================

‎src/roguewavespectrum/_variable_names.py‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -394,6 +394,42 @@
394394
"comment": "Direction the waves are going to, measured anti-clockwise from East",
395395
},
396396
},
397+
"wind_direction": {
398+
"__default__": "mathematical",
399+
"meteorological": {
400+
"units": "degree",
401+
"long_name": "Surface wind direction",
402+
"missing_value": np.nan,
403+
"valid_min": 0,
404+
"valid_max": 360,
405+
"standard_name": "surface_wind_direction_at_ten_meter",
406+
"comment": "Direction the wind is coming from, measured clockwise from North",
407+
},
408+
"oceanographical": {
409+
"units": "degree",
410+
"long_name": "Surface wind direction",
411+
"missing_value": np.nan,
412+
"valid_min": 0,
413+
"valid_max": 360,
414+
"standard_name": "sea_surface_wave_to_direction_at_variance_spectral_density_maximum",
415+
"comment": "Direction the wind is going to, measured clockwise from North",
416+
},
417+
"mathematical": {
418+
"units": "degree",
419+
"long_name": "Surface wind direction",
420+
"missing_value": np.nan,
421+
"valid_min": 0,
422+
"valid_max": 360,
423+
"comment": "Direction the wind is going to, measured anti-clockwise from East",
424+
},
425+
},
426+
"wind_speed": {
427+
"units": "m/s",
428+
"long_name": "Surface wind speed",
429+
"missing_value": np.nan,
430+
"valid_min": 0,
431+
"standard_name": "surface_wind_speed_at_ten_meter",
432+
},
397433
"mean_directional_spread": {
398434
"units": "degree",
399435
"long_name": "Wave directional spread",

‎test/test_data/__init__.py‎

Whitespace-only changes.

‎test/test_data/wind_estimates/__init__.py‎

Whitespace-only changes.

‎test/test_data/wind_estimates/spotter_spectrum.nc‎

339 KB
Binary file not shown.

0 commit comments

Comments
 (0)