Skip to content

Commit 4ec2aff

Browse files
jmcvey3simmsa
jmcvey3
and
simmsa
authored
NDBC Updates (#437)
Opening this to solve Issue #427 I changed the source code to return degrees instead of radians, which was a simple fix. I also went ahead and updated the polar plots so that zero degrees is at the top of the figure, and positive degrees runs clockwise. We should update this in other polar plots in MHKiT too. I also updated units to pass CF conventions, whereafter I noticed that our tests for this code check that the units are the same as what NDBC outputs. I'm not sure that's necessary. I also saw that we're not actually testing the value output from the source code, which is important. --------- Co-authored-by: Simms, Andrew <andrew.simms@nlr.gov>
1 parent 352a99e commit 4ec2aff

8 files changed

Lines changed: 836 additions & 532 deletions

‎examples/directional_waves.ipynb‎

Lines changed: 410 additions & 199 deletions
Large diffs are not rendered by default.

‎examples/environmental_contours_example.ipynb‎

Lines changed: 61 additions & 66 deletions
Large diffs are not rendered by default.

‎examples/extreme_response_MLER_example.ipynb‎

Lines changed: 24 additions & 14 deletions
Large diffs are not rendered by default.

‎examples/extreme_response_contour_example.ipynb‎

Lines changed: 33 additions & 20 deletions
Large diffs are not rendered by default.

‎examples/extreme_response_full_sea_state_example.ipynb‎

Lines changed: 209 additions & 209 deletions
Large diffs are not rendered by default.

‎mhkit/tests/wave/io/test_ndbc.py‎

Lines changed: 68 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,6 @@
1111
import unittest
1212
import os
1313

14-
1514
testdir = dirname(abspath(__file__))
1615
datadir = normpath(join(testdir, "..", "..", "..", "..", "examples", "data", "wave"))
1716

@@ -100,8 +99,27 @@ def tearDownClass(self):
10099
# Realtime data
101100
def test_ndbc_read_realtime_met(self):
102101
data, units = wave.io.ndbc.read_file(join(datadir, "46097.txt"))
102+
cd = np.array(
103+
[
104+
180,
105+
12.0,
106+
np.nan,
107+
np.nan,
108+
np.nan,
109+
np.nan,
110+
np.nan,
111+
1001.0,
112+
6.1,
113+
10.9,
114+
np.nan,
115+
np.nan,
116+
np.nan,
117+
np.nan,
118+
],
119+
)
103120
expected_index0 = datetime(2019, 4, 2, 13, 50)
104121
self.assertSetEqual(set(data.columns), set(self.expected_columns_metRT))
122+
np.testing.assert_equal(data.iloc[-1].astype(float).values, cd)
105123
self.assertEqual(data.index[0], expected_index0)
106124
self.assertEqual(data.shape, (6490, 14))
107125
self.assertEqual(units, self.expected_units_metRT)
@@ -110,8 +128,26 @@ def test_ndbc_read_realtime_met(self):
110128
def test_ndbnc_read_historical_met(self):
111129
# QC'd monthly data, Aug 2019
112130
data, units = wave.io.ndbc.read_file(join(datadir, "46097h201908qc.txt"))
131+
cd = np.array(
132+
[
133+
160.0,
134+
2.7,
135+
np.nan,
136+
np.nan,
137+
np.nan,
138+
np.nan,
139+
np.nan,
140+
1014.9,
141+
14.5,
142+
13.4,
143+
np.nan,
144+
np.nan,
145+
np.nan,
146+
]
147+
)
113148
expected_index0 = datetime(2019, 8, 1, 0, 0)
114149
self.assertSetEqual(set(data.columns), set(self.expected_columns_metH))
150+
np.testing.assert_equal(data.iloc[-1].astype(float).values, cd)
115151
self.assertEqual(data.index[0], expected_index0)
116152
self.assertEqual(data.shape, (4464, 13))
117153
self.assertEqual(units, self.expected_units_metH)
@@ -127,14 +163,34 @@ def test_ndbc_read_spectral(self):
127163
def test_ndbc_read_cwind_no_units(self):
128164
data, units = wave.io.ndbc.read_file(join(datadir, "42a01c2003.txt"))
129165
self.assertEqual(data.shape, (4320, 5))
166+
cd = np.array(
167+
[
168+
[168.0, 6.1, np.nan, np.nan, np.nan],
169+
[153.0, 5.6, np.nan, np.nan, np.nan],
170+
[149.0, 4.8, np.nan, np.nan, np.nan],
171+
[147.0, 4.1, np.nan, np.nan, np.nan],
172+
[151.0, 4.1, 169.0, np.nan, np.nan],
173+
]
174+
)
175+
np.testing.assert_equal(data.tail().values, cd)
130176
self.assertEqual(units, None)
131177

132178
def test_ndbc_read_cwind_units(self):
133179
data, units = wave.io.ndbc.read_file(join(datadir, "46002c2016.txt"))
134180
self.assertEqual(data.shape, (28468, 5))
135-
self.assertEqual(units, wave.io.ndbc.parameter_units("cwind"))
181+
cd = np.array(
182+
[
183+
[303.0, 6.2, np.nan, np.nan, np.nan],
184+
[306.0, 7.8, np.nan, np.nan, np.nan],
185+
[302.0, 7.6, np.nan, np.nan, np.nan],
186+
[293.0, 6.4, np.nan, np.nan, np.nan],
187+
[287.0, 6.3, 307.0, 9.7, 1819.0],
188+
]
189+
)
190+
np.testing.assert_equal(data.tail().values, cd)
136191

137192
def test_ndbc_available_data(self):
193+
# Can't test data values because datafiles are periodically updating
138194
data = wave.io.ndbc.available_data("swden", buoy_number="46029")
139195
cols = data.columns.tolist()
140196
exp_cols = ["id", "year", "filename"]
@@ -158,14 +214,18 @@ def test__ndbc_parse_filenames(self):
158214
self.assertListEqual(fnames, self.filenames)
159215

160216
def test_ndbc_request_data(self):
217+
# Can't test data values because datafiles are periodically updating
161218
filenames = pd.Series(self.filenames[0])
162219
ndbc_data = wave.io.ndbc.request_data("swden", filenames, to_pandas=False)
163220
self.assertTrue(xr.Dataset(self.swden).equals(ndbc_data["1996"]))
164221

165222
def test_ndbc_request_data_from_dataframe(self):
223+
# Can't test data values because datafiles are periodically updating
166224
filenames = pd.DataFrame(pd.Series(data=self.filenames[0]))
167225
ndbc_data = wave.io.ndbc.request_data("swden", filenames)
168-
assert_frame_equal(self.swden, ndbc_data["1996"])
226+
assert_frame_equal(
227+
self.swden, ndbc_data["1996"], check_column_type=False, check_dtype=False
228+
)
169229

170230
def test_ndbc_request_data_filenames_length(self):
171231
with self.assertRaises(ValueError):
@@ -221,9 +281,10 @@ def test_ndbc_date_string_to_datetime(self):
221281
def test_ndbc_parameter_units(self):
222282
parameter = "swden"
223283
units = wave.io.ndbc.parameter_units(parameter)
224-
self.assertEqual(units[parameter], "(m*m)/Hz")
284+
self.assertEqual(units[parameter], "m^2/Hz")
225285

226286
def test_ndbc_request_directional_data(self):
287+
# Can't test data values because datafiles are periodically updating
227288
data = self.directional_data
228289
# correct 5 parameters
229290
self.assertEqual(len(data), 5)
@@ -236,12 +297,14 @@ def test_ndbc_request_directional_data(self):
236297
self.assertEqual(len(data.frequency), 47)
237298

238299
def test_ndbc_create_spread_function(self):
300+
# Can't test data values because datafiles are periodically updating
239301
directions = np.arange(0, 360, 2.0)
240302
spread = wave.io.ndbc.create_spread_function(self.directional_data, directions)
241303
self.assertEqual(spread.shape, (47, 180))
242304
self.assertEqual(spread.units, "1/Hz/deg")
243305

244306
def test_ndbc_create_directional_spectrum(self):
307+
# Can't test data values because datafiles are periodically updating
245308
directions = np.arange(0, 360, 2.0)
246309
spectrum = wave.io.ndbc.create_directional_spectrum(
247310
self.directional_data, directions
@@ -250,6 +313,7 @@ def test_ndbc_create_directional_spectrum(self):
250313
self.assertEqual(spectrum.units, "m^2/Hz/deg")
251314

252315
def test_plot_directional_spectrum(self):
316+
# Can't test data values because datafiles are periodically updating
253317
directions = np.arange(0, 360, 2.0)
254318
spectrum = wave.io.ndbc.create_spread_function(
255319
self.directional_data, directions

‎mhkit/wave/graphics.py‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -935,6 +935,9 @@ def plot_directional_spectrum(
935935

936936
a, f = np.meshgrid(np.deg2rad(spectrum.direction), spectrum.frequency)
937937
_, ax = plt.subplots(subplot_kw=dict(projection="polar"))
938+
# Set 0 degrees to be at the top, and positive angles to be clockwise
939+
ax.set_theta_zero_location("N")
940+
ax.set_theta_direction(-1)
938941
tmp = np.floor(np.min(spectrum.data) * 10) / 10
939942
color_level_min = tmp if (color_level_min is None) else color_level_min
940943
color_level_max = np.ceil(np.max(spectrum.data) * 10) / 10

‎mhkit/wave/io/ndbc.py‎

Lines changed: 28 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -57,7 +57,13 @@ def read_file(file_name, missing_values=["MM", 9999, 999, 99], to_pandas=True):
5757
metadata: dict or None
5858
Dictionary with {column name: units} key value pairs when the NDBC file
5959
contains unit information, otherwise None is returned
60+
61+
Notes
62+
-----
63+
NDBC documentation of measurement descriptions and units can be found here:
64+
https://www.ndbc.noaa.gov/faq/measdes.shtml
6065
"""
66+
6167
if not isinstance(file_name, str):
6268
raise TypeError(f"file_name must be of type str. Got: {type(file_name)}")
6369
if not isinstance(missing_values, list):
@@ -710,9 +716,9 @@ def parameter_units(parameter=""):
710716
}
711717
elif parameter == "cwind":
712718
units = {
713-
"WDIR": "degT",
719+
"WDIR": "degT", # degree True North
714720
"WSPD": "m/s",
715-
"GDR": "degT",
721+
"GDR": "degT", # degree True North
716722
"GST": "m/s",
717723
"GTIME": "hhmm",
718724
}
@@ -764,10 +770,10 @@ def parameter_units(parameter=""):
764770
"WWH": "m",
765771
"WWP": "sec",
766772
"SwD": "-",
767-
"WWD": "degT",
773+
"WWD": "degT", # degree True North
768774
"STEEPNESS": "-",
769775
"APD": "sec",
770-
"MWD": "degT",
776+
"MWD": "degT", # degree True North
771777
}
772778
elif parameter == "srad":
773779
units = {
@@ -777,13 +783,13 @@ def parameter_units(parameter=""):
777783
}
778784
elif parameter == "stdmet":
779785
units = {
780-
"WDIR": "degT",
786+
"WDIR": "degT", # degree True North
781787
"WSPD": "m/s",
782788
"GST": "m/s",
783789
"WVHT": "m",
784790
"DPD": "sec",
785791
"APD": "sec",
786-
"MWD": "degT",
792+
"MWD": "degT", # degree True North
787793
"PRES": "hPa",
788794
"ATMP": "degC",
789795
"WTMP": "degC",
@@ -797,11 +803,11 @@ def parameter_units(parameter=""):
797803
"PRES": "hPa",
798804
"PTIME": "hhmm",
799805
"WSPD": "m/s",
800-
"WDIR": "degT",
806+
"WDIR": "degT", # degree True North
801807
"WTIME": "hhmm",
802808
}
803809
elif parameter == "swden":
804-
units = {"swden": "(m*m)/Hz"}
810+
units = {"swden": "m^2/Hz"}
805811
elif parameter == "swdir":
806812
units = {"swdir": "deg"}
807813
elif parameter == "swdir2":
@@ -812,13 +818,13 @@ def parameter_units(parameter=""):
812818
units = {"swr2": ""}
813819
else:
814820
units = {
815-
"swden": "(m*m)/Hz",
821+
"swden": "m^2/Hz",
816822
"PRES": "hPa",
817823
"PTIME": "hhmm",
818-
"WDIR": "degT",
824+
"WDIR": "degT", # degree True North
819825
"WTIME": "hhmm",
820826
"DPD": "sec",
821-
"MWD": "degT",
827+
"MWD": "degT", # degree True North
822828
"ATMP": "degC",
823829
"WTMP": "degC",
824830
"DEWP": "degC",
@@ -834,7 +840,7 @@ def parameter_units(parameter=""):
834840
"WWH": "m",
835841
"WWP": "sec",
836842
"SwD": "-",
837-
"WWD": "degT",
843+
"WWD": "degT", # degree True North
838844
"STEEPNESS": "-",
839845
"APD": "sec",
840846
"RATE": "mm/h",
@@ -859,7 +865,7 @@ def parameter_units(parameter=""):
859865
"WSPD20": "m/s",
860866
"T": "-",
861867
"HEIGHT": "m",
862-
"GDR": "degT",
868+
"GDR": "degT", # degree True North
863869
"GST": "m/s",
864870
"GTIME": "hhmm",
865871
"DEP01": "m",
@@ -1026,7 +1032,7 @@ def request_directional_data(buoy, year):
10261032

10271033
data_dict["swden"].attrs = {
10281034
"units": "m^2/Hz",
1029-
"long_name": "omnidirecational spectrum",
1035+
"long_name": "omnidirectional spectrum",
10301036
"standard_name": "S",
10311037
"description": "Omnidirectional *sea surface elevation variance (m^2)* spectrum (/Hz).",
10321038
}
@@ -1062,7 +1068,7 @@ def request_directional_data(buoy, year):
10621068
return xr.Dataset(data_dict)
10631069

10641070

1065-
def _create_spectrum(data, frequencies, directions, name, units):
1071+
def _create_spectrum_dataarray(data, frequencies, directions, name, units):
10661072
"""
10671073
Create an xarray.DataArray for storing spectrum data with correct
10681074
dimensions, coordinates, names, and units.
@@ -1166,16 +1172,18 @@ def create_spread_function(data, directions):
11661172

11671173
r1 = data["swr1"].data.reshape(-1, 1)
11681174
r2 = data["swr2"].data.reshape(-1, 1)
1175+
# a1 = degrees CW from N, according to NDBC documentation
11691176
a1 = data["swdir"].data.reshape(-1, 1)
1177+
# a2 degrees CW from N, according to NDBC documentation
11701178
a2 = data["swdir2"].data.reshape(-1, 1)
1171-
a = directions.reshape(1, -1)
1179+
a = directions.reshape(1, -1) # degrees
11721180
spread = (
11731181
1
11741182
/ np.pi
11751183
* (0.5 + r1 * np.cos(np.deg2rad(a - a1)) + r2 * np.cos(2 * np.deg2rad(a - a2)))
11761184
)
1177-
spread = _create_spectrum(
1178-
spread, data.frequency.values, directions, name="Spread", units="1"
1185+
spread = _create_spectrum_dataarray(
1186+
np.rad2deg(spread), data.frequency.values, directions, name="Spread", units="1"
11791187
)
11801188
return spread
11811189

@@ -1209,9 +1217,9 @@ def create_directional_spectrum(data, directions):
12091217
spread = create_spread_function(data, directions).values
12101218
omnidirectional_spectrum = data["swden"].data.reshape(-1, 1)
12111219
spectrum = omnidirectional_spectrum * spread
1212-
spectrum = _create_spectrum(
1220+
spectrum = _create_spectrum_dataarray(
12131221
spectrum,
1214-
data.frequency.values,
1222+
data["frequency"].values,
12151223
directions,
12161224
name="Elevation variance",
12171225
units="m^2",

0 commit comments

Comments
 (0)