Calculate maximum and minimum total SSH for any region on the globe

In [1]:

    

from pylab import *
import netCDF4
import time
from mpl_toolkits.basemap import interp


    

In [2]:

    

# Specify region to extract [lon_min, lon_max, lat_min, lat_max]
# box = [-71.5+360, -63+360, 39.5, 46.0]  # gulf of maine
box = [-75.+360,    -60.+360., 12.,  22.] # puerto rico region
box = [-76.170675+360., -62.976517+360., 33.004829, 41.477867]


    

In [3]:

    

def amp_phase(url):
    nc = netCDF4.Dataset(url).variables
    lon = nc['lon_z'][:]
    lat = nc['lat_z'][:]
    bi = (lon>=box[0]) & (lon<=box[1])
    bj = (lat>=box[2]) & (lat<=box[3])
    zi = nc['hIm'][bi,bj].T    # index dimensions in NetCDF file are reversed from usual [j,i] ordering, so transpose
    zr = nc['hRe'][bi,bj].T
    lon = nc['lon_z'][bi]
    lat = nc['lat_z'][bj]
    # Compute amplitude and phase from real and imaginary parts
    amp = abs(zr+1j*zi)/1000.   # convert millibars to meters              
    phase = (arctan2(-zi,zr)/pi*180)  
    amp = ma.masked_equal(amp,0.)  # mask where amp = 0
    phase = ma.masked_equal(phase,0.)  # mask where phase = 0 
    return lon,lat,amp,phase


    

In [4]:

    

# OPeNDAP Data URLs for TPX08 Data on Geoport THREDDS Server
# cons=['m2','n2','s2','k2','k1','o1','p1','q1','m4']
cons=['m2','n2','s2','k2','k1','o1','p1','q1']
amp={}
for con in cons:
    url='http://geoport.whoi.edu/thredds/dodsC/usgs/data2/rsignell/data/tpx08/nc4/%s_nc4.nc' % con
    print con
    lon1,lat1,amp[con],pha = amp_phase(url)


    

    




m2
n2
s2
k2
k1
o1
p1
q1

In [5]:

    

# add up all the tidal constituents to get the maximum possible amplitude tide
# this is conservative, as this is never likely to happen 
amp_sum=zeros_like(amp['m2'])
for con in cons:
    amp_sum = amp_sum + amp[con]


    

In [6]:

    

# plot up the maximum tidal amplitude
figure(figsize=(20,8))
subplot(111,aspect=1.0/cos(lat1.mean() * pi / 180.0))
pcolormesh(lon1,lat1,amp_sum)
axis(box)
colorbar()
grid()
title('TPX08: maximum tidal height (m)');


    

In [7]:

    

print amp_sum.min()
print amp_sum.max()


    

    




0.102439201204
1.33002871655

In [8]:

    

#http://username%password@host
url='http://aviso-users:grid2010@opendap.aviso.oceanobs.com/thredds/dodsC/dataset-duacs-dt-upd-global-merged-msla-h-daily'
nc = netCDF4.Dataset(url).variables


    

In [9]:

    

nc.keys()


    

    
Out[9]:





[u'LatLon',
 u'NbLatitudes',
 u'NbLongitudes',
 u'LatLonMin',
 u'LatLonStep',
 u'time',
 u'Grid_0001']

In [10]:

    

nc['NbLongitudes'][0:5]


    

    
Out[10]:





array([ 0.        ,  0.33333333,  0.66666667,  1.        ,  1.33333333])

In [11]:

    

lon = nc['NbLongitudes'][:]
lat = nc['NbLatitudes'][:]
bi = (lon>=box[0]) & (lon<=box[1])
bj = (lat>=box[2]) & (lat<=box[3])
lon = nc['NbLongitudes'][bi]
lat = nc['NbLatitudes'][bj]


    

In [12]:

    

# load all 6 hour time steps (6672 of them!) of AVISO SSH from 1992 to 2011
time0 = time.time()
z = nc['Grid_0001'][:,bi,bj].T    # index dimensions in NetCDF file are reversed from usual [j,i] ordering, so transpose
print 'elapsed time %d seconds' % int(time.time() - time0)


    

    




elapsed time 282 seconds

In [13]:

    

times=nc['time']
dates = netCDF4.num2date(times[:],units=times.units)


    

In [14]:

    

dates[0]


    

    
Out[14]:





datetime.datetime(1992, 10, 14, 0, 0)

In [15]:

    

dates[-1]


    

    
Out[15]:





datetime.datetime(2011, 1, 19, 0, 0)

In [16]:

    

shape(z)


    

    
Out[16]:





(32, 40, 6672)

In [17]:

    

# convert cm to m
z = z/100.


    

In [18]:

    

# find the standard deviation of the AVISO SSH 
zstd=std(z,axis=2)


    

In [19]:

    

# find the minimum
zmin=z.min(axis=2)


    

In [20]:

    

# find the maximum
zmax=z.max(axis=2);


    

In [21]:

    

loni,lati = meshgrid(lon1,lat1)
zmin1=interp(zmin,lon,lat,loni,lati)
zmax1=interp(zmax,lon,lat,loni,lati)


    

In [22]:

    

# plot the max SSH
figure(figsize=(20,8))
subplot(111,aspect=1.0/cos(lat.mean() * pi / 180.0))
pcolormesh(lon1,lat1,zmax1)
colorbar()
title('maximum sub-tidal SSH (m): 1992-2011')


    

    
Out[22]:





<matplotlib.text.Text at 0x369b590>






    

In [23]:

    

# plot the min SSH
figure(figsize=(20,8))
subplot(111,aspect=1.0/cos(lat.mean() * pi / 180.0))
pcolormesh(lon1,lat1,zmin1)
colorbar()
title('Minimum sub-tidal SSH (m): 1992-2011')


    

    
Out[23]:





<matplotlib.text.Text at 0x4d9f310>






    

In [24]:

    

# the maximum/mimimum total SSH (subtidal + tide)
ssh_max = amp_sum + zmax1
ssh_min = -amp_sum + zmin1


    

In [25]:

    

# plot the maximum total SSH (subtidal + tide)
figure(figsize=(20,8))
subplot(111,aspect=1.0/cos(lat.mean() * pi / 180.0))
pcolormesh(lon1,lat1,ssh_max)
colorbar()
title('Maximum total SSH (m): 1992-2011')


    

    
Out[25]:





<matplotlib.text.Text at 0x4893e10>






    

In [26]:

    

# plot the mimimum total SSH (subtidal + tide)
figure(figsize=(20,8))
subplot(111,aspect=1.0/cos(lat.mean() * pi / 180.0))
pcolormesh(lon1,lat1,ssh_min)
colorbar()
title('Minimum total SSH (m): 1992-2011')


    

    
Out[26]:





<matplotlib.text.Text at 0x7b4b490>






    

In [27]:

    

# maximum total SSH in entire region
print ssh_max.flatten().max()


    

    




2.01138434433

In [28]:

    

# minimum total SSH in entire region
print ssh_min.flatten().min()


    

    




-2.30636098524

In [28]: