CDIP Station Plot - Period Rose

  • Based on specified CDIP buoy
  • Data from CDIP THREDDS server (archived) and justdar (most recent)

Set Station and Dates

User specifies Station number and timeframe

In [1]:
stn = '198'
startdate = "11/01/2013"
enddate = "11/30/2013"

Import Libraries

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

from itertools import groupby

Open URL as NetCDF

Paste in URL for appropriate buoy from CDIP THREDDS server

In [3]:
# URL - Archive THREDDS
url = '' + stn + 'p1/' + stn + ''

# URL - Realtime THREDDS (Use if Archived URL returns error message)
# url = '' + stn + ''

Open URL as NetCDF file

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

Assign variable names to NetCDF objects

In [5]:
ncTime = nc.variables['waveTime'][:]
Fq = nc.variables['waveFrequency'] # Assign variable name - Wave Frequency
Tp = nc.variables['waveTp'] # Assign variable name - Peak Wave Period
Dp = nc.variables['waveDp'] # Assign variable name - Peak Wave Direction
In [6]:
buoyname = nc.variables['metaStationName'][:]
buoytitle = " ".join(buoyname[:-40])

month_name = calendar.month_name[int(startdate[0:2])]
year_num = (startdate[6:10])

Timestamp Conversions

Define functions to convert between UNIX timestamps and human dates

In [7]:
# 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

Convert user-set dates to UNIX timestamps

In [8]:
unixstart = getUnixTimestamp(startdate,"%m/%d/%Y") 
neareststart = find_nearest(ncTime, unixstart)  # Find the closest unix timestamp
nearIndex = numpy.where(ncTime==neareststart)[0][0]  # Grab the index number of found date

unixend = getUnixTimestamp(enddate,"%m/%d/%Y")
future = find_nearest(ncTime, unixend)  # Find the closest unix timestamp
endIndex = numpy.where(ncTime==future)[0][0]  # Grab the index number of found date

StartDate = getHumanTimestamp(ncTime[nearIndex],"%m/%d/%Y") 

Calculate Peak Period (Tp) frequencies - 16 degree bins

Calculate total length of 'Wave Peak Direction (Dp)' variable

In [9]:
Dptot = len(Dp[nearIndex:endIndex])

Create array of total len of each of the 16 compass-direction bins (aka (number of Dp values in each 22.5-degree bin) / (total number of Dp values in time-chunk))

In [10]:
Dplens = []

dirstart = 11.25
dircount = 22.5
dirend = 348.75

for i in range(len(Dp)):
    Dpcut = (len([i for i, x in enumerate(Dp[nearIndex:endIndex]) if (dirstart < x <= (dirstart+dircount))]))/float(Dptot)
    dirstart = dirstart+dircount
    if dirstart == dirend:
Dplens.append((len([i for i, x in enumerate(Dp[nearIndex:endIndex]) if (348.75 < x or x <= (11.25))]))/float(Dptot))  # Manually append
                                                            # last bin (which crosses 360/0, and cannot be accounted for in for-loop)

Grab maximum Dp length to set range of r-axis (frequency of occurences)

In [11]:
maxlen = max(Dplens)
lenround = round(maxlen,1)

Call array of Tp (peak period) values for specified time-chunk

In [12]:
Tpcut = Tp[nearIndex:endIndex]

Create arrays of index values corresponding to each 22.5-degree bin

  • Index numbers of all Dp values in time-chunk array within each 22.5-degree bin
  • Values of Tp objects that correspond to each Dp value (via array index number)
In [13]:
Dpidxs = []
Tpvals = []

idxstart = 11.25
idxcount = 22.5
idxend = 348.75

for idx in range(len(Dp)):
    idxvals = [idx for idx, x in enumerate(Dp[nearIndex:endIndex]) if (idxstart < x <= (idxstart+idxcount))]
    Tpval = Tpcut[idxvals]
    idxstart = idxstart+idxcount
    if idxstart == idxend:
val15 = [idx for idx, x in enumerate(Dp[nearIndex:endIndex]) if (348.75 < x or x <= (11.25))]
Tp15 = Tpcut[val15]

Tp sub-categories within degree-bins

Separate Tp values within a degree-bin into period-specific categories

In [14]:
Tpcats2 = []

catstart = 0
catstep = 2

for ibin2 in range(len(Tpvals)):
    catslist2 = []
    for val in Tpvals[ibin2]:
        if catstart <= val < (catstart+catstep):
            catval = 1 # Insert '1' for Tp values < 2, to avoid confusion with later step of inserting '0' as spaceholders for non-existent values
            while catstart < val:
                catstart = catstart+catstep
            catval = catstart-2

        catstart = 0

Sort all values within each Tp category ('Tpcats' array) from smallest to largest value

  • Create list of all sorted category values within each row, to be counted
In [15]:
Tpcatsorttot = []

for row in range(len(Tpcats2)):
    sortrow = sorted(Tpcats2[row]) # sort all Tp category valus from smallest to largest

