Create Hindcast Wind

Source Code

create_hindcast_wind.py
 1"""Read wave data from nora3 and generate metocean hindcast readable from SIMA"""
 2
 3import time
 4from pathlib import Path
 5from datetime import datetime
 6import numpy as np
 7import xarray as xr
 8from dmt.dmt_writer import DMTWriter
 9from metocean_api import ts
10import simapy.metocean.hindcast as hc
11
12
13def __create_hindcast(hc_name, hc_values, lat_pos, lon_pos):
14    # Time is given in Unix epoch
15    # convert dates to strings
16    dates = hc_values["time"].astype("datetime64[s]")
17    sdates = np.datetime_as_string(dates, unit="h", timezone="UTC").astype("|S")
18
19    hindcast = hc.Hindcast()
20    hindcast.name = hc_name.replace("-", "_")
21    hindcast.date = sdates
22    hindcast.latitude = lat_pos
23    hindcast.longitude = lon_pos
24
25    speed = hc_values["wind_speed"]
26    # Wind direction in fixed heights above surface,wind_from_direction	degrees clockwise from north (meteorological)
27    # SIMA: North East Down, wind coming from, same as MET
28    direction = hc_values["wind_direction"]
29    heights = hc_values["height"]
30    winds = []
31    for i, height in enumerate(heights):
32        wind = hc.StochasticWind()
33        level = height.values
34        wind.level = level
35        wind.name = f"{level}m"
36        wind.speed = speed[:, i].values
37        wind.direction = direction[:, i].values
38        winds.append(wind)
39
40    hindcast.wind = winds
41    return hindcast
42
43
44if __name__ == "__main__":
45    # https://www.met.no/publikasjoner/met-report  Section Storms in Sulafjord, wind waves and currents
46    # https://www.met.no/publikasjoner/met-report/_/attachment/inline/86f02fd2-a979-43ef-a17e-3bdba201e584:c70eb4b6ffe6f3b7b98016f4a0ebfc5ca501c766/MET-report-03-2024.pdf
47
48    # Storm March 14, 2017
49    sd = datetime(2017, 3, 14, 0, 0)
50    ed = datetime(2017, 3, 14, 23, 0)
51    start_date = sd.strftime("%Y-%m-%d")
52    end_date = ed.strftime("%Y-%m-%d")
53
54    locations = {
55        "Sulesund": {"lat": 62.402865086109195, "lon": 6.028359996993728},
56        "Kvitneset": {"lat": 62.421049661227585, "lon": 6.000482407215768},
57        "BuoyA": {"description": "Sulafjorden", "lat": 62.4263, "lon": 6.0447},
58        "BuoyD": {"description": "Breisundet", "lat": 62.4464, "lon": 5.9336},
59    }
60
61    # Start timing
62    start = time.time()
63    for name, location in locations.items():
64        lat_pos = location["lat"]
65        lon_pos = location["lon"]
66
67        product = "NORA3_wind_sub"
68
69        name = f"hindcast-{name}-{product}-{start_date}-{end_date}"
70
71        nc_file = f"./output/simamet/{name}.nc"
72
73        df_ts = ts.TimeSeries(
74            lon=lon_pos,
75            lat=lat_pos,
76            datafile=nc_file,
77            start_time=start_date,
78            end_time=end_date,
79            product=product,
80        )
81
82        df_ts.import_data(save_csv=False, save_nc=True, use_cache=True)
83
84        with xr.open_dataset(nc_file) as values:
85            values = values.sel(time=slice(sd, ed))
86            hindcast_data = __create_hindcast(
87                name, values, df_ts.lat_data, df_ts.lon_data
88            )
89            path = Path(f"./output/simamet/{name}.h5")
90            path.parent.mkdir(parents=True, exist_ok=True)
91            DMTWriter().write(hindcast_data, path)
92            print(f"Written to {path.resolve()}")
93
94    # End timing and print elapsed time
95    end = time.time()
96    print("Elapsed time: " + str(end - start) + " seconds")