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 numpy.ma 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 = 'http://thredds.cdip.ucsd.edu/thredds/dodsC/cdip/archive/' + stn + 'p1/' + stn + 'p1_historic.nc'
# URL - Realtime THREDDS (Use if Archived URL returns error message)
# url = 'http://thredds.cdip.ucsd.edu/thredds/dodsC/cdip/realtime/' + stn + 'p1_rt.nc'
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
Dplens.append(Dpcut)
if dirstart == dirend:
break
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
Dpidxs.append(idxvals)
Tpvals.append(Tpval)
if idxstart == idxend:
break
val15 = [idx for idx, x in enumerate(Dp[nearIndex:endIndex]) if (348.75 < x or x <= (11.25))]
Tp15 = Tpcut[val15]
Tpvals.append(Tp15)
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
catslist2.append(catval)
else:
while catstart < val:
catstart = catstart+catstep
catval = catstart-2
catslist2.append(catval)
catstart = 0
Tpcats2.append(catslist2)
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
Tpcatsorttot.append(sortrow)
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])]
Tpcatscount.append(catlen)
Tpcatvals.append(catval)
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
continue
else:
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
master2.remove(1)
master2.insert(0,0)
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:
Tpcatscount[c].insert(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]
arr.append(var)
catslist.append(arr)
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 = ax.bar(thetas, catsfreq[0], width=widths, color=cmaplist[0])
for j in xrange(1, catsfreq.shape[0]):
ax.bar(thetas, catsfreq[j], width=widths, color=cmaplist[j], bottom=bottom2[j-1])
if j == ((catsfreq.shape[0])-2):
break
# 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.ax.tick_params(labelsize=14)
cb.set_label('Peak Period (s)', fontsize=20)
