Source code for salishsea_tools.wind_tools

# Copyright 2013-2021 The Salish Sea MEOPAR contributors
# and The University of British Columbia

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#    http://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""Functions for working with wind model results and wind observations
data.
"""
from __future__ import division

from collections import namedtuple
from pathlib import Path
from datetime import datetime
import os

import numpy as np
import pandas as pd
from salishsea_tools import nc_tools, stormtools

# For convenience we import the wind speed conversion factors
# and functions and wind direction manipulation functions so that they
# can be used from either here or the
# :py:mod:`salishsea_tools.unit_conversions` module
from salishsea_tools.unit_conversions import (
    M_PER_S__KM_PER_HR,
    M_PER_S__KNOTS,
    mps_kph,
    mps_knots,
    wind_to_from,
    bearing_heading,
)


__all__ = [
    'calc_wind_avg_at_point',
    'M_PER_S__KM_PER_HR', 'M_PER_S__KNOTS', 'mps_kph', 'mps_knots',
    'wind_to_from', 'bearing_heading',
    'wind_speed_dir',
]


def get_EC_observations(station, start_day, end_day):
    """Gather hourly Environment Canada (EC) weather observations for the
    station and dates indicated.

    This function is a wrapper for stormtools.get_EC_observations, which
    should be migrated to this module at some point

    The hourly EC data is stored in monthly files, so only a single month can
    be downloaded at a time.

    :arg station: Station name (no spaces). e.g. 'PointAtkinson'
    :type station: str

    :arg start_day: Start date in the format '01-Dec-2006'.
    :type start_day: str

    :arg end_day: End date in the format '01-Dec-2006'.
    :type end_day: str

    :returns: wind_speed, wind_dir, temperature, times, lat and lon:
              wind speed is in m/s
              wind_dir is direction wind is blowing to in degrees measured
              counterclockwise from East
              temperature is in Kelvin
              time is UTC
              Also returns latitude and longitude of the station.
    """

    # Call get_EC_observations from stormtools
    wind_spd, wind_dir, temp, times, lat, lon = stormtools.get_EC_observations(
        station, start_day, end_day,
    )

    return wind_spd, wind_dir, temp, times, lat, lon


def parse_DFO_buoy_date(line):
    """Parse the DFO buoy date from a list of data block header strings

    :arg line: list of data block header strings (containing datetime info)

    :returns: parsed datetime object

    :rtype: :py:class:`datetime.datetime`
    """

    # The date format changes throughout the record
    # -- thus the multiple cases
    if (int(line[3]) > 2020) & (int(line[3]) < 202020):
        year, month, day = int(line[3][:4]), int(line[3][4:]), int(line[4])
        HHMM = f'{int(line[5]):04d}'
    elif int(line[3]) > 202020:
        year, month, day = int(line[3][:4]), int(line[3][4:6]), int(line[3][6:])
        HHMM = f'{int(line[4]):04d}'
    elif int(line[4]) > 12:
        year = int(line[3])
        MMDD, HHMM = [f'{int(l):04d}' for l in line[4:6]]
        month, day = int(MMDD[:2]), int(MMDD[2:])
    else:
        year = int(line[3])
        HHMM = f'{int(line[6]):04d}'
        month, day = int(line[4]), int(line[5])
    hour, minute = int(HHMM[:2]), int(HHMM[2:])
    date = datetime(year, month, day, hour, minute, 0)

    return date


def read_DFO_buoy(station, year):
    """Reads the data from a DFO *.fb buoy data file and appends to a timeseries dict

    :arg station: station name string

    :arg year: integer year requested

    :returns: wsp, wdir and time arrays

    :rtype: :py:class:`numpy.ndarray`
    """

    # Station ID dict
    station_ids = {
        'Halibut Bank': 46146,
        'Sentry Shoal': 46131,
    }

    # Data url
    url = 'http://www.meds-sdmm.dfo-mpo.gc.ca/alphapro/wave/waveshare/fbyears'

    # Open the *.zip file from url using Pandas
    ID = f'C{station_ids[station]}'
    file = os.path.join(url, ID, f'{ID.lower()}_{year}.zip')
    csv = pd.read_csv(file, header=None)

    # Initialize parsing booleans
    gotdate, gotwind = True, True

    # Initialize storage dict
    wspd, wdir, time = [], [], []

    # Read file line by line
    for line in csv.values:

        # Parse line to list of strings
        line_parsed = line[0].strip().split()

        # Ignore lines shorter than 4
        if len(line_parsed) > 3:

            # Begin new entry
            if line_parsed[3] == ID:
                gotdate = False

            # Parse date
            elif gotdate == False:
                gotdate, gotwind = True, False
                time.append(parse_DFO_buoy_date(line_parsed))

            # Read wind data
            elif gotwind == False:
                gotwind = True
                wdir.append(float(line_parsed[0].split('W')[0]))
                wspd.append(float(line_parsed[1].split('W')[0]))

    # Transform angle to deg CCW from east
    wdir = 270 - np.array(wdir)
    wdir[wdir < 0] += 360

    return np.array(wspd), wdir, np.array(time)


[docs] def wind_speed_dir(u_wind, v_wind): """Calculate wind speed and direction from u and v wind components. :kbd:`u_wind` and :kbd:`v_wind` may be either scalar numbers or :py:class:`numpy.ndarray` objects, and the elements of the return value will be of the same type. :arg u_wind: u-direction component of wind vector. :arg v_wind: v-direction component of wind vector. :returns: 2-tuple containing the wind speed and direction. The :py:attr:`speed` attribute holds the wind speed(s), and the :py:attr:`dir` attribute holds the wind direction(s). :rtype: :py:class:`collections.namedtuple` """ speed = np.sqrt(u_wind**2 + v_wind**2) dir = np.arctan2(v_wind, u_wind) dir = np.rad2deg(dir + (dir < 0) * 2 * np.pi) speed_dir = namedtuple('speed_dir', 'speed, dir') return speed_dir(speed, dir)
[docs] def calc_wind_avg_at_point(date_time, weather_path, windji, avg_hrs=-4): """Calculate the average wind components for a time period at a single weather dataset grid point. .. note:: The present implementation only supports averaging over a time period ending at:kbd:`date_time` and values for :kbd:`avg_hrs` in the range :kbd:`[-24, 0]`. :arg date_time: Date and time to use as base for the calculation. :type date_time: :py:class:`arrow.arrow.Arrow` :arg str weather_path: The directory where weather dataset files are stored. :arg 2-tuple windji: Indices of weather dataset grid point to calculate the average at as :kbd:`(lon_index, lat_index)`. :arg float avg_hrs: Number of hours to calculate :returns: 2-tuple containing the averaged wind components. The :py:attr:`u` attribute holds the u-direction component, and the :py:attr:`v` attribute holds the v-direction component. :rtype: :py:class:`collections.namedtuple` :raises: :py:exc:`IndexError` if :kbd:`avg_hrs` is outside the range :kbd:`[-24, 0]`. """ weather_filename_tmpl = 'hrdps_y{0.year:4d}m{0.month:02d}d{0.day:02d}.nc' try: weather_file = Path( weather_path, weather_filename_tmpl.format(date_time)) grid_weather = nc_tools.dataset_from_path(weather_file) except IOError: weather_file = Path( weather_path, 'fcst', weather_filename_tmpl.format(date_time)) grid_weather = nc_tools.dataset_from_path(weather_file) wind_u, wind_v, wind_t = nc_tools.uv_wind_timeseries_at_point( grid_weather, *windji) if date_time.hour < abs(avg_hrs): grid_weather = nc_tools.dataset_from_path( weather_file.with_name( weather_filename_tmpl.format(date_time.shift(days=-1)))) wind_prev_day = nc_tools.uv_wind_timeseries_at_point( grid_weather, *windji) wind_u = np.concatenate((wind_prev_day.u, wind_u)) wind_v = np.concatenate((wind_prev_day.v, wind_v)) wind_t = np.concatenate((wind_prev_day.time, wind_t)) i_date_time = np.where(wind_t == date_time.floor('hour'))[0].item() i_date_time_p1 = i_date_time + 1 u_avg = np.mean(wind_u[(i_date_time_p1 + avg_hrs):i_date_time_p1]) v_avg = np.mean(wind_v[(i_date_time_p1 + avg_hrs):i_date_time_p1]) wind_avg = namedtuple('wind_avg', 'u, v') return wind_avg(u_avg, v_avg)