CDIP Station Plot - Timeseries of Directional Displacement (XYZ)

Create subplots of: North/South Displacement (X), East/West Displacement (Y) and Vertical Displacement (Z) using OPeNDAP service from CDIP THREDDS Server.

• See https://cdip.ucsd.edu/themes/cdip?d2=p70&u2=s:067:v:xz_plot:dt:200805210900 for example Timeseries plot

Archive datasets:

• Data for Directional Displacement (XYZ) from previous deployments are only stored in individual Archived Deployment NetCDF files, NOT in the long-term 'historic.nc' file for each station (data are collected on a sub-second level, and are too large to store as one long-term dataset).
• In order to access the Displacement variables, you must set the THREDDS URL to a specific deployment number. Deployment numbers and corresponding dates vary by station.
• All archived data for CDIP stations can be found at this URL: https://thredds.cdip.ucsd.edu/thredds/catalog/cdip/archive/catalog.html. Clicking on an individual station will give you a list of all of its deployments (e.g. 'd17') as .nc files.
• Insert the deployment number into the 'deploy' variable below.

Realtime datasets:

• Realtime data for CDIP stations can be found at this URL: https://thredds.cdip.ucsd.edu/thredds/catalog/cdip/realtime/catalog.html.
• Realtime Directional Displaments (XYZ) data for each station are stored in NetCDF files with the naming format 'XXXp1_xy.nc'.
• The deployment number is not necessary to identify the file as this is the current deployment.
• If a recent deployment has ended (maintenance, offsite or decommisioned) there will not be a realtime dataset. Look to the archive mentioned above for data.

Station Number

Enter CDIP Buoy Station Number, realtime or archive, Deployment Number (archive only), Date and Plot Duration (minutes)

Option to set QC flag level based on the level of qc data you want to plot

• QC Flag Categories: [1,2,3,4,9] = [good,not_evaluated,questionable,bad,missing]

stn = '067'
dataset = 'archive' # Enter 'archive' or 'realtime'
deploy = '17' # If archive dataset, set deployment number from .nc file

start_date = '05/21/2008 09:00' # MM/DD/YYYY HH:MM
duration  = 30 # Set length of timeseries (minutes)

qc_level = 2 # Filter data with qc flags above this number 

Import Python Libraries

This notebook was generated using:

• python==3.8.5
• matplotlib==3.3.2
• netCDF4==1.5.4

import netCDF4
import numpy as np
import matplotlib.pyplot as plt
import datetime
import matplotlib.dates as mpldts
import calendar

Data Access

Data are stored on the CDIP THREDDS server in NetCDF files

Import NetCDF Data from THREDDS URL

XYZ Displacement data from CDIP buoys are available in two seperate locations:

• archive (previous deployments)
• realtime (currently active deployment)

# Archive
data_url = 'http://thredds.cdip.ucsd.edu/thredds/dodsC/cdip/archive/' + stn + 'p1/' + stn + 'p1_d' + deploy + '.nc'
# Realtime
if dataset == 'realtime':
    data_url = 'http://thredds.cdip.ucsd.edu/thredds/dodsC/cdip/realtime/' + stn + 'p1_xy.nc'

Open Remote Dataset from CDIP THREDDS Server

nc = netCDF4.Dataset(data_url)
# Turn off auto masking
nc.set_auto_mask(False)

Read Buoy Variables

Assign variable names to the relevant variables from NetCDF file:

xdisp = nc.variables['xyzXDisplacement'] # Make a numpy array of three directional displacement variables (x, y, z)
ydisp = nc.variables['xyzYDisplacement']
zdisp = nc.variables['xyzZDisplacement']
qc_flag = nc.variables['xyzFlagPrimary']
filter_delay = nc.variables['xyzFilterDelay']
start_time = nc.variables['xyzStartTime'][:] # Variable that gives start time for buoy data collection
sample_rate = nc.variables['xyzSampleRate'][:] # Variable that gives rate (frequency, Hz) of sampling
end_time = start_time + (len(xdisp)/sample_rate) # Calulate end time for buoy data collection

# Get station name and number for plot title
station_name = nc.variables['metaStationName'][:]
station_title = station_name.tobytes().decode().split('\x00',1)[0]

Helper Functions for Dates

NetCDF stores time variables as UNIX dates (values are in seconds since 01-01-1970 00:00:00). The helper functions below allow us to convert between human-format timestamps (e.g. 'MM/DD/YYYY') and UNIX timestamps.

# Convert to unix timestamp
def get_unix_timestamp(human_time,dateFormat):
    unix_timestamp = int(calendar.timegm(datetime.datetime.strptime(human_time, dateFormat).timetuple()))
    return unix_timestamp

# Convert to human readable timestamp
def get_human_timestamp(unix_timestamp, dateFormat):
    human_timestamp = datetime.datetime.utcfromtimestamp(int(unix_timestamp)).strftime(dateFormat)
    return human_timestamp

Check Start and End Dates

View the Start and End dates of all available xyz data

data_start = get_human_timestamp(start_time - filter_delay[0],"%m/%d/%Y %H:%M:%S")
data_end = get_human_timestamp(end_time - filter_delay[0],"%m/%d/%Y %H:%M:%S")
print("data_start: " + data_start)
print("  data_end: " + data_end)

data_start: 04/18/2008 02:59:14
  data_end: 06/19/2009 14:59:15

Time Index Values

• Find the UNIX timestamps for user specified Start and Duration.
• Create a special sample_time array using buoy's Sampling start_time and end_time and (1/sample_rate[0]) as spacing period. Adjust for the filter delay to get the actual displacement time.
• Find the index number in 'sample_time' of the users desired start and end times ('unix_start', 'unix_end')

# Find UNIX timestamps for user human-formatted start/end dates
unix_start = get_unix_timestamp(start_date,"%m/%d/%Y %H:%M") 
unix_end = unix_start + (duration * 60) # Create UNIX end stamp by adding duration to 'unix_start'

# Create specialized array using UNIX Start and End times minus Filter Delay, and Sampling Period (1/sample_rate) 
# to calculate sub-second time values that correspond to Z-Displacement sampling values
sample_time = np.arange((start_time - filter_delay[0]), end_time - filter_delay[0],(1/(sample_rate)))

# Find corresponding start/end date index numbers in 'sample_time' array    
start_index = sample_time.searchsorted(unix_start) 
end_index = sample_time.searchsorted(unix_end)

Make array that only encompasses specified timerange

sample_time_cut = sample_time[start_index:end_index]

Transform 'sample_time_cut' array into Datetime objects so that time component can be plotted as timestamps

• Multiply sample_time_cut*1000 to remove sub-second decimal place and incorporate millisecond steps
• Call 'sample_time_cut' datetime objects as '[ms]'
• Convert to datetime objects for plotting

sample_time_cut *= 1000
sample_t_cut_ms = sample_time_cut.astype('datetime64[ms]').astype(datetime.datetime)

Plot Directional Displacements (XYZ)

Plot three directional displacements as separate subplots

• Test x,y,z, datasets for QC Primary Flags
• Mask data that do not pass QC flag threshold (set at top)
• Specify axis ranges, labels and grids for each subplot

# Specify figure size
fig = plt.figure(figsize=(15,15))

# Limit data to date/times
x = xdisp[start_index:end_index]
y = ydisp[start_index:end_index]
z = zdisp[start_index:end_index]
qc = qc_flag[start_index:end_index]

# Filter out by quality control level
x = np.ma.masked_where(qc>qc_level,x)
y = np.ma.masked_where(qc>qc_level,y)
z = np.ma.masked_where(qc>qc_level,z)

# Create 3 stacked subplots for three Directional Displacement Parameters (xyz)
plt_x = plt.subplot(3,1,1)
plt_x.plot(sample_t_cut_ms,x,'b')
plt_y = plt.subplot(3,1,2, sharex=plt_x)
plt_y.plot(sample_t_cut_ms,y,'b')
plt_z = plt.subplot(3,1,3, sharex=plt_x)
plt_z.plot(sample_t_cut_ms,z,'g')

# Set titles
plt.suptitle(station_title + "\n" + "Time: " + start_date, fontsize=22, y=0.97)

# Set x-axis tick format to "HH:MM:SS" and tick interval to every 5 minutes
#days = mpldts.MinuteLocator(interval=5) 
daysFmt = mpldts.DateFormatter('%H:%M:%S')
#plt.gca().xaxis.set_major_locator(days)
plt.gca().xaxis.set_major_formatter(daysFmt)

ymin=np.floor(min(min(x), min(y), min(z)))
ymax=np.ceil(max(max(x), max(y), max(z)))

# Set y-axis limits for each plot
plt_x.set_ylim(ymin,ymax)
plt_y.set_ylim(ymin,ymax)
plt_z.set_ylim(ymin,ymax)

# Label each subplot title
plt_x.set_title('North/South Displacement', fontsize=14)
plt_y.set_title('East/West Displacement', fontsize=14)
plt_z.set_title('Vertical Displacement', fontsize=14,y=1)

# Label each y-axis
plt_x.set_ylabel('X Displacement (m)', fontsize=14)
plt_y.set_ylabel('Y Displacement (m)', fontsize=14)
plt_z.set_ylabel('Z Displacement (m)', fontsize=14)

# Label x-axis
plt.xlabel('Time (UTC)', fontsize=18)

# Plot dashed gridlines
plt_x.grid(axis='y', which='major', color='b', linestyle='-', alpha=0.25)
plt_y.grid(axis='y', which='major', color='b', linestyle='-', alpha=0.25)
plt_z.grid(axis='y', which='major', color='g', linestyle='-', alpha=0.25)