In [59]:

    

# import required packages
import netCDF4
#import matplotlib, pylab # already in pylab
import openearthtools.io.netcdf as nc


    

In [60]:

    

# define file, use r'name' to be able to  use windows-style slashes in name
ncfile = r'd:\DATA\Resmon4\nc\MERIS2WAQ_2003_2011_tsm_dd_grid_TIM_log_1.nc'


    

In [61]:

    

# open file
grid = netCDF4.Dataset(ncfile)
# see what variables are in there
grid.variables


    

    
Out[61]:





OrderedDict([(u'time', <netCDF4.Variable object at 0x04490F30>), (u'crs', <netCDF4.Variable object at 0x044B78A0>), (u'lon', <netCDF4.Variable object at 0x044B79E0>), (u'lat', <netCDF4.Variable object at 0x044B7990>), (u'lonbounds', <netCDF4.Variable object at 0x044B7A30>), (u'latbounds', <netCDF4.Variable object at 0x044B7A80>), (u'SPM', <netCDF4.Variable object at 0x044B7AD0>), (u'mask', <netCDF4.Variable object at 0x044B7B20>)])

In [62]:

    

# load coordinates
# the [:,:] are required as grid.variables['var'] is only a pointer
lon = grid.variables['lon'][:,:]
lat = grid.variables['lat'][:,:]


    

In [63]:

    

# load the timesteps
tim = netCDF4.num2date(grid.variables['time'], units=grid.variables['time'].units)
# see what times are in there; tim[1,-1] ???
tim


    

    
Out[63]:





array([datetime.datetime(2003, 1, 3, 12, 0),
       datetime.datetime(2003, 1, 4, 12, 0),
       datetime.datetime(2003, 1, 5, 12, 0), ...,
       datetime.datetime(2011, 12, 19, 12, 0),
       datetime.datetime(2011, 12, 22, 11, 0),
       datetime.datetime(2011, 12, 27, 12, 0)], dtype=object)

In [64]:

    

# see how many times are in there, and select one
grid.variables['SPM'].shape


    

    
Out[64]:





(2320, 235, 243)

In [65]:

    

# pick a time and load that slice
it  = 56
val = grid.variables['SPM'][it,:,:]


    

In [68]:

    

# load coastline annotation from opendap.deltares.nl
ldbfile = r'd:\opendap.deltares.nl\thredds\dodsC\opendap\deltares\landboundaries\northsea.nc'
tmp = netCDF4.Dataset(ldbfile); tmp.variables; L = {}
L['lon'] = tmp.variables['lon'][:]
L['lat'] = tmp.variables['lat'][:]


    

In [69]:

    

matplotlib.pyplot.pcolor(lon,lat,10**val, vmin=0, vmax=1.5) # pcolormesh won't work here, it interprets data at centersmatplotlib.pyplot.colorbar()
matplotlib.pyplot.axis('tight')
matplotlib.pyplot.axis([-2,9,50,57])
matplotlib.pyplot.title(str(tim[it]))
gca().set_aspect(1./cos(mean(ylim())/180*pi),adjustable='box') # see Matlab axislat()
colorbar()
plot(L['lon'],L['lat'],'k')


    

    
Out[69]:





[<matplotlib.lines.Line2D at 0x8a63f90>]






    

In [ ]: