-
Notifications
You must be signed in to change notification settings - Fork 16
/
IOOS_Models-Sediment.py
179 lines (138 loc) · 4.99 KB
/
IOOS_Models-Sediment.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
# -*- coding: utf-8 -*-
# <nbformat>3.0</nbformat>
# <headingcell level=1>
# Using Iris to access data from US-IOOS models
# <codecell>
import datetime as dt
import time
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
# <markdowncell>
# Note: `iris` is not a default package in Wakari or Anaconda, [but installation is easy](https://github.com/ioos/conda-recipes/issues/11).
# <codecell>
# use contraints to select nearest time
mytime=dt.datetime(2015,1,29,17) #specified time...
#mytime=dt.datetime.utcnow() # .... or now
print mytime
# <codecell>
# [lon_min, lon_max, lat_min, lat_max] for subsetting and plotting
bbox=[-72,-63,39,46]
# <codecell>
import iris
# <codecell>
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
# <codecell>
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]
# <codecell>
def minmax(v):
return np.min(v), np.max(v)
# <codecell>
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
# <codecell>
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));
# <codecell>
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)
# <codecell>
myplot(slice,model=ds_name,bbox=bbox,ptype='log10',crange=[-3,-1])
# <codecell>
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)
# <codecell>
myplot(slice,model=model,bbox=bbox,ptype='log10',crange=[-12,-1])
# <codecell>
print slice
# <codecell>