# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # IMOVE Time-Series and Big-Data Workshop
# In this workshop we will use the time-series functionality of [Pandas](https://pandas.pydata.org/docs/) and [Xarray](http://xarray.pydata.org/en/stable/) to explore some data collected by the [Ocean Observatories Intiative](https://oceanobservatories.org/). Hopefully we will also get a chance to explore [Dask](https://docs.dask.org/en/latest/) and [Dask Delayed](https://docs.dask.org/en/latest/delayed.html) functions to parallelize a data analysis workflow in the cloud. We will be working on the [OOICloud](https://www.ooicloud.org/) [Pangeo](http://pangeo.io/) instance. Further information on using Python to analyze Earth science datasets can be found in the book [Earth and Environmental Data Science](https://earth-env-data-science.github.io/intro) which I have been using to teach Research Computing in the Earth Sciences this semester.

# ## Bottom pressure data at Axial Seamount
# Here we find data using the new OOI Data Explorer and use Pandas [read_csv()](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.read_csv.html) and [time-series functionality](https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html) to plot a smoothed representation of the bottom pressure at Axial Seamount. 

# ### First import some required packages

import pandas as pd
from matplotlib import pyplot as plt
# %matplotlib inline
plt.rcParams['figure.figsize'] = (14, 8)

from IPython import display
display.set_matplotlib_formats('retina')

# ### Next find some data on the OOI Data Explorer

url = 'http://erddap.dataexplorer.oceanobservatories.org/erddap/tabledap/ooi-rs03ccal-mj03f-05-botpta301.csv?time%2Cbotpres%2Cbotpres_qc_agg%2Cz&time%3E%3D2014-08-29T20%3A59%3A00Z&time%3C%3D2020-11-21T06%3A00%3A00Z'

url

botpt = pd.read_csv(url, parse_dates=True, usecols = ['time','botpres'], index_col='time', skiprows=[1])
botpt.head()

botpt = botpt.rename(columns={'botpres':'pressure'})

botpt.head()

len(botpt)

botpt[0:20000].pressure.plot(ylabel='pressure')

botpt['2014'].plot()

botpt['2019-04-28':'2019-07-12'].plot()

start = pd.Timestamp('2020-08')
botpt[start:start+pd.Timedelta(days=30)].plot()

botpt.plot()

botpt_rolling = botpt.rolling('14D', min_periods = 60*24*7).pressure.mean()
botpt_rolling.plot(ylabel='pressure')

botpt_rolling = botpt.rolling(60*24*14, win_type='hamming', min_periods = 60*24*7).pressure.mean()
botpt_rolling.plot(ylabel='pressure')

import hvplot.pandas

botpt_rolling[::60].hvplot(ylabel='pressure')

# ## Earthquake catalog from the OOI seismic array at Axial Seamount
# Here we parse and plot Axial Seamount earthquake catalog data from [<NAME>'s near-real-time automated earthquake location system](http://axial.ocean.washington.edu/). The data we will use is a text file in they HYPO71 output format located here: http://axial.ocean.washington.edu/hypo71.dat.

eqs_url = 'http://axial.ocean.washington.edu/hypo71.dat'

col_names = ['ymd', 'hm', 's', 'lat_deg', 'lat_min', 'lon_deg', 'lon_min',
        'depth', 'MW', 'NWR', 'GAP', 'DMIN',  'RMS',  'ERH', 'ERZ', 'ID', 'PMom', 'SMom']

eqs = pd.read_csv(eqs_url, sep='\s+', header=0, names=col_names)

eqs.head()

from datetime import datetime
def parse_hypo_date(ymd, hm, s):
    hour = int(hm.zfill(4)[0:2])
    minute = int(hm.zfill(4)[2:])
    second = float(s)
    if second == 60:
        second = 0
        minute += 1
    if minute == 60:
        minute=0
        hour +=1
    eq_date_str = ('%s%02.0f%02.0f%05.2f' % (ymd, hour, minute, second))
    return datetime.strptime(eq_date_str, '%Y%m%d%H%M%S.%f')


eqs = pd.read_csv(eqs_url, sep='\s+', header=0, names=col_names, parse_dates=[[0,1,2]], date_parser=parse_hypo_date)

eqs.head()

eqs['lat'] = eqs.lat_deg+eqs.lat_min/60
eqs['lon'] = -(eqs.lon_deg+eqs.lon_min/60)

eqs.head()

eqs.rename(columns={'ymd_hm_s': 'time', 'MW': 'mw'}, inplace=True)

eqs.set_index('time', inplace=True)

eqs.head()

eqs = eqs[['lat', 'lon', 'depth', 'mw']]

eqs.head()

len(eqs)

eqs.mw.plot(marker='.', linestyle='', markersize=1)

daily_count = eqs.mw.resample('1D').agg('count')
daily_count.head()

fig, ax1 = plt.subplots()
ax1.bar(daily_count.index, daily_count.values, width=5)
ax1.set_ylim(ymax=3000)

fig, ax1 = plt.subplots()
ax1.bar(daily_count.index, daily_count.values, width=5)
ax1.set_ylim(ymax=2500)
ax2 = ax1.twinx()
ax2.plot(botpt_rolling, color='r')

# +
fig, ax1 = plt.subplots()
ax1.bar(daily_count['2015'].index, daily_count['2015'].values, width=1)
ax1.set_ylim(ymax=2500)
ax1.set_ylabel('count')

