<jupyter_start><jupyter_text>Development of Algorithms for ERDDAP Dataset from 2nd Narrows HADCP ObservationsFiguring out how to transform 2nd Narrows horizontal ADCP observations fromAIS CSV files into a netCDF file that is part of an ERDDAP dataset:* Read CSV file into `pandas` dataframe* Filter to include only 2nd Narrow HADCP rows* Convert to `xarray.Dataset`* Add metadata* **BONUS:** Add observations from an hourly file to an existing netCDF fileso that we can have daily or monthly netCDF instead of hourly<jupyter_code>from pathlib import Path

import arrow
import pandas
import numpy
import xarray
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Inspect a Sample `csv` File<jupyter_code>!wc -l /opp/observations/AISDATA/20180928T120000Z-20180928T130000Z.csv
!head /opp/observations/AISDATA/20180928T120000Z-20180928T130000Z.csv
!tail /opp/observations/AISDATA/20180928T120000Z-20180928T130000Z.csv
!wc -l /opp/observations/AISDATA/20180928T120000Z-20180928T130000Z.csv<jupyter_output>229 /opp/observations/AISDATA/20180928T120000Z-20180928T130000Z.csv<jupyter_text>Read `csv` File into a `panada.Dataframe`<jupyter_code>all_df = pandas.read_csv('/opp/observations/AISDATA/20180928T120000Z-20180928T130000Z.csv', skiprows=3)
all_df<jupyter_output><empty_output><jupyter_text>Isolate 2nd Narrows Observations* Drop duplicated rows* Use the `Time` column as the index<jupyter_code>hadcp_df = all_df.loc[all_df.MMSI == 3160171].drop_duplicates().set_index('Time')
hadcp_df.index = pandas.to_datetime(hadcp_df.index, format='%d/%m/%Y %H:%M')
hadcp_df<jupyter_output><empty_output><jupyter_text>Isolate the Seawater Current Observations* Drop the columns that are full of `N/A`s* Drop other columns we're not interested in by name* Rename the index to `time`* Simplify the current speed and direction column names<jupyter_code>xr_df = hadcp_df \
    .dropna(axis='columns', how='all') \
    .drop(['Name', 'MMSI', 'Air Temp.', 'Water Temp.', 'Water Level'], axis='columns')
xr_df.index.name = 'time'
xr_df.columns = ('speed', 'direction')
xr_df<jupyter_output><empty_output><jupyter_text>Store Observations in netCDF4 File* transform the `pandas.Dataframe` into an `xarray.Dataset`* define encoding for netCDF4 serialization* write to file<jupyter_code>ds = xarray.Dataset.from_dataframe(xr_df)
ds<jupyter_output><empty_output><jupyter_text>Compress the speed and direction observation values by scaling themto that they are stored with a precision of 0.1 in 16-bit integersand with `NaN`s converted to `-9999`.Use the default `proleptic_gregorian` calendar for times andstore them as seconds since the Unix epoch.Set the `time` dimension as unlimited to make concatenation easy.<jupyter_code>encoding={
    'speed': {'dtype': 'int16', 'scale_factor': 0.1, '_FillValue': -9999},
    'direction': {'dtype': 'int16', 'scale_factor': 0.1, '_FillValue': -9999},
    'time': {'dtype': 'float', 'units': 'seconds since 1970-01-01T00:00:00Z'},
}
ds.to_netcdf('/tmp/foo.nc', mode='w', encoding=encoding, unlimited_dims='time')<jupyter_output><empty_output><jupyter_text>Inspect the netCDF4 File<jupyter_code>!ls -lh /tmp/foo.nc
!ncdump -cst /tmp/foo.nc
currents = xarray.open_dataset('/tmp/foo.nc')
currents
currents.speed.plot(marker='o', linestyle='None')
currents.direction.plot(marker='o', linestyle='None')<jupyter_output><empty_output><jupyter_text>Concatenate Observations from Another `csv` to the netCDF4 File<jupyter_code>def csv_to_dataset(csv_file):
    df = pandas.read_csv(csv_file, skiprows=3)
    df = df.loc[df.Name == 'VFPA CN RAIL BRIDGE'].drop_duplicates().set_index('Time')
    df = df.dropna(axis='columns', how='all')
    df = df.drop(['Name', 'MMSI', 'Air Temp.', 'Water Temp.', 'Water Level'], axis='columns')
    df.index = pandas.to_datetime(df.index, format='%d/%m/%Y %H:%M')
    df.index.name = 'time'
    df.columns = ('speed', 'direction')
    return xarray.Dataset.from_dataframe(df)
