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 http://cdip.ucsd.edu/themes/cdip?pb=1&bl=cdip?pb=1&d2=p1&d2=p70&u3=s:093:st:1:v:xy_plot for example Timeseries plots

NOTE:

  • Data for Directional Displacement (XYZ) are only stored in individual Archived Deployment NetCDF files, NOT in the long-term 'historic.nc' 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: http://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. 'd41') as .nc files.
  • You can either select a specific deployment and copy its OPeNDAP URL, or insert the deployment number into the 'deploy' variable below.

Station Number

Enter CDIP Buoy Station Number, Deployment Number, Date and Plot Duration to plot averages

  • Option to set QC flag level, based on the level of data you want to plot
    • QC Flag Categories: [1,2,3,4,9] = [good,not_evaluated,questionable,bad,missing]
In [1]:
stn = '067'
deploy = '17' # Set deployment number from .nc file

startdate = "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 Libraries

Libraries are Python packages that contain a subset of specific functions

In [2]:
import netCDF4
import numpy as np
import matplotlib.pyplot as plt
import datetime
import matplotlib.dates as mpldts
import time
import calendar

Data Access

  • Data are stored on the CDIP THREDDS server as NetCDF files

Import Data as NetCDF from THREDDS URL

XYZ Displacement data from CDIP buoys are available only via deployment-specific archived NetCDF files.

In [3]:
# CDIP Archived Individual Deployment Dataset URL
data_url = 'http://thredds.cdip.ucsd.edu/thredds/dodsC/cdip/archive/' + stn + 'p1/' + stn + 'p1_d' + deploy + '.nc'

Open Remote Dataset from CDIP THREDDS Server

In [4]:
nc = netCDF4.Dataset(data_url)

Read Buoy Variables

Assign variable names to the relevant variables from NetCDF file:

In [5]:
ncTime = nc.variables['waveTime'][:]
xdisp = nc.variables['xyzXDisplacement'] # Make a numpy array of three directional displacement variables (x, y, z)
ydisp = nc.variables['xyzYDisplacement']
zdisp = nc.variables['xyzZDisplacement']
qcflag = nc.variables['xyzFlagPrimary']
starttime = nc.variables['xyzStartTime'][:] # Variable that gives start time for buoy data collection
samplerate = nc.variables['xyzSampleRate'][:] # Variable that gives rate (frequency, Hz) of sampling

# Specify variables to automatically grab station name and number
buoyname = nc.variables['metaStationName'][:]
buoytitle = " ".join(buoyname[:-42])

Local Indexing Functions

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

In [6]:
# Find nearest value in numpy array
def find_nearest(array,value):
    idx = (np.abs(array-value)).argmin()
    return array[idx]

# Convert to unix timestamp
def getUnixTimestamp(humanTime,dateFormat):
    unixTimestamp = int(calendar.timegm(datetime.datetime.strptime(humanTime, dateFormat).timetuple()))
    return unixTimestamp

# Convert to human readable timestamp
def getHumanTimestamp(unixTimestamp, dateFormat):
    humanTimestamp = datetime.datetime.utcfromtimestamp(int(unixTimestamp)).strftime(dateFormat)
    return humanTimestamp

Check Deployment Start/End Dates

View the Start and End dates of the specified deployment file, in order to set user-specified Start/End dates above

In [7]:
StartDate = getHumanTimestamp(ncTime[0],"%m/%d/%Y")
EndDate = getHumanTimestamp(ncTime[-1],"%m/%d/%Y")
print StartDate
print EndDate
04/18/2008
06/19/2009

Time Index Values

  • Find the UNIX values that correspond to specified Timeseries Start and End times.
  • Create a special time array using buoy Sampling Start Time ('starttime' variable), 'unixend' as end UNIX stamp, and (1/samplerate[0]) as spacing period.
  • Find the index number in 'sampletime' of the unix start/end times ('unixstart'/'unixend')
In [8]:
# Find UNIX timestamps for human-formatted start/end dates
unixstart = getUnixTimestamp(startdate,"%m/%d/%Y %H:%M") 
unixend = unixstart + (duration * 60) # Create UNIX end stamp by adding 30 minutes (30 min * 60 sec) to 'unixstart'

# Create specialized array using UNIX Start and End times minus Filter Delay, and Sampling Period (1/samplerate) 
# to calculate sub-second time values that correspond to Z-Displacement sampling values
sampletime = np.arange((starttime),(ncTime[-1]),(1/(samplerate)))

# Find corresponding start/end date index numbers in 'Sampletime' array    
startindex = sampletime.searchsorted(unixstart) 
endindex = sampletime.searchsorted(unixend)

Call variable that only encompasses specified timerange

In [9]:
sampletimecut = sampletime[startindex:endindex]

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

  • Multiply sampletimecut*1000 to remove sub-second decimal place and incorporate millisecond steps
  • Call 'sampletimecut' datetime objects as '[ms]'
  • Convert to datetime objects for plotting
In [10]:
sampletimecut *= 1000
sampledtcutms = sampletimecut.astype('datetime64[ms]').astype(datetime.datetime)

Plot data on graphs

  • Plot three directional displacements as separate subplots, using specified subset timerange
    • 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
In [11]:
# Specify figure size
fig = plt.figure(figsize=(15,15))

# Limit data to date/times
x = xdisp[startindex:endindex]
y = ydisp[startindex:endindex]
z = zdisp[startindex:endindex]
qc = qcflag[startindex:endindex]

# 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)
ppltxdisp = plt.subplot(3,1,1)
ppltxdisp.plot(sampledtcutms,x,'b')
ppltydisp = plt.subplot(3,1,2, sharex=ppltxdisp)
ppltydisp.plot(sampledtcutms,y,'b')
ppltzdisp = plt.subplot(3,1,3, sharex=ppltxdisp)
ppltzdisp.plot(sampledtcutms,z,'b')

# Set titles
plt.suptitle(buoytitle + "\n" + "Time: " + startdate, 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=floor(min(min(x), min(y), min(z)))
ymax=ceil(max(max(x), max(y), max(z)))

# Set y-axis limits for each plot
ppltxdisp.set_ylim(ymin,ymax)
ppltydisp.set_ylim(ymin,ymax)
ppltzdisp.set_ylim(ymin,ymax)

# Label each subplot title
ppltxdisp.set_title('North/South Displacement', fontsize=14)
ppltydisp.set_title('East/West Displacement', fontsize=14)
ppltzdisp.set_title('Vertical Displacement', fontsize=14,y=1)

# Label each y-axis
ppltxdisp.set_ylabel('X Displacement (m)', fontsize=14)
ppltydisp.set_ylabel('Y Displacement (m)', fontsize=14)
ppltzdisp.set_ylabel('Z Displacement (m)', fontsize=14)

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

# Plot dashed gridlines
ppltxdisp.grid(axis='y', which='major', color='b', linestyle='-', alpha=0.25)
ppltydisp.grid(axis='y', which='major', color='b', linestyle='-', alpha=0.25)
ppltzdisp.grid(axis='y', which='major', color='b', linestyle='-', alpha=0.25)