ax2 = ax1.twinx()
ax2.plot(botpt_rolling['2015'], color='r')
ax2.set_ylabel('pressure')
# -

# ### Mapping eq data
# Let's make some maps just because we can.

import cartopy.crs as ccrs
import cartopy
import cartopy.feature as cfeature
import numpy as np

caldera = pd.read_csv('caldera.csv', sep=',')

now = pd.Timestamp('now')

eqs_sub = eqs[(now-pd.Timedelta(weeks=8)):]

eqs_sub

eqs_sub.index.dayofyear

# +
ax = plt.axes(projection=ccrs.Robinson(central_longitude=-130))
ax.plot(caldera.lon, caldera.lat,transform=ccrs.Geodetic(), c=(0.8, 0.8, 0.8))
sc = ax.scatter(eqs_sub.lon, eqs_sub.lat,
                s=0.00001*(eqs_sub.mw+3)**11, c=eqs_sub.index.dayofyear,
                edgecolor='k', cmap='viridis',
                transform=ccrs.PlateCarree())
plt.colorbar(sc, label='Day of Year')
ax.gridlines()
ax.set_title('Eqs from last 8 weeks');

extent = [-130.07, -129.95, 45.90, 46.02]
ax.set_extent(extent)
# -

eqs_sub = eqs['2016']

# +
ax = plt.axes(projection=ccrs.Robinson(central_longitude=-130))
ax.plot(caldera.lon, caldera.lat,transform=ccrs.PlateCarree())
sc = ax.scatter(eqs_sub.lon, eqs_sub.lat,
                c=eqs_sub.mw, s=40,
                edgecolor='k', cmap='Reds',
                transform=ccrs.PlateCarree())
plt.colorbar(sc, label='magnitude')
ax.gridlines()
ax.set_title('Eqs from 2016');

extent = [-130.07, -129.95, 45.90, 46.02]
ax.set_extent(extent)
# -

# ## OOI Seafloor Camera Data
# Now let's look at some video data from the [OOI Seafloor Camera](https://oceanobservatories.org/instrument-class/camhd/) system deployed at Axial Volcano on the Juan de Fuca Ridge. We will make use of the [Pycamhd](https://github.com/tjcrone/pycamhd) library, which can be used to extract frames from the ProRes encoded Quicktime files. These data are hosted on Microsoft's [Azure Open Datasets](https://azure.microsoft.com/en-us/services/open-datasets/).

import pycamhd as camhd
import numpy as np
from ipywidgets import interact
from ipywidgets import IntSlider

dbcamhd_url = 'https://ooiopendata.blob.core.windows.net/camhd/dbcamhd.json'

dbcamhd = pd.read_json(dbcamhd_url, orient="records", lines=True)

dbcamhd.tail()

dbcamhd.frame_count.sum()

mov = dbcamhd.iloc[7000]
mov


def show_image(frame_number):
    plt.rc('figure', figsize=(12, 6))
    plt.rcParams.update({'font.size': 8})
    frame = camhd.get_frame(mov.url, frame_number)
    fig, ax = plt.subplots();
    im1 = ax.imshow(frame);
    plt.yticks(np.arange(0,1081,270))
    plt.xticks(np.arange(0,1921,480))
    plt.title('Deployment: %s    File: %s    Frame: %s' % (mov.deployment, mov['name'], frame_number));


# %matplotlib inline
# %config InlineBackend.figure_format = 'svg'
initial_frame = 8060
frame_slider = IntSlider(min=0, max=mov.frame_count-1, step=1, value=initial_frame, continuous_update=False)
interact(show_image, frame_number=frame_slider);





# ## Scratch scratch

df = pd.DataFrame(index=test.groupby(test.index.date).count().index)
df['count'] = test.groupby(test.index.date).count().lat

counts = pd.DataFrame(index=test.resample('1D').agg('count').index)
#counts['count'] = test.resample('1D').agg('count')

test.nlargest(100, 'MW').plot.scatter(x='lon', y='lat', s='MW')

test['2020-01-01':].plot.scatter(x = 'lon', y = 'lat', c = 'depth', s='MW', cmap='magma')

test.hvplot.scatter(x = 'lon', y = 'lat', c='depth', datashade=True, dynspread=True)

eqs = eqs.rename(columns={''}
botpt = botpt.rename(columns={'botpres':'pressure'})

eqs = pd.read_csv(eq_url, sep='\s+', header=0, names=col_names, parse_dates=[[0,1,2]], date_parser=parse_hypo_date)

url = 'http://erddap.dataexplorer.oceanobservatories.org/erddap/tabledap/ooi-rs03int1-mj03c-10-trhpha301.csv?time%2Ctrhphte_abs%2Ctrhphte_abs_qc_agg%2Cz&time%3E%3D2014-09-27T06%3A46%3A00Z&time%3C%3D2020-11-21T05%3A47%3A00Z'

test2 = pd.read_csv(url, parse_dates=True, index_col='time', skiprows=[1])

test2.nsmallest(1000, 'trhphte_abs')

test2.iloc[-200000:].hvplot(y='trhphte_abs', marker='.')

url = 'http://erddap.dataexplorer.oceanobservatories.org/erddap/tabledap/ooi-rs03int2-mj03d-06-botpta303.nc?time%2Cbotpres%2Cbotpres_qc_agg%2Cz&time%3E%3D2014-08-29T23%3A17%3A00Z&time%3C%3D2020-11-21T06%3A00%3A00Z'

StarcoderdataJupyterScriptsDedupFiltered