ds = csv_to_dataset('/opp/observations/AISDATA/20180928T130000Z-20180928T140000Z.csv')
ds
currents = xarray.open_dataset('/tmp/foo.nc')
extended_currents = xarray.concat((currents, ds), dim='time')
extended_currents
extended_currents.to_netcdf('/tmp/bar.nc', mode='w', encoding=encoding, unlimited_dims='time')<jupyter_output><empty_output><jupyter_text>File only increases in size by 120 bytes!!That bodes well for the idea of storing month-long netCDF4 files for ERDDAP.<jupyter_code>!ls -l /tmp/bar.nc
!ls -l /tmp/foo.nc
!ncdump -cst /tmp/bar.nc
currents = xarray.open_dataset('/tmp/bar.nc')
currents
currents.speed.plot(marker='o', linestyle='None')
currents.direction.plot(marker='o', linestyle='None')<jupyter_output><empty_output><jupyter_text>Assemble a Month-long netCDF4 File<jupyter_code>def csv_to_dataset(csv_file):
    df = pandas.read_csv(csv_file, skiprows=3)
    df = df.loc[df.MMSI == 3160171].drop_duplicates().set_index('Time')
    df = df.dropna(axis='columns', how='all')
    for col in ['Name', 'MMSI', 'Air Temp.', 'Water Temp.', 'Water Level']:
        try:
            df = df.drop([col], axis='columns')
        except KeyError:
            pass
    df.index = pandas.to_datetime(df.index, format='%d/%m/%Y %H:%M')
    df.index.name = 'time'
    df.columns = ('speed', 'direction')
    return xarray.Dataset.from_dataframe(df)
!rm /tmp/baz.nc
attrs = {
    'speed': {
        'name': 'speed',
        'units': 'm/s',
        'ioos_category': 'currents',
        'standard_name': 'sea_water_speed',
        'long_name': 'Current Speed',
    },
    'direction': {
        'name': 'direction',
        'units': 'degree',
        'ioos_category': 'currents',
        'standard_name': 'sea_water_direction',
        'long_name': 'Current Direction',
    },
}
ais_path = Path('/opp/observations/AISDATA')
datetime_fmt = 'YYYYMMDDTHHmmss'
utc_start_hr = arrow.get('2018-09-01 00:00')
utc_offset = utc_start_hr.to('local').utcoffset()
csv_filename = f'{utc_start_hr.format(datetime_fmt)}Z-{utc_start_hr.shift(hours=+1).format(datetime_fmt)}Z.csv'
print(csv_filename)
ds = csv_to_dataset(ais_path/csv_filename)
ds['time'] -= pandas.to_timedelta(utc_offset)
for var in attrs:
    ds[var].attrs = attrs[var]
ds.to_netcdf('/tmp/baz.nc', mode='w', encoding=encoding, unlimited_dims='time')
ds.close()
for hr in arrow.Arrow.range('hour', start=utc_start_hr.shift(hours=+1), end=utc_start_hr.shift(months=+1, hours=-1)):
    with xarray.open_dataset('/tmp/baz.nc') as stored_obs:
        csv_filename = f'{hr.format(datetime_fmt)}Z-{hr.shift(hours=+1).format(datetime_fmt)}Z.csv'
        new_obs = csv_to_dataset(ais_path/csv_filename)
        new_obs['time'] -= pandas.to_timedelta(utc_offset)
        for var in attrs:
            new_obs[var].attrs = attrs[var]
        extended_obs = xarray.concat((stored_obs, new_obs), dim='time')
    # Drop repeated times because some files contain hh:00 to hh+1:00 instead of ending at hh:58
    _, index = numpy.unique(extended_obs.time.values, return_index=True)
    extended_obs = extended_obs.isel(time=index)
    extended_obs.to_netcdf('/tmp/baz.nc', mode='w', encoding=encoding, unlimited_dims='time')<jupyter_output><empty_output><jupyter_text>Inspect the netCDF File<jupyter_code>!ls -l /tmp/baz.nc
!ls -lh /tmp/baz.nc
with xarray.open_dataset('/tmp/baz.nc') as currents:
    print(currents)
    print(currents.time.min())
    print(currents.time.max())
    print(currents.speed)
    print(currents.direction)
    currents.speed.plot(marker='o', linestyle='None')<jupyter_output><xarray.Dataset>
Dimensions:    (time: 21073)
Coordinates:
  * time       (time) datetime64[ns] 2018-09-01 ... 2018-09-30T23:58:00
Data variables:
    speed      (time) float32 ...
    direction  (time) float32 ...
<xarray.DataArray 'time' ()>
array('2018-09-01T00:00:00.000000000', dtype='datetime64[ns]')
<xarray.DataArray 'time' ()>
array('2018-09-30T23:58:00.000000000', dtype='datetime64[ns]')
<xarray.DataArray 'speed' (time: 21073)>
array([1.2, 1.2, 1.1, ..., 0.5, 0.6, 0.6], dtype=float32)
Coordinates:
  * time     (time) datetime64[ns] 2018-09-01 ... 2018-09-30T23:58:00
Attributes:
    name:           speed
    units:          m/s
    ioos_category:  currents
    standard_name:  sea_water_speed
    long_name:      Current Speed
<xarray.DataArray 'direction' (time: 21073)>
array([85., 83., 83., ..., 91., 92., 93.], dtype=float32)
Coordinates:
  * time     (time) datetime64[ns] 2018-09-01 ... 2018-09-30T23:58:00
Attributes:
    name:           direction
    units:          degree
   [...]<jupyter_text>Initialize ERDDAP Dataset with 2017-11-20 14:00 Hour Observations2017-11-20 14:00 Pacific time is when Michael started the hourly storage of `csv` filesfrom the Vancouver Fraser Port Authority (VFPA) AIS feed.<jupyter_code>ais_path = Path('/opp/observations/AISDATA')
datetime_fmt = 'YYYYMMDDTHHmmss'

utc_start_hr = arrow.get('2017-11-20 22:00')
utc_offset = utc_start_hr.to('local').utcoffset()
netcdf_filename = f'VFPA_2ND_NARROWS_HADCP_2s_{utc_start_hr.format("YYYYMM")}.nc'
var_attrs = {
    'speed': {
        'name': 'speed',
        'units': 'm/s',
        'ioos_category': 'currents',
        'standard_name': 'sea_water_speed',
        'long_name': 'Current Speed',
    },
    'direction': {
        'name': 'direction',
        'units': 'degree',
        'ioos_category': 'currents',
        'standard_name': 'sea_water_direction',
        'long_name': 'Current Direction',
    },
}
dataset_attrs = {
    'cdm_data_type': 'TimeSeries',
    'cdm_timeseries_variables': 'speed, direction',
    'infoUrl': 'https://salishsea-meopar-tools.readthedocs.org/en/latest/results_server/',
    'acknowledgement': 'MEOPAR, DFO IOS, VFPA',
    'institution': 'UBC EOAS, DFO IOS, VFPA',
}
!rm '/tmp/{netcdf_filename}'
csv_filename = f'{utc_start_hr.format(datetime_fmt)}Z-{utc_start_hr.shift(hours=+1).format(datetime_fmt)}Z.csv'
print(csv_filename)
ds = csv_to_dataset(ais_path/csv_filename)
ds['time'] -= pandas.to_timedelta(utc_offset)
ds.coords['latitude'] = 49.2938889
ds.coords['longitude'] = -123.0247222
for var in var_attrs:
    ds[var].attrs = attrs[var]
ds.attrs = dataset_attrs
ds.to_netcdf(f'/tmp/{netcdf_filename}', mode='w', encoding=encoding, unlimited_dims='time')
ds.close()
!ls -lh /tmp/VFPA_2ND_NARROWS_HADCP_2s_201711.nc
with xarray.open_dataset(f'/tmp/{netcdf_filename}') as currents:
    print(currents)
    print(currents.time.min())
    print(currents.time.max())
    print(currents.speed)
    print(currents.direction)
    currents.speed.plot(marker='o', linestyle='None')<jupyter_output><xarray.Dataset>
Dimensions:    (time: 31)
Coordinates:
  * time       (time) datetime64[ns] 2017-11-20T22:00:00 ... 2017-11-20T23:00:00
    latitude   float64 ...
    longitude  float64 ...
Data variables:
    speed      (time) float32 ...
    direction  (time) float32 ...
Attributes:
    institution:  UBC EOAS, DFO IOS, VFPA
