CDIP Station Plot - Polar Spectrum (Energy Density) - CDIP MEM 2D Spectrum

  • Based on specified CDIP buoy
  • Polar plot of Energy Spectrum, based on Frequency bands (64) and Directions (0-360)
  • Data from CDIP Maximum Entropy Method (MEM) 2D Energy Spectrum (http://cdip.ucsd.edu/data_access/MEM_2dspectra.cdip)
  • CDIP buoys record data every 30 minutes (48 times/day)

Specify station and date

Enter CDIP Buoy Station Number and start date for plot

In [1]:
stn = '071'
startdate = "10/17/2012 16:00"

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 as dt
import urllib
import time
import calendar

Data Access

THREDDS NetCDF

Import THREDDS NetCDF file to grab Frequency and Time variables

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

In [4]:
timevar = nc.variables['waveTime'][:]

Convert to UNIX timestamp to find nearest date

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

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

Data - MEM 2D Spectrum

Import data as URL from CDIP MEM data access page

In [7]:
url = 'http://cdip.ucsd.edu/data_access/MEM_2dspectra.cdip?sp' + stn + '01' + neardate

Open data as text file and split into individual rows

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

  • Split lines into individual objects
  • Remove 'pre' characters and filter array for empty objects
In [9]:
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

  • Append the first colum of 'Ed_array' to the end of 'Ed_array' so that dimensions are (64L,73L), and the last column connects back to the first column. This allows the polar plot to call all 72 (5-degree bin) columns of real data when plotting the colormesh.
In [10]:
Edarray = np.asarray(datas2, dtype=object)
Ednew = np.append(Edarray,Edarray[:,0:1],1)

Dmean and Fq Variables

Create Degrees ('Dmean') and Frequency ('Fq') variables to use in plotting Energy Density

  • Extract Frequency variable from NetCDF, to use as 'Fq' array
  • Create manually-specified Degrees ('Dmean') list, of 5-degree intervals from 0-360
In [11]:
# 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']
In [12]:
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]
In [13]:
Edmax = np.amax(Ednew)
Edfloat = float(Edmax)

Plot Polar Spectrum

  • Specify radial and angular axes of polar plot
    • 'radii' = Fq array
    • 'thetas' = Dmean_rad array
  • Define colormesh range
  • Change 'thetas' alignment so that 0 Degrees is at the top
  • Apply a colormesh (contourf) of Energy Density data to polar plot
  • Add colorbar and titles
In [14]:
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)

Plot Polar Spectrum - inverted Frequency scale

Set smallest frequency (longest period) on outer edge, and largest frequency (smallest period) on inner edge

In [15]:
def mapr(radiir):
   """Remap the radial axis."""
   return 0.579998 - radiir
In [16]:
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)