Year-long Avg of Ocean Temperature - CDIP buoys

  • Based on specified CDIP buoy
  • Data from CDIP THREDDS server (archived) and justdar (most recent)
  • CDIP buoys record data every 30 minutes (48 times/day)

User label changes

In [1]:
# Label which CDIP buoy is being used
stn_num = '093'
Buoy = 'CDIP Mission Bay, CA'


# Paste in URL for appropriate buoy from CDIP THREDDS server
urlarc = 'http://thredds.cdip.ucsd.edu/thredds/dodsC/cdip/archive/' + stn_num + 'p1/' + stn_num + 'p1_historic.nc'


# URL to most recent CDIP temperature data via justdar
today = datetime.datetime.today()
todaynew = datetime.datetime.strptime(str(today),"%Y-%m-%d %H:%M:%S.%f").strftime('%Y%m%d')
urlcurr= 'http://cdip.ucsd.edu/data_access/justdar.cdip?093+dt+20140709+' + todaynew

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 urllib
import pandas as pd
import datetime
import time

from numpy import *
from time import mktime

Open URL as NetCDF

In [3]:
# Archived CDIP THREDDS data - through 7/9/2014
nc = netCDF4.Dataset(urlarc)

Assign Variables - archived data

In [4]:
ncTime = nc.variables['sstTime'][:]
timeall = [datetime.datetime.fromtimestamp(t) for t in ncTime] # Convert ncTime variable to datetime stamps
npTime = np.asarray(ncTime)  # Make a numpy array of ncTime variable
lat = nc.variables['gpsLatitude']
lon = nc.variables['gpsLongitude']
nctemp = nc.variables['sstSeaSurfaceTemperature']

Buoy variables - recent data from justdar

In [5]:
data = urllib.urlopen(urlcurr) # Assign variable to text file of recent data from justdar
readdata = data.read() # Read text file of recent data

datas = readdata.split("\n") # Split text file into individual rows

datas.remove('</pre>')
In [6]:
# Create a new data array of individual data rows

datas2 = []

for item in datas:
     line = item.strip().split()
     datas2.append(line) 
In [7]:
# Create data arrays for the 'time' and 'temperature' variables from recent data

timecurr = []
tempcurr = []
count = 0

for each in datas2:
    timevar = datas2[count][0]
    tempvar = datas2[count][1]
    count = count+1
    timecurr.append(timevar)
    tempcurr.append(tempvar)
In [8]:
## Array for newest Temp values

# Create an array of strings of the temperature data values
tempcurrarr = np.array(tempcurr)
tempcurrarr = tempcurrarr[2:]

# Change array of strings to array of floats, for plotting
tempcurrfloat = tempcurrarr.astype(np.float)


## Array for newest Time values

# Create an array of strings of time data values
timecurrarr = np.array(timecurr)
timecurrarr = timecurrarr[2:]

# Change array of strings to array of floats, for plotting
timecurrfloat = timecurrarr.astype(np.float)

Change UNIX timestamps to human dates

In [9]:
# Define 3 timestamp functions

# 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(time.mktime(datetime.datetime.strptime(humanTime, dateFormat).timetuple()))
    return unixTimestamp

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


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

daySteps = 365*48 #Number of measurements per day (48) multiplied by year-days (365)
In [10]:
# Check Start and End dates of buoy data
StartDate = getHumanTimestamp(ncTime[0],"%m/%d/%Y")
EndDate = getHumanTimestamp(ncTime[-1],"%m/%d/%Y")
print StartDate
print EndDate
10/23/1997
07/09/2014

Create individual year-long variables

In [11]:
yearsall = [datetime.datetime.fromtimestamp(t).year for t in ncTime] # extract 'year' components from each timestamp
years = list(set(yearsall)) # Organize all years into chronological list, and remove duplicates
print years
[1997, 1998, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014]

In [12]:
from time import mktime

# Create empty arrays for index numbers for 'start' and 'end' of each year
timeindex_start = []
timeindex_end = []
yearcount = 0


for yeari in years:
    ys = "01/05/" + str(years[yearcount]) # Get unix timestamp of Jan. 1 and Dec. 31 for each year
    ye = "12/28/" + str(years[yearcount])

    unixstart = getUnixTimestamp(ys,"%m/%d/%Y")
    nearest_date = find_nearest(npTime, unixstart)  # Find the closest unix timestamp
    near_index = numpy.where(npTime==nearest_date)[0][0]  # Grab the index number of found date
    unixend = getUnixTimestamp(ye,"%m/%d/%Y")
    future_date = find_nearest(npTime, unixend)  # Find the closest unix timestamp
    future_index = numpy.where(npTime==future_date)[0][0]  # Grab the index number of found date

    start = getHumanTimestamp(ncTime[near_index],"%m/%d/%Y")
    end = getHumanTimestamp(ncTime[future_index],"%m/%d/%Y")
    yearcount = yearcount+1
    timeindex_start.append(near_index)
    timeindex_end.append(future_index)
