In [1]:

    

# In ArcMap, with the multidimensional supplemental tools installed, go to: 
# Geoprocessing=>Python  
# and type these commands into the python command window


    

In [2]:

    

import netCDF4
url='http://www.smast.umassd.edu:8080/thredds/dodsC/fvcom/hindcasts/30yr_gom3/mean'
nc = netCDF4.Dataset(url)


    

In [3]:

    

nc.variables.keys()


    

    
Out[3]:





[u'lat',
 u'latc',
 u'lon',
 u'lonc',
 u'z0b',
 u'nprocs',
 u'partition',
 u'x',
 u'y',
 u'xc',
 u'yc',
 u'siglay',
 u'siglev',
 u'h',
 u'nv',
 u'nbe',
 u'ntsn',
 u'nbsn',
 u'ntve',
 u'nbve',
 u'a1u',
 u'a2u',
 u'aw0',
 u'awx',
 u'awy',
 u'art2',
 u'art1',
 u'iint',
 u'time',
 u'Itime',
 u'Itime2',
 u'Times',
 u'zeta',
 u'file_date',
 u'u',
 u'v',
 u'omega',
 u'ww',
 u'ua',
 u'va',
 u'temp',
 u'salinity',
 u'km',
 u'kh',
 u'kq',
 u'q2',
 u'q2l',
 u'l',
 u'short_wave',
 u'net_heat_flux',
 u'uwind_stress',
 u'vwind_stress',
 u'fvcom_mesh']

In [4]:

    

temp = nc.variables['temp']
print temp


    

    




<type 'netCDF4.Variable'>
float32 temp(u'time', u'siglay', u'node')
    long_name: temperature
    standard_name: sea_water_temperature
    units: degrees_C
    coordinates: time siglay lat lon
    type: data
    mesh: fvcom_mesh
    location: node
unlimited dimensions = ()
current size = (396, 45, 48451)

In [5]:

    

timevar = nc.variables['time']
print timevar


    

    




<type 'netCDF4.Variable'>
float32 time(u'time',)
    long_name: time
    units: days since 1858-11-17 00:00:00
    format: modified julian day (MJD)
    time_zone: UTC
unlimited dimensions = ()
current size = (396,)

In [6]:

    

# get first (0) and last (-1) time values and convert to datestamps
start = netCDF4.num2date(timevar[0],timevar.units)
stop = netCDF4.num2date(timevar[-1],timevar.units)

print('start: %s' % start.strftime('%Y-%b-%d %H:%M:%S'))
print(' stop: %s' % stop.strftime('%Y-%b-%d %H:%M:%S'))


    

    




start: 1978-Jan-16 12:00:00
 stop: 2010-Dec-16 10:01:53

In [7]:

    

siglay = nc.variables['siglay']  # fraction of total depth
h=nc.variables['h']    # water depth
print siglay
print h


    

    




<type 'netCDF4.Variable'>
float32 siglay(u'siglay', u'node')
    long_name: Sigma Layers
    standard_name: ocean_sigma_coordinate
    positive: up
    valid_min: -1.0
    valid_max: 0.0
    formula_terms: sigma: siglay eta: zeta depth: h
unlimited dimensions = ()
current size = (45, 48451)

<type 'netCDF4.Variable'>
float32 h(u'node',)
    long_name: Bathymetry
    standard_name: sea_floor_depth_below_geoid
    units: m
    coordinates: lat lon
    type: data
    mesh: fvcom_mesh
    location: node
unlimited dimensions = ()
current size = (48451,)

In [8]:

    

# z(meters) of layer 20 at node 2000
node = 2000
lev = 20
z = siglay[lev,node]*h[node]
print z


    

    




-851.425

In [9]:

    

# plot sigma values at specified node
import matplotlib.pyplot as plt
plt.plot(siglay[:,node],'k-o')
plt.grid()
plt.xlabel('sigma level')
plt.ylabel('fraction of total depth')
plt.show()


    

In [9]: