In [1]:

    

from pylab import *
import netCDF4


    

In [2]:

    

url = 'http://testbedapps-dev.sura.org/thredds/dodsC/alldata/Shelf_Hypoxia/tamu/roms/tamu_roms.nc'
#url = 'http://tds.ve.ismar.cnr.it:8080/thredds/dodsC/ismar/model/field2/run1/Field_2km_1108_out30min_his_0724.nc'
#url='http://tds.ve.ismar.cnr.it:8080/thredds/dodsC/field2_test/run1/his'
#url='http://ourocean.jpl.nasa.gov:8080/thredds/dodsC/SCBfcst/scb_latest_fcst_roms.nc'
#url='http://geoport.whoi.edu/thredds/dodsC/examples/bora_feb.nc'
#####################################################################################

nc = netCDF4.Dataset(url)
mask = nc.variables['mask_rho'][:]

# read longitude, latitude
lon_rho = nc.variables['lon_rho'][:]
lat_rho = nc.variables['lat_rho'][:]

# read water depth
depth = nc.variables['h'][:,:]


'''
z(n,k,j,i) = eta(n,j,i)*(1+s(k)) + depth_c*s(k) +
             (depth(j,i)-depth_c)*C(k)

  C(k) = (1-b)*sinh(a*s(k))/sinh(a) + 
         b*[tanh(a*(s(k)+0.5))/(2*tanh(0.5*a)) - 0.5]

formula_terms: s: s_rho eta: zeta depth: h a: theta_s b: theta_b depth_c: hc
'''


    

    
Out[2]:





'\nz(n,k,j,i) = eta(n,j,i)*(1+s(k)) + depth_c*s(k) +\n             (depth(j,i)-depth_c)*C(k)\n\n  C(k) = (1-b)*sinh(a*s(k))/sinh(a) + \n         b*[tanh(a*(s(k)+0.5))/(2*tanh(0.5*a)) - 0.5]\n\nformula_terms: s: s_rho eta: zeta depth: h a: theta_s b: theta_b depth_c: hc\n'

In [3]:

    

s = nc.variables['s_rho'][:]
a = nc.variables['theta_s'][:]
b = nc.variables['theta_b'][:]
depth_c = nc.variables['hc'][:]

C = (1-b)*sinh(a*s)/sinh(a) + b*[tanh(a*(s+0.5))/(2*tanh(0.5*a)) - 0.5]


    

    




---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-3-19806ecf7bc0> in <module>()
      1 s = nc.variables['s_rho'][:]
----> 2 a = nc.variables['theta_s'][:]
      3 b = nc.variables['theta_b'][:]
      4 depth_c = nc.variables['hc'][:]
      5 

KeyError: 'theta_s'

In [7]:

    

nc.variables.keys()
print nc.variables['s_rho']


    

    




<type 'netCDF4.Variable'>
float64 s_rho(u's_rho',)
    long_name: S-coordinate at RHO-points
    valid_min: -1.0
    valid_max: 0.0
    positive: up
    standard_name: ocean_s_coordinate_g1
    formula_terms: s: s_rho C: Cs_r eta: zeta depth: h depth_c: hc
    field: s_rho, scalar
    _CoordinateTransformType: Vertical
    _CoordinateAxes: s_rho
    _CoordinateAxisType: GeoZ
    _CoordinateZisPositive: up
unlimited dimensions = ()
current size = (20,)

In [4]:

    

# Reshape 1D vertical variables so we can broadcast
C.shape = (np.size(C), 1, 1)
s.shape = (np.size(s), 1, 1)


    

    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-4-6acb5541b674> in <module>()
      1 # Reshape 1D vertical variables so we can broadcast
----> 2 C.shape = (np.size(C), 1, 1)
      3 s.shape = (np.size(s), 1, 1)

NameError: name 'C' is not defined

In [7]:

    

tidx = -1       # just get the final time step, for now.
# read a 3D temperature field at specified time step
temp = nc.variables['temp'][tidx, :, :, :]
# read a 2D water level (height of ocean surface) at specified time step
eta = nc.variables['zeta'][tidx, :, :]
# calculate the 3D field of z values (vertical coordinate) at this time step
z = eta*(1+s) + depth_c*s + (depth-depth_c)*C


    

In [18]:

    

shape(s)


    

    
Out[18]:





(20, 1, 1)

In [9]:

    

# plot surface temperature
pcolormesh(lon_rho,lat_rho,temp[-1,:,:])


    

    
Out[9]:





<matplotlib.collections.QuadMesh at 0x31ad190>






    

In [10]:

    

shape(temp)


    

    
Out[10]:





(20, 60, 160)

In [11]:

    

shape(z)


    

    
Out[11]:





(20, 60, 160)

In [12]:

    

shape(lon_rho)


    

    
Out[12]:





(60, 160)

In [13]:

    

# are there tools in EPD to plot an isosurface of this 
# curvilinear, stretched vertical coordinate temperature field?


    

In [19]:

    

lon3d=lon_rho*ones(20).shape(20,1,1)


    

    




---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-19-d56be457d3fe> in <module>()
----> 1 lon3d=lon_rho*ones(20).shape(20,1,1)

TypeError: 'tuple' object is not callable

In [20]:

    

onez=ones(20)
onez.shape=(20,1,1)


    

In [21]:

    

shape(onez)


    

    
Out[21]:





(20, 1, 1)

In [22]:

    

lon3d = lon_rho*onez


    

In [23]:

    

shape(lon3d)


    

    
Out[23]:





(20, 60, 160)

In [24]:

    

lon3d[:,30,20]


    

    
Out[24]:





array([ 12.86812732,  12.86812732,  12.86812732,  12.86812732,
        12.86812732,  12.86812732,  12.86812732,  12.86812732,
        12.86812732,  12.86812732,  12.86812732,  12.86812732,
        12.86812732,  12.86812732,  12.86812732,  12.86812732,
        12.86812732,  12.86812732,  12.86812732,  12.86812732])

In [ ]: