MHKiT-Software / MHKiT-Python

MHKiT-Python provides the marine renewable energy (MRE) community tools for data processing, visualization, quality control, resource assessment, and device performance.
https://mhkit-software.github.io/MHKiT/
BSD 3-Clause "New" or "Revised" License
51 stars 45 forks source link

xarray please! #78

Closed ryancoe closed 1 year ago

ryancoe commented 4 years ago

@ssolson knows that I'm a big xarray evangelist... I need to use some NDBC data and I think we should update the ndbc.py module to output xarray DataArray objects.

I know @ssolson has already implemented parts of the snippet below (e.g., replacing 999.00 with NaN) in ndbc.py, but I wasn't quite sure where to insert things so I just made a MWE. Should be fairly straightforward to incorporate from here.

import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

# save some data from NDBC as 'data.txt', in this case I used NDBC46050
# https://www.ndbc.noaa.gov/download_data.php?filename=46050w1996.txt.gz&dir=data/historical/swden/

# import data
df = pd.read_csv('data.txt', 
                 sep='\s+', 
                 parse_dates=[['YY', 'MM', 'DD', 'hh']], 
                 index_col=0)
da = df.to_xarray().to_array(dim='Frequency',
                             name='Spectral energy density')
da.attrs["units"] = "m$^2$/Hz"
da = da.rename({'YY_MM_DD_hh':'Date'})
da['Frequency'] = da['Frequency'].astype(np.float)
da['Frequency'].attrs['units'] = 'Hz'
da = da.where(da != 999.00)

# make a plot for a single day
fig, ax = plt.subplots(nrows=1, sharex=True, figsize=(7,4))
mdate = '1996-07-17'
dtp = da.sel(Date=slice(mdate,mdate))
dtp.plot.contourf()
dtp.plot.contour(colors='w', linewidths=0.5)
ax.set_title('Spectrogram for {}'.format(mdate))
fig.tight_layout()

# mean spectrum for that day
fig, ax = plt.subplots(nrows=1, sharex=True, figsize=(7,4))
pm = dtp.mean('Date').plot(ax=ax, color='b', linewidth=2, label='Mean')
ph = [x[1].plot(ax=ax, color='b', linewidth=0.25, alpha=0.5, label='Hourly') for x in dtp.groupby('Date')]
ax.set_ylim(bottom=0)
ax.set_title('Mean wave spectrum for {}'.format(mdate))
plt.autoscale(enable=True, axis='x', tight=True)
ax.legend(handles=[pm[0], ph[0][0]])

image image

<xarray.DataArray 'Spectral energy density' (Frequency: 38, Date: 5592)>
array([[      nan, 1.400e-01, 1.100e-01, ...,       nan, 1.000e-02,
        0.000e+00],
       [      nan, 8.800e-01, 4.600e-01, ...,       nan, 1.000e-02,
        0.000e+00],
       [      nan, 1.175e+01, 6.730e+00, ...,       nan, 2.000e-02,
        2.000e-02],
       ...,
       [      nan, 4.000e-02, 4.000e-02, ...,       nan, 2.000e-02,
        3.000e-02],
       [      nan, 3.000e-02, 3.000e-02, ...,       nan, 2.000e-02,
        2.000e-02],
       [      nan, 3.000e-02, 2.000e-02, ...,       nan, 2.000e-02,
        3.000e-02]])
Coordinates:
  * Date       (Date) datetime64[ns] 1996-01-01 ... 1996-08-21T23:00:00
  * Frequency  (Frequency) float64 0.03 0.04 0.05 0.06 ... 0.37 0.38 0.39 0.4
Attributes:
    units:    m$^2$/Hz
Alex-McVey commented 2 years ago

Issue #57 is a potential example of where this may be applicable

akeeste commented 1 year ago

See discussion in #260