-
Notifications
You must be signed in to change notification settings - Fork 6
/
read_glider_data.py
481 lines (384 loc) · 15.2 KB
/
read_glider_data.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
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#%%
def read_glider_data_thredds_server(url_thredds,var_name,scatter_plot,**kwargs):
"""
Created on Tue Feb 5 10:05:37 2019
@author: aristizabal
This function reads glider data from the IOOS Data Assembly Center (DAC).
Inputs:
url_thredds: url address or directory on local computer where the netcdf
file with the glider data resides. Example:
url_glider = 'https://data.ioos.us/thredds/dodsC/deployments/rutgers/ru33-20180801T1323/ru33-20180801T1323.nc3.nc'
var_name: variable to plot. Ex: 'temperature', 'salinity'. Make sure
to use the same name as defined in the netcdf file
scatter_plot: if equal to 'yes' then a scatter plot
of the glider transect is plotted
kwargs = dict(date_ini='2018/09/01/00',date_end='2018/09/10/00')
date_ini: initial date the user wish to visualize the data.
This function uses the data format '%Y/%m/%d/%H'.
Examaple: date_ini = '2018/09/01/00'
if it is not passed date_ini is the initial time of deployment
date_end: final date the user wish to visualize the data.
This function uses the data format '%Y/%m/%d/%H'.
Examaple: date_ini = '2018/09/10/00'
if it is not passed date_ini is the initial time of deployment
Outputs:
varg: all the glider profiles for the variable chosen within the user defined time window
latg: latitude within the user defined time window
long: longitude within the user defined time window
timeg: user defined time window
depthg: depth vector for all profiles
"""
import xarray as xr
import datetime
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import cmocean
date_ini = kwargs.get('date_ini', None)
date_end = kwargs.get('date_end', None)
gdata = xr.open_dataset(url_thredds+'#fillmismatch')#,decode_times=False)
dataset_id = gdata.id.split('_')[0]
variable = np.asarray(gdata.variables[var_name][0][:])
latitude = np.asarray(gdata.latitude[0])
longitude = np.asarray(gdata.longitude[0])
depth = np.asarray(gdata.depth[0])
time = np.asarray(gdata.time[0])
time = np.asarray(mdates.num2date(mdates.date2num(time)))
#time = netCDF4.num2date(time,time.units)
#timestamp = np.asarray(gdata.time[0])
# Find time window of interest
if date_ini==None:
tti = time[0]
else:
tti = datetime.datetime.strptime(date_ini,'%Y/%m/%d/%H')
if date_end==None:
tte = time[-1]
else:
tte = datetime.datetime.strptime(date_end,'%Y/%m/%d/%H')
oktimeg = np.logical_and(mdates.date2num(time) >= mdates.date2num(tti),\
mdates.date2num(time) <= mdates.date2num(tte))
# Fields within time window
varg = variable[oktimeg,:].T
latg = latitude[oktimeg]
long = longitude[oktimeg]
depthg = depth[oktimeg,:].T
timeg = time[oktimeg]
# Scatter plot
if scatter_plot == 'yes':
if var_name == 'temperature':
color_map = cmocean.cm.thermal
clabel = var_name[0].upper()+var_name[1:] + ' ($^oC$)'
else:
if var_name == 'salinity':
color_map = cmocean.cm.haline
clabel = var_name[0].upper()+var_name[1:]
else:
color_map = 'RdBu_r'
timeg_matrix = np.tile(timeg.T,(depthg.shape[0],1))
ttg = np.ravel(timeg_matrix)
dg = np.ravel(depthg)
teg = np.ravel(varg)
kw = dict(c=teg, marker='*', edgecolor='none')
fig, ax = plt.subplots(figsize=(10, 3))
cs = ax.scatter(ttg,-dg,cmap=color_map,**kw)
ax.set_xlim(timeg[0], timeg[-1])
ax.set_ylabel('Depth (m)',fontsize=14)
cbar = plt.colorbar(cs)
cbar.ax.set_ylabel(clabel,fontsize=14)
ax.set_title(dataset_id,fontsize=16)
xfmt = mdates.DateFormatter('%H:%Mh\n%d-%b')
ax.xaxis.set_major_formatter(xfmt)
plt.ylim([-np.nanmax(dg),0])
return varg, timeg, latg, long, depthg, dataset_id
#%%
def retrieve_dataset_id_erddap_server(url_erddap,lat_lim,lon_lim,date_ini,date_end):
"""
Created on Tue Feb 5 10:05:37 2019
@author: aristizabal
This function retrieves glider ids from the IOOS
Data Assembly Center (DAC).
Inputs:
url_erddap: url address of erddap server
Example: 'https://data.ioos.us/gliders/erddap'
lat_lim: latitude limits for the search.
Example, lat_lim = [38.0,40.0]
lon_lim: longitude limits for the search.
Example, lon_lim = [-75.0,-72.0]
date_ini: initial date of time window.
This function accepts the data formats '%Y-%m-%d T %H:%M:%S Z' and '%Y/%m/%d/%H'.
Examaple: date_ini = '2018-08-02T00:00:00Z'
date_end: initial date of time window.
This function accepts the data formats '%Y-%m-%d T %H:%M:%S Z' and '%Y/%m/%d/%H'.
Examaple: date_ini = '2018-08-10T00:00:00Z'
Outputs:
gliders: list of gliders ids that fall within the lat, lon and
time constraints
"""
from erddapy import ERDDAP
import pandas as pd
e = ERDDAP(server = url_erddap)
# Search constraints
kw = {
'min_lon': lon_lim[0],
'max_lon': lon_lim[1],
'min_lat': lat_lim[0],
'max_lat': lat_lim[1],
'min_time': date_ini,
'max_time': date_end,
}
search_url = e.get_search_url(response='csv', **kw)
search = pd.read_csv(search_url)
# Extract the IDs
gliders = search['Dataset ID'].values
return gliders
#%%
def retrieve_variable_names_erddap_server(url_erddap,dataset_id):
"""
Created on Tue Nov 3 11:26:05 2020
@author: aristizabal
This function retrieves the variable names from the IOOS
and Rutgers erddapp glider servers.
Inputs:
url_erddap: url address of erddap server
Example: 'https://data.ioos.us/gliders/erddap'
dataset_id: Example: 'ng231-20190901T0000'
Outputs:
variables: list of variables for the requested dataset_id
"""
from erddapy import ERDDAP
e = ERDDAP(
server=url_erddap,
protocol='tabledap',
response='nc')
e.dataset_id = dataset_id
df = e.to_pandas()
variable_names = [var for var in df.columns]
print('List of available variables ')
print(variable_names)
return variable_names
#%%
def read_glider_variables_erddap_server(url_erddap,dataset_id,\
lat_lim,lon_lim,\
variable_names=['time'],
**kwargs):
"""
Created on Tue Nov 3 11:26:05 2020
@author: aristizabal
This function reads glider variables from the IOOS
and Rutgers erddapp glider servers.
Inputs:
url_erddap: url address of erddap server
Example: 'https://data.ioos.us/gliders/erddap'
dataset_id: Example: 'ng231-20190901T0000'
variable_names: list of variable names.
Example:
variable_names = ['depth',
'latitude',
'longitude',
'time',
'temperature',
'salinity']
The default value is variable_names=['time']
lat_lim: latitude limits for the search.
Example, lat_lim = [38.0,40.0]
lon_lim: longitude limits for the search.
Example, lon_lim = [-75.0,-72.0]
date_ini: initial date of time window.
This function accepts the data formats '%Y-%m-%d T %H:%M:%S Z' and '%Y/%m/%d/%H'.
Examaple: date_ini = '2018-08-02T00:00:00Z' or '2018/08/02/00'
date_end: initial date of time window.
This function uses the data format '%Y-%m-%d T %H:%M:%S Z'.
Examaple: date_ini = '2018-08-10T00:00:00Z' and '2018/08/10/00'
Outputs:
df: Pandas data frame with all the variables requested as vectors
"""
from erddapy import ERDDAP
import numpy as np
date_ini = kwargs.get('date_ini', None)
date_end = kwargs.get('date_end', None)
# Find time window of interest
if np.logical_or(date_ini==None,date_end==None):
constraints = {
'latitude>=': lat_lim[0],
'latitude<=': lat_lim[1],
'longitude>=': lon_lim[0],
'longitude<=': lon_lim[1],
}
else:
constraints = {
'time>=': date_ini,
'time<=': date_end,
'latitude>=': lat_lim[0],
'latitude<=': lat_lim[1],
'longitude>=': lon_lim[0],
'longitude<=': lon_lim[1],
}
e = ERDDAP(
server=url_erddap,
protocol='tabledap',
response='nc'
)
e.dataset_id = dataset_id
e.constraints = constraints
e.variables = variable_names
# Converting glider data to data frame
# Cheching that data frame has data
df = e.to_pandas()
if len(df) > 3:
df = e.to_pandas(
parse_dates=True)
return df
#%%
def read_glider_data_erddap_server(url_erddap,dataset_id,\
lat_lim,lon_lim,scatter_plot,**kwargs):
"""
Created on Tue Feb 5 10:05:37 2019
@author: aristizabal
This function reads glider data from the IOOS
Data Assembly Center (DAC).
Inputs:
url_erddap: url address of thredds server
Example: 'https://data.ioos.us/gliders/erddap'
dataset_id: this id is retrieved from the glider DAC using the
function "retrieve_glider_id_erddap_server".
Example: 'ru30-20180705T1825'
lat_lim: latitude limits for the search.
Example, lat_lim = [38.0,40.0]
lon_lim: longitude limits for the search.
Example, lon_lim = [-75.0,-72.0]
date_ini: initial date of time window.
This function accepts the data formats '%Y-%m-%d T %H:%M:%S Z' and '%Y/%m/%d/%H'.
Examaple: date_ini = '2018-08-02T00:00:00Z' or '2018/08/02/00'
date_end: initial date of time window.
This function uses the data format '%Y-%m-%d T %H:%M:%S Z'.
Examaple: date_ini = '2018-08-10T00:00:00Z' and '2018/08/10/00'
scatter_plot: if equal to 'yes' then a scatter plot
of the glider transect is plotted
Outputs:
tempg: all the glider profiles of temperature within the user defined time window
saltg: all the glider profiles of salinity within the user defined time window
latg: latitude within the user defined time window
long: longitude within the user defined time window
timeg: user defined time window
depthg: depth vector for all profiles
"""
from erddapy import ERDDAP
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import cmocean
import numpy as np
date_ini = kwargs.get('date_ini', None)
date_end = kwargs.get('date_end', None)
# Find time window of interest
if np.logical_or(date_ini==None,date_end==None):
constraints = {
'latitude>=': lat_lim[0],
'latitude<=': lat_lim[1],
'longitude>=': lon_lim[0],
'longitude<=': lon_lim[1],
}
else:
constraints = {
'time>=': date_ini,
'time<=': date_end,
'latitude>=': lat_lim[0],
'latitude<=': lat_lim[1],
'longitude>=': lon_lim[0],
'longitude<=': lon_lim[1],
}
variables = [
'depth',
'latitude',
'longitude',
'time',
'temperature',
'salinity'
]
e = ERDDAP(
server=url_erddap,
protocol='tabledap',
response='nc'
)
e.dataset_id = dataset_id
e.constraints = constraints
e.variables = variables
# Converting glider data to data frame
# Cheching that data frame has data
df = e.to_pandas()
if len(df) > 3:
df = e.to_pandas(
index_col='time (UTC)',
parse_dates=True,
skiprows=(1,) # units information can be dropped.
).dropna()
# Coverting glider vectors into arrays
timeg, ind = np.unique(df.index.values,return_index=True)
latg = df['latitude (degrees_north)'].values[ind]
long = df['longitude (degrees_east)'].values[ind]
dg = df['depth (m)'].values
vg1 = df[df.columns[3]].values
vg2 = df[df.columns[4]].values
zn = np.int(np.max(np.diff(np.hstack([ind,len(dg)]))))
depthg = np.empty((zn,len(timeg)))
depthg[:] = np.nan
tempg = np.empty((zn,len(timeg)))
tempg[:] = np.nan
saltg = np.empty((zn,len(timeg)))
saltg[:] = np.nan
for i,ii in enumerate(ind):
if i < len(timeg)-1:
depthg[0:len(dg[ind[i]:ind[i+1]]),i] = dg[ind[i]:ind[i+1]]
tempg[0:len(vg1[ind[i]:ind[i+1]]),i] = vg1[ind[i]:ind[i+1]]
saltg[0:len(vg2[ind[i]:ind[i+1]]),i] = vg2[ind[i]:ind[i+1]]
else:
depthg[0:len(dg[ind[i]:len(dg)]),i] = dg[ind[i]:len(dg)]
tempg[0:len(vg1[ind[i]:len(vg1)]),i] = vg1[ind[i]:len(vg1)]
saltg[0:len(vg2[ind[i]:len(vg2)]),i] = vg2[ind[i]:len(vg2)]
# Scatter plot
if scatter_plot == 'yes':
color_map = cmocean.cm.thermal
varg = tempg
timeg_matrix = np.tile(timeg.T,(depthg.shape[0],1))
ttg = np.ravel(timeg_matrix)
dg = np.ravel(depthg)
teg = np.ravel(varg)
kw = dict(c=teg, marker='*', edgecolor='none')
fig, ax = plt.subplots(figsize=(10, 3))
cs = ax.scatter(ttg,-dg,cmap=color_map,**kw)
#fig.colorbar(cs)
ax.set_xlim(timeg[0], timeg[-1])
ax.set_ylabel('Depth (m)',fontsize=14)
cbar = plt.colorbar(cs)
cbar.ax.set_ylabel('Temperature ($^oC$)',fontsize=14)
ax.set_title(dataset_id,fontsize=16)
xfmt = mdates.DateFormatter('%H:%Mh\n%d-%b')
ax.xaxis.set_major_formatter(xfmt)
plt.ylim([-np.nanmax(dg),0])
color_map = cmocean.cm.haline
varg = saltg
timeg_matrix = np.tile(timeg.T,(depthg.shape[0],1))
ttg = np.ravel(timeg_matrix)
dg = np.ravel(depthg)
teg = np.ravel(varg)
kw = dict(c=teg, marker='*', edgecolor='none')
fig, ax = plt.subplots(figsize=(10, 3))
cs = ax.scatter(ttg,-dg,cmap=color_map,**kw)
#fig.colorbar(cs)
ax.set_xlim(timeg[0], timeg[-1])
ax.set_ylabel('Depth (m)',fontsize=14)
cbar = plt.colorbar(cs)
cbar.ax.set_ylabel('Salinity',fontsize=14)
ax.set_title(dataset_id,fontsize=16)
xfmt = mdates.DateFormatter('%H:%Mh\n%d-%b')
ax.xaxis.set_major_formatter(xfmt)
plt.ylim([-np.nanmax(dg),0])
else:
tempg = np.nan
saltg = np.nan
timeg = np.nan
latg = np.nan
long = np.nan
depthg = np.nan
return tempg, saltg, timeg, latg, long, depthg