import datetime as dt
import time
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
Note: iris
is not a default package in Wakari or Anaconda, but installation is easy.
# use contraints to select nearest time
mytime=dt.datetime(2015,1,29,17) #specified time...
#mytime=dt.datetime.utcnow() # .... or now
print mytime
2015-01-29 17:00:00
# [lon_min, lon_max, lat_min, lat_max] for subsetting and plotting
bbox=[-72,-63,39,46]
import iris
def find_timevar(cube):
"""Return the variable attached to
time axis and rename it to time."""
try:
cube.coord(axis='T').rename('time')
print('Renaming {} to time'.format(cube.coord('time').var_name))
except:
pass
timevar = cube.coord('time')
return timevar
def time_near(cube, start):
"""Return the nearest time to `start`.
TODO: Adapt to the new slice syntax"""
timevar = find_timevar(cube)
try:
time1 = timevar.units.date2num(start)
itime = timevar.nearest_neighbour_index(time1)
except IndexError:
itime = -1
return timevar.points[itime]
def minmax(v):
return np.min(v), np.max(v)
def var_lev_date(url=None,var=None,mytime=None,lev=0,subsample=1,bbox=None):
'''
specify lev=None if the variable does not have layers
'''
time0= time.time()
# cube = iris.load_cube(url,iris.Constraint(name=var.strip()))[0]
cube = iris.load_cube(url,var)
# cube = iris.load(url,var)[0]
# print cube.coord('time')
try:
cube.coord(axis='T').rename('time')
except:
pass
slice = cube.extract(iris.Constraint(time=time_near(cube,mytime)))
if bbox is None:
imin=0
jmin=0
imax=-2
jmax=-2
else:
lats=slice.coord(axis='Y').points
lons=slice.coord(axis='X').points
inregion = np.logical_and(np.logical_and(lons > bbox[0], lons < bbox[1]),
np.logical_and(lats > bbox[2], lats < bbox[3]))
# extract the rectangular subarray containing the whole valid region ...
region_inds = np.where(inregion)
jmin, jmax = minmax(region_inds[0])
imin, imax = minmax(region_inds[1])
if lev is None:
slice = slice[jmin:jmax+1:subsample,imin:imax+1:subsample]
else:
slice = slice[lev,jmin:jmax+1:subsample,imin:imax+1:subsample]
print 'slice retrieved in %f seconds' % (time.time()-time0)
return slice
def myplot(slice,model=None,crange=None,bbox=None,ptype='linear'):
"""
bbox = [lon_min, lon_max, lat_min, lat_max]
crange = [vmin, vmax]
ptype = 'linear' or 'log10'
model = dataset name (string), e.g. 'COAWST US-EAST'
"""
fig=plt.figure(figsize=(12,8))
lat=slice.coord(axis='Y').points
lon=slice.coord(axis='X').points
time=slice.coord('time')[0]
ax1 = plt.subplot(111,aspect=(1.0/cos(mean(lat)*pi/180.0)))
if ptype is 'linear':
data = ma.masked_invalid(slice.data)
elif ptype is 'log10':
data = np.log10(ma.masked_invalid(slice.data))
else:
raise
if crange is None:
im = ax1.pcolormesh(lon,lat,data);
else:
im = ax1.pcolormesh(lon,lat,data,vmin=crange[0],vmax=crange[-1]);
if bbox is not None:
ax1.set_xlim(bbox[:2])
ax1.set_ylim(bbox[2:])
fig.colorbar(im,ax=ax1)
ax1.grid()
date=time.units.num2date(time.points)
date_str=date[0].strftime('%Y-%m-%d %H:%M:%S %Z')
plt.title('%s: %s: %s' % (model,slice.long_name,date_str));
ds_name='VIIRS Data'
url='http://www.star.nesdis.noaa.gov/thredds//dodsC/CoastWatch/VIIRS/Rrs672/Daily/NE05'
var='Remote sensing reflectance at 672 nm'
lev=0
slice=var_lev_date(url=url,var=var, mytime=mytime, lev=lev, subsample=2,bbox=bbox)
Renaming None to time slice retrieved in 12.877612 seconds
myplot(slice,model=ds_name,bbox=bbox,ptype='log10',crange=[-3,-1])
/home/usgs/anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.py:1039: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal if aspect == 'normal': /home/usgs/anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.py:1044: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal elif aspect in ('equal', 'auto'): -c:17: RuntimeWarning: divide by zero encountered in log10
model='USGS/COAWST Model'
url='http://geoport.whoi.edu/thredds/dodsC/coawst_4/use/fmrc/coawst_4_use_best.ncd'
var='suspended noncohesive sediment, size class 06'
lev=-1
#slice=var_lev_date(url=url,var=var, mytime=mytime, lev=lev, subsample=1,bbox=[-72,-63,39,46])
# if we issue the above command to slice out data in a specified bounding box,
# it actually takes longer (34 seconds) because we have a curvlinear grids
# and the whole 2D lon,lat fields must be downloaded to determine the index range to extract
# of course if were doing this again, we could just use the already calculated index ranges, and
# it would be faster. But since we just need this one slice, we here just download the whole thing
slice=var_lev_date(url=url,var=var, mytime=mytime, lev=lev, subsample=1)
/home/usgs/anaconda/lib/python2.7/site-packages/iris/fileformats/cf.py:1038: UserWarning: Ignoring formula terms variable u'h' referenced by data variable u'v' via variable u's_rho': Dimensions (u'eta_rho', u'xi_rho') do not span (u'time', u's_rho', u'eta_v', u'xi_v') warnings.warn(msg) /home/usgs/anaconda/lib/python2.7/site-packages/iris/fileformats/cf.py:1038: UserWarning: Ignoring formula terms variable u'zeta' referenced by data variable u'v' via variable u's_rho': Dimensions (u'time', u'eta_rho', u'xi_rho') do not span (u'time', u's_rho', u'eta_v', u'xi_v') warnings.warn(msg) /home/usgs/anaconda/lib/python2.7/site-packages/iris/fileformats/cf.py:1038: UserWarning: Ignoring formula terms variable u'h' referenced by data variable u'u' via variable u's_rho': Dimensions (u'eta_rho', u'xi_rho') do not span (u'time', u's_rho', u'eta_u', u'xi_u') warnings.warn(msg) /home/usgs/anaconda/lib/python2.7/site-packages/iris/fileformats/cf.py:1038: UserWarning: Ignoring formula terms variable u'zeta' referenced by data variable u'u' via variable u's_rho': Dimensions (u'time', u'eta_rho', u'xi_rho') do not span (u'time', u's_rho', u'eta_u', u'xi_u') warnings.warn(msg) /home/usgs/anaconda/lib/python2.7/site-packages/iris/fileformats/_pyke_rules/compiled_krb/fc_rules_cf_fc.py:1291: UserWarning: Gracefully filling 'time' dimension coordinate masked points warnings.warn(msg.format(str(cf_coord_var.cf_name)))
slice retrieved in 13.558762 seconds
myplot(slice,model=model,bbox=bbox,ptype='log10',crange=[-12,-1])
-c:17: RuntimeWarning: invalid value encountered in log10 /home/usgs/anaconda/lib/python2.7/site-packages/numpy/ma/core.py:797: RuntimeWarning: invalid value encountered in less_equal return umath.less_equal(x, self.critical_value)
print slice
suspended noncohesive sediment, size class 06 / (kilogram meter-3) (-- : 335; -- : 895) Auxiliary coordinates: latitude x x longitude x x sea_floor_depth x x water_surface_height_above_reference_datum x x Scalar coordinates: S-coordinate at RHO-points: -0.03125 S-coordinate parameter, critical depth: 5.0 meter S-coordinate stretching curves at RHO-points: -0.00225488671278 forecast_reference_time: 2015-01-29 13:00:00 time: 2015-01-29 17:00:00 Attributes: CPP_options: COAWST, ANA_BPFLUX, ANA_BSFLUX, ANA_BTFLUX, ANA_FSOBC, ANA_M2OBC, ANA_NUDGCOEF,... Conventions: CF-1.4, _Coordinates DODS_EXTRA.Unlimited_Dimension: ocean_time NLM_LBC: EDGE: WEST SOUTH EAST NORTH zeta: Clo Cha ... _ChunkSize: [ 1 4 112 299] _CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention ana_file: ROMS/Functionals/ana_btflux.h, ROMS/Functionals/ana_fsobc.h, ROMS/Functionals/ana_m2obc.h,... bry_file_01: ./forcings3/USE_coawst_bdy.nc cdm_data_type: GRID clm_file_01: ./forcings3/USE_coawst_clm.nc code_dir: /raid1/jcwarner/Models/COAWST_forecast compiler_command: /usr/local/mpi/bin/mpif90 compiler_flags: -fastsse -Mipa=fast -tp k8-64 -Mfree compiler_system: pgi cpu: x86_64 featureType: GRID field: sand_06, scalar, series file: ./Output/coawst_use_his_0001.nc format: netCDF-3 64bit offset file frc_file_01: ./forcings/cloud_coawst.nc frc_file_02: ./forcings/Pair_coawst.nc frc_file_03: ./forcings/Qair_coawst.nc frc_file_04: ./forcings/rain_coawst.nc frc_file_05: ./forcings/swrad_coawst.nc frc_file_06: ./forcings/Tair_coawst.nc frc_file_07: ./forcings/lwrad_coawst.nc frc_file_08: ./forcings/Big_grd2_coawstwind.nc frc_file_09: ./grids/tide_forc_USeast_grd16_osu_rev2.nc grd_file: ./grids/USeast_grd19.nc header_dir: /raid1/jcwarner/Models/COAWST_forecast/Projects/coawst header_file: coawst.h his_base: ./Output/coawst_use_his history: ROMS/TOMS, Version 3.7, Thursday - February 5, 2015 - 3:12:36 AM ; FMRC... ini_file: ./restart/ocean_use_rst.nc location: Proto fmrc:coawst_4_use os: Linux rst_file: ./Output/ocean_use_rst.nc script_file: size_class: 3.1250E-02 millimeter summary: Experimental forecast model product from the USGS Coupled Ocean Atmosphere... svn_rev: exported svn_url: https:://myroms.org/svn/src tiling: 010x004 time: ocean_time title: COAWST Forecast System : USGS : US East Coast and Gulf of Mexico (Expe... type: ROMS/TOMS history file