Count number of objects in each category within each row (using 'Tpcatsorttot)

In [16]:
Tpcatscount = []
Tpcatvals = []

for cat in range(len(Tpcatsorttot)):
    catlen = [len(list(group)) for key, group in groupby(Tpcatsorttot[cat])]
    catval = [key for key, group in groupby(Tpcatsorttot[cat])]

Create master list of expected categorical bin values (every 2 counts, but replacing 0 with 1. E.g. categories = 1,2,4,6,8,10...)

  • Compare 'Tpcatvals' list to master list
  • Insert 0 values into missing locations in Tpcatvals (presence of each expected length-2 categorical bin)
In [17]:
master = range(2,26,2) # Create list of values from 2-26, every 2 steps
master.insert(0,1) # add '1' as first value in the array

for val in range(len(Tpcatvals)):
    for i in range(len(master)):
        if master[i] in Tpcatvals[val]: # compare Tpcatsvals value to master list. If values match up, move on to next element
            Tpcatvals[val].insert(i,0) # If Tpcatsval value does not match master list value, insert '0' placeholder value

Create copy of master list to set colorbar range

  • Remove '1' bin category value and replace with '0' (to create '0-2 ft' range, rather than '1-2 ft' for colorbar)
In [18]:
master2 = master

Compare 'Tpcatscount' list to Tpcatvals

  • Insert 0 values into missing locations in Tpcatscount (number of counts of each category instance)
In [19]:
for c in range(len(Tpcatscount)):
    for e in range(len(Tpcatvals[c])):
        if Tpcatvals[c][e] == 0:

Category-specific array

  • Iterate over each element in a list, across all 16 lists
In [20]:
catslist = []

for el in range(len(Tpcatscount[2])):
    arr = []
    for count3 in range(len(Tpcatscount)):
        var = Tpcatscount[count3][el]

Call 'catslist' list as array, for plotting

In [21]:
catsarray = np.asarray(catslist)

Divide each element in 'catsarray' by Dptot (total count-length for specified time-period) to create 'frequencies of occurrences' array

In [22]:
Dpfloat = float(Dptot)
catsfreq = catsarray/Dpfloat

Get number of instances of Tp range (e.g. 0-2 ft) across each 22.5-degree bin

  • Find length of row in Tpcatscount
    • If length < total number of count-categories (e.g. len=2), then
In [23]:
bottom = np.cumsum(catsarray, axis=0) # Create array to call "bottom" array when plotting stacked barplot (below)

Get number of instances of Tp range (e.g. 0-2 ft) across each 22.5-degree bin in 'Frequencies of Occurrencs' array

In [24]:
bottom2 = np.cumsum(catsfreq, axis=0)

Dp Directional bins

  • Create variable of 16 mid-point compass directions
  • Convert directions from degrees to radians
In [25]:
radconvert = np.pi/180

Dpdiralt = ([18,40.5,63,85.5,108,130.5,153,175.5,198,220.5,243,265.5,288,310.5,333,355.5]) # Specify radial locations to correspond to mid-barplots
Dpradalt = [d*radconvert for d in Dpdiralt]

Plot Period Rose

  • Create Polar plot using Degree bings data ('Dpradalt') as thetas, and 'Dplens' as radii
  • Plot bottom row of barplot data using catsfreq[0] - 'Frequencies of Occurrence' array
  • Use for-loop to loop through successive rows of 'catsfreq' array
In [26]:
# Create figure and specify size
fig = plt.figure(figsize=(10,10))

# Set radial/degree variables
thetas = Dpradalt
radii = Dplens[:]

# Specify barplot parameters
widths = 0.17
cmap = plt.get_cmap('gist_ncar',12)
cmaplist = [cmap(i2) for i2 in range(cmap.N)]

# Create subplot for barplot
ax = plt.subplot(111, polar=True) # Specify polar axes
ax.set_theta_direction(-1) # Reverse direction of degrees (CW)
ax.set_theta_zero_location("N") # Specify 0-degrees as North

## Plot Tp count data
# Plot bottom row of barplot data using 'Frequencies of Occurrence' array ('catsfreq'), then loop through successive rows of 'catsarray'
bars =, catsfreq[0], width=widths, color=cmaplist[0])

for j in xrange(1, catsfreq.shape[0]):, catsfreq[j], width=widths, color=cmaplist[j], bottom=bottom2[j-1])
    if j == ((catsfreq.shape[0])-2):

# Set plot titles
plt.suptitle(buoytitle + "\n" + "Stn. " + stn + " - Period Rose", fontsize=30, y=1.12)
title(startdate + '-' + enddate, fontsize=24, y=1.08)

# Set degree gridlines to plot every 22.5 degrees
degrange = arange(0,360,22.5)
lines, labels = plt.thetagrids(degrange, labels=None, fontsize=14)

# Create colorbar
bounds = master2
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
ax2 = fig.add_axes([0.1, 0.025, 0.8, 0.03])
cb = mpl.colorbar.ColorbarBase(ax2, cmap=cmap, norm=norm, spacing='proportional', ticks=bounds, orientation='horizontal',boundaries=bounds, format='%1i')
cb.set_label('Peak Period (s)', fontsize=20)