Enter CDIP Buoy Station Number and start date for plot
stn = '071'
startdate = "10/17/2012 16:00"
Libraries are Python packages that contain a subset of specific functions
import netCDF4
import numpy as np
import matplotlib.pyplot as plt
import datetime as dt
import urllib
import time
import calendar
Import THREDDS NetCDF file to grab Frequency and Time variables
urlarc = 'http://thredds.cdip.ucsd.edu/thredds/dodsC/cdip/archive/' + stn + 'p1/' + stn + 'p1_historic.nc'
nc = netCDF4.Dataset(urlarc)
Extract waveTime variable to find nearest buoy time to user-input time
timevar = nc.variables['waveTime'][:]
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.
# 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
Change user-inputted date to datestring that can be attached to MEM URL
unixtimestamp = getUnixTimestamp(startdate, "%m/%d/%Y %H:%M")
unix_nearest = find_nearest(timevar, unixtimestamp) # Find the closest unix timestamp
neardate = getHumanTimestamp(unix_nearest, '%Y%m%d%H%M') # Convert unix timestamp to string format to attach to URL
Import data as URL from CDIP MEM data access page
url = 'http://cdip.ucsd.edu/data_access/MEM_2dspectra.cdip?sp' + stn + '01' + neardate
Open data as text file and split into individual rows
data = urllib.urlopen(url)
readdata = data.read() # Read text file of recent data
datas = readdata.split("\n") # Split text file into individual rows
Create a new array
datas2 = []
for item in datas:
line = item.strip().split()
datas2.append(line)
datas2[0].remove('<pre>')
datas2[64].remove('</pre>')
datas2 = filter(None, datas2)
Change 'datas2' list to 'Ed_array' array
Edarray = np.asarray(datas2, dtype=object)
Ednew = np.append(Edarray,Edarray[:,0:1],1)
Create Degrees ('Dmean') and Frequency ('Fq') variables to use in plotting Energy Density
# Fq = np.asarray(np.arange(0.03,0.59,0.01)) # Use this Frequency range option when calling the 'even' option for a station
Fq = nc.variables['waveFrequency']
Dmean = np.asarray(np.arange(2.5,367.5,5))
# Convert 'Dmean' from degrees to radians
rad_convert = (np.pi/180)
Dmean_rad = [d*rad_convert for d in Dmean]
Edmax = np.amax(Ednew)
Edfloat = float(Edmax)
fig = plt.figure(figsize=(11,11))
radii = Fq[0:35] # Only call frequency bands > 0.3 Hz
thetas = Dmean_rad[:]
## Color-scale
# contour_levels = np.arange(0.00,0.1,0.0001) # Manually set colorbar min/max values
contour_levels = np.arange(Edfloat/1000,Edfloat/2,0.0001) # Automatically set colorbar min/max values based on 'Ed'
ax = plt.subplot(111, polar=True)
ax.set_theta_direction(-1)
ax.set_theta_zero_location("N")
ylabels = ([20,10,6.7,5,4])
ax.set_yticklabels(ylabels)
colorax = ax.contourf(thetas, radii, Ednew[0:35],contour_levels)
## Set titles and colorbar
suptitle('Polar Spectrum at Stn. ' + stn, fontsize=22, y=0.95, x=0.45)
title(startdate, fontsize=16, y=1.11)
cbar = fig.colorbar(colorax)
cbar.set_label('Energy Density (m*m/Hz/deg)', rotation=270, fontsize=16)
cbar.ax.get_yaxis().labelpad = 30
degrange = range(0,360,30)
lines, labels = plt.thetagrids(degrange, labels=None, frac = 1.07)
Set smallest frequency (longest period) on outer edge, and largest frequency (smallest period) on inner edge
def mapr(radiir):
"""Remap the radial axis."""
return 0.579998 - radiir
figr = plt.figure(figsize=(11,11))
radiir = Fq[0:35]
thetasr = Dmean_rad
## Color-scale
# contour_levelsr = np.arange(0,0.3,0.0001) # Manually set colorbar min/max values
contour_levelsr = np.arange(Edfloat/1000,Edfloat/2,0.0001) # Automatically set colorbar min/max values based on 'Ed'
axr = plt.subplot(111, polar=True)
axr.set_theta_direction(-1)
axr.set_theta_zero_location("N")
axr.set_yticks(arange(0.33,0.58,0.05))
axr.set_yticklabels(ylabels[::-1])
coloraxr = axr.contourf(thetasr, mapr(radiir), Ednew[0:35], contour_levelsr)
## Set titles and colorbar
suptitle('Polar Spectrum at Stn. ' + stn, fontsize=22, y=0.95, x=0.45)
title(startdate, fontsize=16, y=1.11)
cbarr = figr.colorbar(coloraxr)
cbarr.set_label('Energy Density (m*m/Hz/deg)', rotation=270, fontsize=16)
cbarr.ax.get_yaxis().labelpad = 30
degranger = range(0,360,30)
lines, labels = plt.thetagrids(degranger, labels=None, frac = 1.07)