In [13]:
# Combine 'start' and 'end' indices into one time index
timeindex = timeindex_start + timeindex_end
timeindexsort = sorted(timeindex, key=int)
In [14]:
timeallcop = []

suball = 0

for each in timeall:
    eachdate = timeall[suball]
    timecop = eachdate.replace(year=2008)
    suball = suball+1
    timeallcop.append(timecop)

2014 Continuing Temp

In [15]:
# Create a masked array of most recent time data (from justdar), so that it is shifted forward to match up with archived data

timecurrdate = []

sub14c = 0

for timevar in timecurrarr:
    timecurrsub = timecurrarr[sub14c]
    timenewform = datetime.datetime.strptime(timecurrsub, '%Y%m%d%H%M')
    sub14c = sub14c+1
    timecurrdate.append(timenewform)
In [16]:
## Create new time array to replace all year dates with '2008' - this allows all temp-series to be plotted on same graph axes.
timecopcurr = []

subcurr = 0

for each in timecurrdate:
    datecurr = timecurrdate[subcurr]
    timecurrcop = datecurr.replace(year=2008)
    subcurr = subcurr+1
    timecopcurr.append(timecurrcop)
In [17]:
time2014place = numpy.zeros(shape=(11752))

 # This number comes from 17520-'time2013', aka the total expected length of a year-long temp dataset (17520) minus the length of the archived 2013 dataset. Set this number so that you don't include 2014 data in the 2013 temp series
time2014comb = np.hstack([time2014place,timecopcurr])
time2014comb = ma.masked_where(time2014comb==0,time2014comb)
time2014combstr = time2014comb.tolist()
In [18]:
# Create a masked array of most recent temp data so that it is shifted forward to match up with archived data
# NOTE: this array does NOT combine with the archived data from 2014 in 'temp2014'

# Set length of masked array to length of temp2014, aka archived data from beginning of year (so that current data will be shifted forward that many steps)
temp2014place = numpy.zeros(shape=(11752))

# This number comes from 17520-'temp2013', aka the total expected length of a year-long temp dataset (17520) minus the length of the archived 2014 dataset. 
temp2014comb = np.hstack([temp2014place,tempcurrfloat])
temp2014comb = ma.masked_where(temp2014comb==0,temp2014comb)
temp2014combstr = str(temp2014comb)

Plot Temperature Time-Series

In [19]:
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111)
i = 0

for plotindex in range(len(timeindexsort)):
    pyear = ax.plot(timeallcop[timeindexsort[i]:timeindexsort[i+1]], nctemp[timeindexsort[i]:timeindexsort[i+1]], 'k', alpha=0.15)
    i = i+2    
    if i == 22:
      break
        
# Add current data manually 
p14curr, = ax.plot(time2014combstr,temp2014comb, 'r')

# Create manual variable for legend
p97, = ax.plot(timeallcop[timeindexsort[0]:timeindexsort[1]], nctemp[timeindexsort[0]:timeindexsort[1]], 'k', alpha=0.15)
     
    
yearrange = str(years[0]) + "-" + str(years[-2])

plt.ylim(0,30)
ax.set_xticklabels(['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'])  
title('Daily SST at ' + Buoy + ' ' + yearrange, fontsize=30)
plt.ylabel(u'Temperature (°C)', fontsize=18)
plt.xlabel('Month', fontsize = 16)
ax.legend((p97,p14curr,),(yearrange,"2014",))
Out[19]:
<matplotlib.legend.Legend at 0x10e27c150>
In [20]:
fig.savefig('CDIP/MB_SST.png')
---------------------------------------------------------------------------
IOError                                   Traceback (most recent call last)
<ipython-input-20-d3e1bea0c9f0> in <module>()
----> 1 fig.savefig('CDIP/MB_SST.png')

/usr/local/lib/python2.7/site-packages/matplotlib/figure.pyc in savefig(self, *args, **kwargs)
   1468             self.set_frameon(frameon)
   1469 
-> 1470         self.canvas.print_figure(*args, **kwargs)
   1471 
   1472         if frameon:

/usr/local/lib/python2.7/site-packages/matplotlib/backend_bases.pyc in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, **kwargs)
   2190                 orientation=orientation,
   2191                 bbox_inches_restore=_bbox_inches_restore,
-> 2192                 **kwargs)
   2193         finally:
   2194             if bbox_inches and restore_bbox:

/usr/local/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in print_png(self, filename_or_obj, *args, **kwargs)
    516         renderer.dpi = self.figure.dpi
    517         if is_string_like(filename_or_obj):
--> 518             filename_or_obj = open(filename_or_obj, 'wb')
    519             close = True
    520         else:

IOError: [Errno 2] No such file or directory: 'CDIP/MB_SST.png'
In []:
 
Window size: x
Viewport size: x