<xarray.DataArray 'time' ()>
array('2017-11-20T22:00:00.000000000', dtype='datetime64[ns]')
Coordinates:
    latitude   float64 ...
    longitude  float64 ...
<xarray.DataArray 'time' ()>
array('2017-11-20T23:00:00.000000000', dtype='datetime64[ns]')
Coordinates:
    latitude   float64 ...
    longitude  float64 ...
<xarray.DataArray 'speed' (time: 31)>
array([0.4, 0.5, 0.5, 0.5, 0.6, 0.6, 0.6, 0.6, 0.7, 0.7, 0.7, 0.8, 0.8, 0.8,
       0.9, 0.9, 1. , 1. , 1. , 1.2, 1.1, 1.1, 1.2, 1.1, 1.1, 1.2, 1.2, 1.2,
       1.3, 1.3, 1.3], dtype=float32)
Coordinates:
  * time       (time) datetime64[ns] 2017-11-20T22:00:00 ... 2017-11-20T23:00:00
    latitude   float64 ...
  [...]<jupyter_text>Confirm that we can append "new" observations to the initial file:<jupyter_code>for hr in arrow.Arrow.range('hour', start=utc_start_hr.shift(hours=+1), end=utc_start_hr.shift(days=+0, hours=+4)):
    with xarray.open_dataset(f'/tmp/{netcdf_filename}') as stored_obs:
        csv_filename = f'{hr.format(datetime_fmt)}Z-{hr.shift(hours=+1).format(datetime_fmt)}Z.csv'
        new_obs = csv_to_dataset(ais_path/csv_filename)
        new_obs['time'] -= pandas.to_timedelta(utc_offset)
        new_ds.coords['latitude'] = 49.2938889
        new_ds.coords['longitude'] = -123.0247222
        for var in attrs:
            new_obs[var].attrs = attrs[var]
        new_ds.attrs = dataset_attrs
        extended_obs = xarray.concat((stored_obs, new_obs), dim='time')
    # Drop repeated times because early files contain hh:00 to hh+1:00 instead of ending at hh:58
    _, index = numpy.unique(extended_obs.time.values, return_index=True)
    extended_obs = extended_obs.isel(time=index)
    extended_obs.to_netcdf(f'/tmp/{netcdf_filename}', mode='w', encoding=encoding, unlimited_dims='time')
with xarray.open_dataset(f'/tmp/{netcdf_filename}') as currents:
    print(currents)
    print(currents.time.min())
    print(currents.time.max())
    print(currents.speed)
    print(currents.direction)
    currents.speed.plot(marker='o', linestyle='None')<jupyter_output><xarray.Dataset>
Dimensions:    (time: 151)
Coordinates:
  * time       (time) datetime64[ns] 2017-11-20T22:00:00 ... 2017-11-21T03:00:00
Data variables:
    speed      (time) float32 ...
    direction  (time) float32 ...
<xarray.DataArray 'time' ()>
array('2017-11-20T22:00:00.000000000', dtype='datetime64[ns]')
<xarray.DataArray 'time' ()>
array('2017-11-21T03:00:00.000000000', dtype='datetime64[ns]')
<xarray.DataArray 'speed' (time: 151)>
array([0.4, 0.5, 0.5, 0.5, 0.6, 0.6, 0.6, 0.6, 0.7, 0.7, 0.7, 0.8, 0.8, 0.8,
       0.9, 0.9, 1. , 1. , 1. , 1.2, 1.1, 1.1, 1.2, 1.1, 1.1, 1.2, 1.2, 1.2,
       1.3, 1.3, 1.3, 1.3, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2,
       1.2, 1.2, 1.1, 1.1, 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 0.8,
       0.9, 0.9, 0.9, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.7, 0.7, 0.7,
       0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.6, 0.6,
       0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.5, 0.5, 0.6, 0.4, 0.4, 0.4,
       0.4, 0.4,[...]
{'fraction_non_alphanumeric': 0.16204337488291665, 'fraction_numerical': 0.07399668564017581, 'mean_word_length': 4.2023988005997, 'pattern_counts': {'":': 0, '<': 66, '<?xml version=': 0, '>': 66, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 6, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '51818', 'license': 'Apache-2.0', 'n_tokens_mistral': 5989, 'n_tokens_neox': 5328, 'n_words': 1045, 'path': 'notebooks/2ndNarrowsHADPCtoERDDAP.ipynb', 'repo_name': 'SalishSeaCast/analysis-doug', 'size': '13879'}