wavedata.py

#! /usr/bin/env python
# -*- coding: utf-8 -*-
# vim:fenc=utf-8
#
# =============================================================
#  Copyright © 2017 Daniel Santiago <dpelaez@cicese.edu.mx>
#  Distributed under terms of the GNU/GPL license.
# =============================================================

"""
File:     _wavedata.py
Created:  2017-08-01 11:19
Modified: 2017-08-01 11:19

"""

# --- import libs ---
import numpy as np
import netCDF4 as nc
import datetime as dt
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
plt.ion()


def main(date, lon0, lat0):

    # --- get opendap data ---
    url = f"http://oos.soest.hawaii.edu/thredds/dodsC/hioos/model/wav/ww3/" + \
          f"runs/WaveWatch_III_Global_Wave_Model_RUN_{date}T00:00:00Z"

    data = nc.Dataset(url)
    lat = data.variables["lat"][:]
    lon = data.variables["lon"][:]

    lon0 = 360 + lon0 if lon0 < 0 else lon0
    ix_lat = (lat - lat0 >= 0).nonzero()[0][0]
    ix_lon = (lon - lon0 >= 0).nonzero()[0][0]

    time = nc.num2date(data.variables['time'][:], data.variables['time'].units)
    time -= dt.timedelta(hours=8)

    Hs  = np.squeeze(data.variables["Thgt"][:,:,ix_lat, ix_lon])
    Tp  = np.squeeze(data.variables["Tper"][:,:,ix_lat, ix_lon])
    Dr  = np.squeeze(data.variables["Tdir"][:,:,ix_lat, ix_lon])

    sHs = np.squeeze(data.variables["shgt"][:,:,ix_lat, ix_lon])
    sTp = np.squeeze(data.variables["sper"][:,:,ix_lat, ix_lon])
    sDr = np.squeeze(data.variables["sdir"][:,:,ix_lat, ix_lon])

    wHs = np.squeeze(data.variables["whgt"][:,:,ix_lat, ix_lon])
    wTp = np.squeeze(data.variables["wper"][:,:,ix_lat, ix_lon])
    wDr = np.squeeze(data.variables["wdir"][:,:,ix_lat, ix_lon])


    # --- plot time series ---
    fig = plt.figure(figsize=(6.2, 6.2))
    ax1 = fig.add_subplot(311)
    ax2 = fig.add_subplot(312, sharex=ax1)
    ax3 = fig.add_subplot(313, sharex=ax1)
    #
    ax1.plot(time,  Hs, lw=1.5, c="#0080FF", label="Total")
    ax1.plot(time, sHs, lw=1.5, c="#B18904", label="Swell")
    ax1.plot(time, wHs, lw=1.5, c="#61210B", label="Wind-sea")
    #
    ax2.plot(time,  Tp, '.', c="#0080FF")
    ax2.plot(time, sTp, '.', c="#B18904")
    ax2.plot(time, wTp, '.', c="#61210B")
    #
    ax3.plot(time,  Dr, '.', c="#0080FF")
    ax3.plot(time, sDr, '.', c="#B18904")
    ax3.plot(time, wDr, '.', c="#61210B")
    ax3.set_yticks([0, 90, 180, 270, 360])
    ax3.set_yticklabels(["W", "N", "E", "S", "W"])

    ax1.legend(loc=1, ncol=3, bbox_to_anchor=(1,1.25))
    ax1.set_ylabel("$H_s$ [m]")
    ax2.set_ylabel("$T_p$ [s]")
    ax3.set_ylabel("$\\theta$ [deg]")
    plt.setp(ax1.get_xticklabels(), visible=False)
    plt.setp(ax2.get_xticklabels(), visible=False)
    ax3.xaxis.set_major_locator(mdates.DayLocator())
    ax3.xaxis.set_minor_locator(mdates.HourLocator(interval=1))
    ax3.xaxis.set_major_formatter(mdates.DateFormatter("%d %h"))
    #
    fig.savefig("wave_{date}.pdf")


if __name__ == '__main__':

    # --- define params ---
    date = "2017-07-03"
    lat0 = 31.8
    lon0 = -117.
    main(date, lon0, lat0)