forked from dankelley/oce
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathargo.R
1882 lines (1838 loc) · 92.6 KB
/
argo.R
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
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# vim:textwidth=128:expandtab:shiftwidth=4:softtabstop=4
#' Class to Store Argo Data
#'
#' This class stores data from Argo floats.
#'
#' An `argo` object may be read with [read.argo()] or
#' created with [as.argo()]. Argo data can be gridded to constant
#' pressures with [argoGrid()] or subsetted with
#' [subset,argo-method()]. Plots can be made with
#' [plot,argo-method()], while [summary,argo-method()]
#' produces statistical summaries and `show` produces overviews.
#'
#' See \url{http://www.argo.ucsd.edu/Gridded_fields.html} for some
#' argo-related datasets that may be useful in a wider context.
#'
#' @templateVar class argo
#'
#' @templateVar dataExample The key items stored in this slot include equal-length vectors `time`, `longitude`, `latitude` and equal-dimension matrices `pressure`, `salinity`, and `temperature`.
#'
#' @templateVar metadataExample Examples that are of common interest include `id`, a vector of ID codes for the profiles, and `dataMode`, a vector of strings indicating whether the profile is in archived mode (`"A"`), realtime mode (`"R"`), or delayed mode (`"D"`).
#'
#' @template slot_summary
#'
#' @template slot_put
#'
#' @template slot_get
#'
#' @author Dan Kelley and Clark Richards
#'
#' @family classes provided by oce
#' @family things related to argo data
setClass("argo", contains="oce")
#' ARGO float dataset
#'
#' This holds data from ARGO 6900388 in the North Atlantic.
#'
#' To quote Argo's website: "These data were collected and made freely
#' available by the International Argo Program and the national programs
#' that contribute to it. (http//www.argo.ucsd.edu,
#' http://argo.jcommops.org). The Argo Program is part of the
#' Global Ocean Observing System."
#'
#' Below is the official citation (note that this DOI has web links for
#' downloads):
#' Argo (2017). Argo float data and metadata from Global Data Assembly Centre
#' (Argo GDAC) - Snapshot of Argo GDAC of July, 8st 2017. SEANOE.
#' \url{http://doi.org/10.17882/42182#50865}
#'
#' @name argo
#' @docType data
#'
#' @examples
#'\donttest{
#' library(oce)
#' data(argo)
#' summary(argo)
#' data(coastlineWorld)
#' plot(argo, which="trajectory")
#'}
#'
#' @source This file was downloaded using the unix command
#'\preformatted{
#' ftp ftp://ftp.ifremer.fr/ifremer/argo/dac/bodc/6900388/6900388_prof.nc
#'} issued on 2017 July 7.
#'
#' @family datasets provided with oce
#' @family things related to argo data
NULL
#' Extract Something From an Argo Object
#'
#' @param x an [argo-class] object.
#'
#' @templateVar class argo
#'
#' @section Details of the specialized `argo` method:
#'
#' There are several possibilities, depending on the nature of `i`.
#' Note that all of these calculations are done with
#' `salinityAdjusted`, if that is present, or with `salinity`
#' otherwise, and similar for temperature and pressure.
#'
#' * If `i` is `"profile"` and `j` is an integer vector,
#' then an argo object is returned, as specified by `j`. For example,
#' `argo[["profile", 2:5]]` is equivalent to
#' `subset(argo, profile \%in\% 2:5)`.
#'
#' * If `i` is `"CT"`, then
#' Conservative Temperature is returned, as computed with
#' [`gsw::gsw_CT_from_t`]`(SA,t,p)`, where
#' first `SA` is computed as explained
#' in the next item, `t` is in-situ temperature,
#' and `p` is pressure.
#'
#' * If `i` is `"N2"`, then the square of buoyancy is returned,
#' as computed with [swN2()].
#'
#' * If `i` is `"SA"`, then
#' Absolute Salinity is returned, as computed with
#' [gsw::gsw_SA_from_SP()].
#'
#' * If `i` is `"sigmaTheta"`, then
#' potential density anomaly (referenced to zero
#' pressure) is computed, with [swSigmaTheta()], where the
#' equation of state is taken to be
#' [`getOption`]`("oceEOS",default="gsw")`.
#'
#' * If `i` is `"theta"`, then
#' potential temperature (referenced to zero
#' pressure) is computed, with [swTheta()], where the
#' equation of state is taken to be
#' [`getOption`]`("oceEOS",default="gsw")`.
#'
#' * If `i` is `"depth"`, then
#' a matrix of depths is returned.
#'
#' * If `i` is in the `data` slot of `x`,
#' then it is returned, otherwise if it is in the `metadata` slot,
#' then that is returned, otherwise `NULL` is returned.
#'
#' @template sub_subTemplate
#'
#' @examples
#' data(argo)
#' # 1. show that dataset has 223 profiles, each with 56 levels
#' dim(argo[['temperature']])
#'
#' # 2. show importance of focussing on data flagged 'good'
#' fivenum(argo[["salinity"]],na.rm=TRUE)
#' fivenum(argo[["salinity"]][argo[["salinityFlag"]]==1],na.rm=TRUE)
#'
#' @family things related to argo data
#' @author Dan Kelley
setMethod(f="[[",
signature(x="argo", i="ANY", j="ANY"),
definition=function(x, i, j, ...) {
res <- NULL
if (i == "profile") {
## This assignment to profile is merely to prevent a NOTE from
## the syntax checker. It is needed because of issues with non-standard
## evaluation in subset() calls. This is a problem that many
## package authors have encountered; see e.g.
## https://stackoverflow.com/questions/23475309/in-r-is-it-possible-to-suppress-note-no-visible-binding-for-global-variable?noredirect=1&lq=1
## https://stackoverflow.com/questions/9439256/how-can-i-handle-r-cmd-check-no-visible-binding-for-global-variable-notes-when
profile <- NULL # does *not* affect the subset() call that follows
if (missing(j))
stop("must provide an integer vector to access, e.g. argo[[\"profile\", 1:3]]")
return(subset(x, profile %in% j))
}
if (i %in% c("CT", "N2", "SA", "sigmaTheta", "theta")) {
## FIXME: should we prefer e.g. salinityAdjusted or salinity?
names <- names(x@data)
salinity <- x@data[[if ("salinityAdjusted" %in% names) "salinityAdjusted" else "salinity"]]
pressure <- x@data[[if ("pressureAdjusted" %in% names) "pressureAdjusted" else "pressure"]]
temperature <- x@data[[if ("temperatureAdjusted" %in% names) "temperatureAdjusted" else "temperature"]]
dim <- dim(salinity)
## Do not need longitude and latitude if eos="unesco", but retain for code clarity
longitude <- rep(x@data$longitude, each=dim[1])
latitude <- rep(x@data$latitude, each=dim[1])
if (i == "CT") {
res <- gsw_CT_from_t(x[["SA"]], temperature, pressure)
} else if (i == "N2") {
##nprofile <- dim[2]
res <- array(NA_real_, dim=dim)
for (i in seq_len(dim[2])) {
##message("i=",i, ", nprofile=", nprofile)
##if (i == 14) browser()
if (sum(!is.na(pressure[,i])) > 2) {
ctd <- as.ctd(salinity=salinity[,i],
temperature=temperature[,i],
pressure=pressure[,i],
longitude=x@data$longitude[i],
latitude=x@data$latitude[i])
res[,i] <- swN2(ctd)
} else {
res[,i] <- rep(NA, length(salinity[,i]))
}
}
} else if (i == "SA") {
res <- gsw_SA_from_SP(salinity, pressure, longitude=longitude, latitude=latitude)
} else if (i == "sigmaTheta") {
res <- swSigmaTheta(salinity, temperature=temperature, pressure=pressure,
referencePressure=0, longitude=longitude, latitude=latitude,
eos=getOption("oceEOS", default="gsw"))
} else if (i == "theta") {
res <- swTheta(salinity, temperature=temperature, pressure=pressure,
referencePressure=0, longitude=longitude, latitude=latitude,
eos=getOption("oceEOS", default="gsw"))
} else {
stop("coding error: unknown item '", i, "'")
}
dim(res) <- dim
} else if (i == "depth") {
## This accessor added for issue 1333. Note that the
## fix for that issue was sometimes calling with
## vector-form argo object. I don't know how that vector
## form is arising, but it is likely an index without
## a drop=FALSE condition ... if I find it, I'll fix it,
## but the following works fine, so I don't really care too
## much.
if (is.matrix(x@data$pressure)) {
n <- dim(x@data$pressure)[1]
latitude <- matrix(rep(x@data$latitude, each=n),
nrow=n, byrow=TRUE)
res <- swDepth(x@data$pressure, latitude)
##. print("matrix ... lat and then pres... and the depth...")
##. print(latitude[1:3, 1:3])
##. print(x@data$pressure[1:3, 1:3])
##. print(res[1:3, 1:3])
} else {
res <- swDepth(x@data$pressure, x@data$latitude)
}
} else if (i == "latitude") {
res <- x@data$latitude
} else if (i == "longitude") {
res <- x@data$longitude
} else {
##message("FIXME: [[,argo-method calling next method")
res <- callNextMethod() # [[ defined in R/AllClass.R
}
res
})
#' Replace Parts of an Argo Object
#'
#' @param x an [argo-class] object.
#'
#' @template sub_subsetTemplate
#'
#' @family things related to argo data
setMethod(f="[[<-",
signature(x="argo", i="ANY", j="ANY"),
definition=function(x, i, j, ..., value) {
callNextMethod(x=x, i=i, j=j, ...=..., value=value) # [[<-
})
setMethod(f="initialize",
signature="argo",
definition=function(.Object, time, id, longitude, latitude, salinity, temperature, pressure, filename, dataMode, ...) {
.Object <- callNextMethod(.Object, ...)
if (!missing(time)) .Object@data$time <- time
if (!missing(id)) .Object@metadata$id <- id
if (!missing(longitude)) .Object@data$longitude <- longitude
if (!missing(latitude)) .Object@data$latitude <- latitude
if (!missing(salinity)) .Object@data$salinity <- salinity
if (!missing(temperature)) .Object@data$temperature <-temperature
if (!missing(pressure)) .Object@data$pressure <- pressure
.Object@metadata$filename <- if (missing(filename)) "" else filename
.Object@metadata$dataMode <- if (missing(dataMode)) "" else dataMode
.Object@processingLog$time <- presentTime()
.Object@processingLog$value <- "create 'argo' object"
return(.Object)
})
## a local function -- no need to pollute namespace with it
maybeLC <- function(s, lower)
if (lower) tolower(s) else s
## a local function -- no need to pollute namespace with it
getData <- function(file, name, quiet=FALSE)
{
res <- NA
## wrap in capture.output to silence what seems like direct printing to stdout()
## or stderr() by ncvar_get().
capture.output(
{
res <- try(ncdf4::ncvar_get(file, name), silent=TRUE)
})
if (inherits(res, "try-error")) {
if (!quiet)
warning(file$filename, " has no variable named '", name, "'\n", sep='')
res <- NULL
}
if (is.array(res) && 1 == length(dim(res))) res <- matrix(res) else res
}
#' Convert Argo Data Name to Oce Name
#'
#' This function is used internally by [read.argo()] to convert Argo-convention
#' data names to oce-convention names. Users should not call this directly, since
#' its return value may be changed at any moment (e.g. to include units as well
#' as names).
#'
#'
#' The inference of names was done
#' by inspection of some data files, based on reference 1. It should be noted,
#' however, that the data files examined contain some names that are not
#' undocumented in reference 1, and others that are listed only in its changelog,
#' with no actual definitions being given. For example, the files had six distinct
#' variable names that seem to relate to phase in the oxygen sensor, but
#' these are not translated by the present function because these
#' variable names are not defined in reference 1, or not defined uniquely
#' in reference 2.
#'
#' The names are converted with
#' [gsub()], using the `ignore.case` argument of the present
#' function.
#' The procedure
#' is to first handle the items listed in the following table, with string
#' searches anchored to the start of the string. After that,
#' the qualifiers
#' `_ADJUSTED`, `_ERROR` and `_QC`,
#' are translated to `Adjusted`, `Error`, and `QC`, respectively.
#'
#' \tabular{ll}{
#' **Argo name** \tab **oce name**\cr
#' `BBP` \tab `bbp`\cr
#' `BETA_BACKSCATTERING` \tab `betaBackscattering`\cr
#' `BPHASE_OXY` \tab `bphaseOxygen`\cr
#' `CDOM` \tab `CDOM`\cr
#' `CNDC` \tab `conductivity`\cr
#' `CHLA` \tab `chlorophyllA`\cr
#' `CP` \tab `beamAttenuation`\cr
#' `CYCLE_NUMBER` \tab `cycleNumber`\cr
#' `DATA_CENTRE` \tab `dataCentre`\cr
#' `DATA_MODE` \tab `dataMode`\cr
#' `DATA_STATE_INDICATOR` \tab `dataStateIndicator`\cr
#' `DC_REFERENCE` \tab `DCReference`\cr
#' `DIRECTION` \tab `direction`\cr
#' `DOWN_IRRADIANCE` \tab `downwellingIrradiance`\cr
#' `DOWNWELLING_PAR` \tab `downwellingPAR`\cr
#' `FIRMWARE_VERSION` \tab `firmwareVersion`\cr
#' `FIT_ERROR_NITRATE` \tab `fitErrorNitrate`\cr
#' `FLUORESCENCE_CDOM` \tab `fluorescenceCDOM`\cr
#' `FLUORESCENCE_CHLA` \tab `fluorescenceChlorophyllA`\cr
#' `INST_REFERENCE` \tab `instReference`\cr
#' `JULD` \tab `juld` (and used to compute `time`)\cr
#' `JULD_QC_LOCATION` \tab `juldQCLocation`\cr
#' `LATITUDE` \tab `latitude`\cr
#' `LONGITUDE` \tab `longitude`\cr
#' `MOLAR_DOXY` \tab `oxygenUncompensated`\cr
#' `PH_IN_SITU_FREE` \tab `pHFree`\cr
#' `PH_IN_SITU_TOTAL` \tab `pH`\cr
#' `PI_NAME` \tab `PIName`\cr
#' `PLATFORM_NUMBER` \tab `id`\cr
#' `POSITION_ACCURACY` \tab `positionAccuracy`\cr
#' `POSITIONING_SYSTEM` \tab `positioningSystem`\cr
#' `PROFILE` \tab `profile`\cr
#' `PROJECT_NAME` \tab `projectName`\cr
#' `RAW_DOWNWELLING_IRRADIANCE` \tab `rawDownwellingIrradiance`\cr
#' `RAW_DOWNWELLING_PAR` \tab `rawDownwellingPAR`\cr
#' `RAW_UPWELLING_RADIANCE` \tab `rawUpwellingRadiance`\cr
#' `STATION_PARAMETERS` \tab `stationParameters`\cr
#' `TEMP` \tab `temperature`\cr
#' `TEMP_CPU_CHLA` \tab `temperatureCPUChlorophyllA`\cr
#' `TEMP_DOXY` \tab `temperatureOxygen`\cr
#' `TEMP_NITRATE` \tab `temperatureNitrate`\cr
#' `TEMP_PH` \tab `temperaturePH`\cr
#' `TEMP_SPECTROPHOTOMETER_NITRATE` \tab `temperatureSpectrophotometerNitrate`\cr
#' `TILT` \tab `tilt`\cr
#' `TURBIDITY` \tab `turbidity`\cr
#' `UP_RADIANCE` \tab `upwellingRadiance`\cr
#' `UV_INTENSITY` \tab `UVIntensity`\cr
#' `UV_INTENSITY_DARK_NITRATE` \tab `UVIntensityDarkNitrate`\cr
#' `UV_INTENSITY_NITRATE` \tab `UVIntensityNitrate`\cr
#' `VRS_PH` \tab `VRSpH`\cr
#' `WMO_INST_TYPE` \tab `WMOInstType`\cr
#'}
#'
#' @param names vector of character strings containing names in the Argo convention.
#'
#' @param ignore.case a logical value passed to [gsub()], indicating whether to
#' ignore the case of input strings. The default is set to `TRUE` because some data
#' files use lower-case names, despite the fact that the Argo documentation specifies
#' upper-case.
#'
#' @return A character vector of the same length as `names`, but with
#' replacements having been made for all known quantities.
#'
#' @references
#' 1. Argo User's Manual Version 3.3, Nov 89th, 2019, available at
#' \url{https://archimer.ifremer.fr/doc/00187/29825/}
#'
#' 2. Argo list of parameters in an excel spreadsheet, available at
#' \url{http://www.argodatamgt.org/content/download/27444/187206/file/argo-parameters-list-core-and-b.xlsx}
#'
#' @family things related to argo data
argoNames2oceNames <- function(names, ignore.case=TRUE)
{
## do NOT change the order below, because we are working with partial strings.
names <- gsub("^BBP([0-9_]*)", "BBP\\1", names, ignore.case=ignore.case)
names <- gsub("^BETA_BACKSCATTERING([0-9_]*)", "betaBackscattering\\1", names, ignore.case=ignore.case)
names <- gsub("^BPHASE_DOXY", "bphaseOxygen", names, ignore.case=ignore.case)
names <- gsub("^CHLA", "chlorophyllA", names, ignore.case=ignore.case)
names <- gsub("^CDOM", "CDOM", names, ignore.case=ignore.case)
names <- gsub("^CNDC([0-9_]*)", "conductivity\\1", names, ignore.case=ignore.case)
names <- gsub("^CP([0-9_]*)", "beamAttenuation\\1", names, ignore.case=ignore.case)
names <- gsub("^CYCLE_NUMBER", "cycleNumber", names, ignore.case=ignore.case)
names <- gsub("^DOWN_IRRADIANCE", "downwellingIrradiance", names, ignore.case=ignore.case)
names <- gsub("^DOWNWELLING_PAR", "downwellingPAR", names, ignore.case=ignore.case)
names <- gsub("^FIT_ERROR_NITRATE", "fitErrorNitrate", names, ignore.case=ignore.case) # put before CHLA
names <- gsub("^FLUORESCENCE_CDOM", "fluorescenceCDOM", names, ignore.case=ignore.case) # put before CHLA
names <- gsub("^FLUORESCENCE_CHLA", "fluorescenceChlorophyllA", names, ignore.case=ignore.case) # put before CHLA
names <- gsub("^MOLAR_DOXY", "oxygenUncompensated", names, ignore.case=ignore.case)
names <- gsub("^PH_IN_SITU_FREE", "pHFree", names, ignore.case=ignore.case)
names <- gsub("^PH_IN_SITU_TOTAL", "pH", names, ignore.case=ignore.case)
names <- gsub("^TEMP_DOXY", "temperatureOxygen", names, ignore.case=ignore.case)
names <- gsub("^TEMP_NITRATE", "temperatureNitrate", names, ignore.case=ignore.case)
names <- gsub("^TEMP_PH", "temperaturePH", names, ignore.case=ignore.case)
names <- gsub("^TEMP_SPECTROPHOTOMETER_NITRATE", "temperatureSpectrophotometerNitrate", names, ignore.case=ignore.case)
names <- gsub("^TEMP_CPU_CHLA", "temperatureCPUChlA", names, ignore.case=ignore.case)
names <- gsub("^TEMP_VOLTAGE_DOXY", "temperatureVoltageOxygen", names, ignore.case=ignore.case)
names <- gsub("^TEMP_", "temperature_", names, ignore.case=ignore.case)
names <- gsub("^POSITION_ACCURACY", "positionAccuracy", names, ignore.case=ignore.case)
names <- gsub("^NITRATE", "nitrate", names, ignore.case=ignore.case)
names <- gsub("^DOXY", "oxygen", names, ignore.case=ignore.case)
names <- gsub("^PRES", "pressure", names, ignore.case=ignore.case)
names <- gsub("^PSAL", "salinity", names, ignore.case=ignore.case)
names <- gsub("^RAW_DOWNWELLING_IRRADIANCE", "rawDownwellingIrradiance", names, ignore.case=ignore.case)
names <- gsub("^RAW_DOWNWELLING_PAR", "rawDownwellingPAR", names, ignore.case=ignore.case)
names <- gsub("^RAW_UPWELLING_RADIANCE", "rawUpwellingRadiance", names, ignore.case=ignore.case)
names <- gsub("^TEMP([0-9_]*)$", "temperature\\1", names, ignore.case=ignore.case)
names <- gsub("^TILT([0-9_]*)$", "tilt\\1", names, ignore.case=ignore.case)
names <- gsub("^TRANSMITTANCE_PARTICLE_BEAM_ATTENUATION([0-9_]*)$",
"transmittanceParticleBeamAttenuation\\1", names, ignore.case=ignore.case)
names <- gsub("^TURBIDITY([0-9_]*)$", "turbidity\\1", names, ignore.case=ignore.case)
names <- gsub("^UP_RADIANCE", "upwellingRadiance", names, ignore.case=ignore.case)
names <- gsub("^UV_INTENSITY_DARK_NITRATE", "UVIntensityDarkNitrate", names, ignore.case=ignore.case)
names <- gsub("^UV_INTENSITY_NITRATE", "UVIntensityNitrate", names, ignore.case=ignore.case)
names <- gsub("^UV_INTENSITY", "UVIntensity", names, ignore.case=ignore.case)
names <- gsub("^VRS_PH", "VRSpH", names, ignore.case=ignore.case)
names <- gsub("_ADJUSTED", "Adjusted", names, ignore.case=ignore.case)
names <- gsub("_QC", "QC", names, ignore.case=ignore.case)
names <- gsub("_ERROR", "Error", names, ignore.case=ignore.case)
names
}
#' Subset an Argo Object
#'
#' Subset an argo object, either by selecting just the "adjusted" data
#' or by subsetting by pressure or other variables.
#'
#' @details
#' If `subset` is the string `"adjusted"`, then `subset`
#' replaces the station variables with their adjusted counterparts. In
#' the argo notation, e.g. `PSAL` is replaced with `PSAL_ADJUSTED`;
#' in the present notation, this means that `salinity` in the `data`
#' slot is replaced with `salinityAdjusted`, and the latter is deleted.
#' Similar replacements are also done with the flags stored in the `metadata`
#' slot.
#'
#' If `subset` is an expression, then the action is somewhat similar
#' to other `subset` functions, but with the restriction that
#' only one independent variable may be
#' used in in any call to the function, so that
#' repeated calls will be necessary to subset based on more than one
#' independent variable. Subsetting may be done by anything
#' stored in the data, e.g. `time`,
#' `latitude`, `longitude`, `profile`, `dataMode`,
#' or `pressure` or by `profile` (a made-up variable)
#' or `id` (from the `metadata` slot). Note that subsetting by `pressure`
#' preserves matrix shape, by setting discarded values to `NA`, as opposed
#' to dropping data (as is the case with `time`, for example).
#'
#' @param x an [argo-class] object.
#'
#' @param subset An expression indicating how to subset `x`.
#'
#' @param ... optional arguments, of which only the first is examined. The
#' only possibility is `within`, a polygon enclosing data to be
#' retained. This must be either a list or data frame, containing items
#' named either `x` and `y` or `longitude` and
#' `latitude`; see Example 4. If `within` is given,
#' then `subset` is ignored.
#'
#' @return An argo object.
#'
#' @examples
#' library(oce)
#' data(argo)
#'
#' # Example 1: buset by time, longitude, and pressure
#' par(mfrow=c(2,2))
#' plot(argo)
#' plot(subset(argo, time > mean(time)))
#' plot(subset(argo, longitude > mean(longitude)))
#' plot(subset(argoGrid(argo), pressure > 500 & pressure < 1000), which=5)
#'
#' # Example 2: restrict attention to delayed-mode profiles.
#' par(mfrow=c(1, 1))
#' plot(subset(argo, dataMode == "D"))
#'
#' # Example 3: contrast corrected and uncorrected data
#' par(mfrow=c(1,2))
#' plotTS(argo)
#' plotTS(subset(argo, "adjusted"))
#'
#' # Example 4. Subset by a polygon determined with locator()
#'\dontrun{
#' par(mfrow=c(2, 1))
#' plot(argo, which="map")
#' bdy <- locator(4) # Click the mouse on 4 boundary points
#' argoSubset <- subset(argo, within=bdy)
#' plot(argoSubset, which="map")
#'}
#'
#' @author Dan Kelley
#'
#' @family things related to argo data
#' @family functions that subset oce objects
#' @aliases subset.argo
setMethod(f="subset",
signature="argo",
definition=function(x, subset, ...) {
subsetString <- paste(deparse(substitute(subset)), collapse=" ")
res <- x
dots <- list(...)
dotsNames <- names(dots)
withinGiven <- length(dots) && ("within" %in% dotsNames)
debug <- getOption("oceDebug")
if (length(dots) && ("debug" %in% names(dots)))
debug <- dots$debug
if (withinGiven) {
oceDebug(debug, "subsetting with 'within' method")
polygon <- dots$within
if (!is.data.frame(polygon) && !is.list(polygon))
stop("'within' must be a data frame or a polygon")
polygonNames <- names(polygon)
lonp <- if ("x" %in% polygonNames) {
polygon$x
} else if ("longitude" %in% polygonNames) {
polygon$longitude
} else {
stop("'within' must contain either 'x' or 'longitude'")
}
latp <- if ("y" %in% polygonNames) {
polygon$y
} else if ("latitude" %in% polygonNames) {
polygon$latitude
} else {
stop("'within' must contain either 'y' or 'latitude'")
}
lon <- x[["longitude", "byStation"]]
lat <- x[["latitude", "byStation"]]
if (requireNamespace("sp", quietly=TRUE)) {
keep <- 1==sp::point.in.polygon(lon, lat, lonp, latp)
} else {
stop("cannot use 'within' because the 'sp' package is not installed")
}
## Metadata
for (name in names(x@metadata)) {
oceDebug(debug, "subsetting metadata item named '", name, "'\n", sep="")
## Pass oce-generated things through directly.
if (name %in% c("units", "flags", "filename", "flagScheme", "dataNamesOriginal")) {
##.message(" ... special case, so passed directly")
next
}
item <- x@metadata[[name]]
## Handle things that are encoded as characters in a string,
## namely 'direction', 'juldQc', and 'positionQc'.
if (is.character(item) && length(item) == 1) {
##.message("'", name, "' is character-encoded")
res@metadata[[name]] <- paste(strsplit(item,"")[[1]][keep],collapse="")
##.message(" ... ok")
} else if (is.vector(name)) {
##.message("'", name, "' is a vector")
res@metadata[[name]] <- item[keep]
##.message(" ... ok")
} else if (is.matrix(name)) {
##.message("'", name, "' is a matrix")
res@metadata[[name]] <- item[, keep]
##.message(" ... ok")
} else {
stop("cannot subset metadata item named '", name, "' because it is not a length-one string, a vector, or a matrix")
}
}
## Data
for (name in names(x@data)) {
oceDebug(debug, "subsetting data item named '", name, "'\n", sep="")
item <- x@data[[name]]
if ("time" == name) {
##.message("'", name, "' is time (not a vector)")
res@data$time <- item[keep]
##.message(" ... ok")
} else if (is.vector(item)) {
##.message("'", name, "' is vector")
res@data[[name]] <- item[keep]
##.message(" ... ok")
} else if (is.matrix(item)) {
##.message("'", name, "' is matrix")
res@data[[name]] <- item[, keep]
##.message(" ... ok")
} else {
stop("argo object has data item '", name, "' that is neither a vector nor a matrix, so we cannot subset it")
}
}
res@processingLog <- processingLogAppend(res@processingLog,
paste("subset(x, within) kept ", sum(keep), " of ",
length(keep), " stations", sep=""))
} else {
if (is.character(substitute(subset))) {
if (subset != "adjusted")
stop("if subset is a string, it must be \"adjusted\"")
dataNames <- names(x@data)
## Seek 'Adjusted' names
adjustedIndices <- grep(".*Adjusted$", dataNames)
for (i in adjustedIndices) {
adjusted <- dataNames[i]
base <- gsub("Adjusted$", "", adjusted)
adjustedError <- paste(adjusted, "Error", sep="")
##> message(" base: ", base)
##> message(" adjusted: ", adjusted)
##> message(" adjustedError: ", adjustedError)
res@data[[base]] <- res@data[[adjusted]]
res@data[[adjusted]] <- NULL
res@data[[adjustedError]] <- NULL
}
flagNames <- names(x@metadata$flags)
adjustedIndices <- grep("Adjusted", flagNames)
##> message("FLAGS")
##> message("flagNames...");print(flagNames)
##> message("adjustedIndices");print(adjustedIndices)
for (i in adjustedIndices) {
adjusted <- flagNames[i]
base <- gsub("Adjusted", "", adjusted)
adjustedError <- paste(adjusted, "Error", sep="")
##> message(" base: ", base)
##> message(" adjusted: ", adjusted)
##> message(" adjustedError: ", adjustedError)
res@metadata$flags[[base]] <- res@metadata$flags[[adjusted]]
res@metadata$flags[[adjusted]] <- NULL
res@metadata$flags[[adjustedError]] <- NULL
}
res@processingLog <- processingLogAppend(res@processingLog,
paste("subset.argo(x, subset=\"",
as.character(subset), "\")", sep=""))
} else {
subsetString <- paste(deparse(substitute(subset)), collapse=" ")
res <- x
if (length(grep("time", subsetString)) ||
length(grep("longitude", subsetString)) || length(grep("latitude", subsetString))) {
keep <- eval(substitute(subset), x@data, parent.frame(2))
} else if (length(grep("id", subsetString))) {
## add id into the data, then do as usual
tmp <- x@data
tmp$id <- x@metadata$id
keep <- eval(substitute(subset), tmp, parent.frame(2))
rm(tmp)
} else if (length(grep("profile", subsetString))) {
## add profile into the data, then do as usual
tmp <- x@data
tmp$profile <- seq_along(x@data$time)
keep <- eval(substitute(subset), tmp, parent.frame(2))
rm(tmp)
} else if (length(grep("pressure", subsetString))) {
## issue1628 ## check that it is a "gridded" argo
## issue1628 gridded <- ifelse(all(apply(x@data$pressure, 1, diff) == 0, na.rm=TRUE), TRUE, FALSE)
## issue1628 if (gridded) {
## issue1628 x@data$pressure <- x@data$pressure[, 1] ## FIXME: have to convert pressure to vector
## issue1628 keep <- eval(substitute(subset), x@data, parent.frame(2))
## issue1628 x@data$pressure <- res@data$pressure ## FIXME: convert back to original for subsetting below
## issue1628 } else {
## issue1628 stop("cannot subset ungridded argo by pressure -- use argoGrid() first", call.=FALSE)
## issue1628 }
keep <- eval(substitute(subset), x@data, parent.frame(2))
} else if (length(grep("dataMode", subsetString))) {
keep <- eval(substitute(subset), x@metadata, parent.frame(2))
} else {
stop("can only subset by time, longitude, latitude, pressure, dataMode, and not by combinations", call.=FALSE)
}
if (length(grep("pressure", subsetString))) {
## Now do the subset. Note that we preserve matrix dimensions, by setting
## discarded values to NA.
fieldname <- names(x@data)
for (field in fieldname) {
if (field != 'time' & field != 'longitude' & field != 'latitude') { # DEBUG: see issue 1327
ifield <- which(field == fieldname)
##debug message("ifield=", ifield, ", field=", field,
##debug "\n\tlength(keep)=", length(keep),
##debug "\n\tsum(keep)=", sum(keep))
if (is.matrix(res@data[[ifield]])) {
res@data[[ifield]][!keep] <- NA
} else {
res@data[[ifield]][!keep] <- NA
}
}
}
fieldname <- names(x@metadata$flags)
for (field in fieldname) {
ifield <- which(field == fieldname)
res@metadata$flags[[ifield]][!keep] <- NA
}
## res@data$salinity <- x@data$salinity[keep, ]
## res@data$temperature <- x@data$temperature[keep, ]
## res@data$pressure <- x@data$pressure[keep, ]
res@processingLog <- processingLogAppend(res@processingLog, paste("subset.argo(x, subset=", subsetString, ")", sep=""))
} else {
res@data$time <- x@data$time[keep]
res@data$longitude <- x@data$longitude[keep]
res@data$latitude <- x@data$latitude[keep]
res@data$profile <- x@data$profile[keep]
res@metadata$dataMode <- x@metadata$dataMode[keep]
fieldname <- names(x@data)
for (field in fieldname) {
if (field != 'time' && field != 'longitude' && field != 'latitude' && field != 'profile') {
ifield <- which(field == fieldname)
res@data[[ifield]] <- if (is.matrix(x@data[[ifield]]))
x@data[[ifield]][, keep, drop=FALSE] else x@data[[ifield]][keep]
}
}
fieldname <- names(x@metadata$flags)
for (field in fieldname) {
ifield <- which(field == fieldname)
res@metadata$flags[[ifield]] <- res@metadata$flags[[ifield]][, keep]
}
#if (sum(keep) < 1) warning("In subset.argo() :\n removed all profiles", call.=FALSE)
## res@data$salinity <- x@data$salinity[, keep]
## res@data$temperature <- x@data$temperature[, keep]
## res@data$pressure <- x@data$pressure[, keep]
}
res@processingLog <- processingLogAppend(res@processingLog, paste("subset.argo(x, subset=", subsetString, ")", sep=""))
}
}
res
})
#' Summarize an Argo Object
#'
#' @description Summarizes some of the data in an `argo` object.
#'
#' @details Pertinent summary information is presented.
#' @param object}{an object of class `"argo"`, usually, a result of a
#' call to [read.argo()].
#' @param ... Further arguments passed to or from other methods.
#'
#' @return A matrix containing statistics of the elements of the `data` slot.
#' @examples
#' library(oce)
#' data(argo)
#' summary(argo)
#'
#' @author Dan Kelley
#' @family things related to argo data
#' @aliases summary.argo
setMethod(f="summary",
signature="argo",
definition=function(object, ...) {
cat("Argo Summary\n------------\n\n")
showMetadataItem(object, "filename", "Source: ", quote=TRUE)
nid <- length(unique(object@metadata$id))
if (1 == nid)
cat("* id: \"", object@metadata$id[1], "\"\n", sep="")
else cat("* id list: \"", object@metadata$id[1], "\", \"", object@metadata$id[2], "\", ...\n", sep="")
nD <- sum(object@metadata$dataMode == "D")
nA <- sum(object@metadata$dataMode == "A")
nR <- sum(object@metadata$dataMode == "R")
cat("* Profiles: ", nD, " delayed; ", nA, " adjusted; ", nR, " realtime", "\n", sep="")
invisible(callNextMethod()) # summary
})
ncdfFixMatrix <- function(x)
{
if (length(dim(x)) == 1)
x <- as.vector(x)
x
}
#' Grid Argo float data
#'
#' Grid an Argo float, by interpolating to fixed pressure levels.
#' The gridding is done with [approx()]. If there is
#' sufficient user demand, other methods may be added, by analogy to
#' [sectionGrid()].
#'
#' @template flagDeletionTemplate
#'
#' @param argo A `argo` object to be gridded.
#'
#' @param p Optional indication of the pressure levels to which interpolation
#' should be done. If this is not supplied, the pressure levels will be
#' calculated based on the existing values, using medians. If `p="levitus"`,
#' then pressures will be set to be those of the Levitus atlas, given by
#' [standardDepths()], trimmed to the maximum pressure in `argo`.
#' If `p` is a single numerical value, it is taken as the number of
#' subdivisions to use in a call to [seq()] that has range from 0 to the
#' maximum pressure in `argo`. Finally, if a vector numerical values is
#' provided, then it is used as is.
#'
#' @param debug A flag that turns on debugging. Higher values provide deeper
#' debugging.
#'
#' @param ... Optional arguments to [approx()], which is used to do the
#' gridding.
#'
#' @return x an [argo-class] object.
#'
#' @examples
#' library(oce)
#' data(argo)
#' g <- argoGrid(argo, p=seq(0, 100, 1))
#' par(mfrow=c(2,1))
#' t <- g[["time"]]
#' z <- -g[["pressure"]][,1]
#' ## Set zlim because of spurious temperatures.
#' imagep(t, z, t(g[['temperature']]), ylim=c(-100,0), zlim=c(0,20))
#' imagep(t, z, t(g[['salinity']]), ylim=c(-100,0))
#'
#' @family things related to argo data
#' @author Dan Kelley and Clark Richards
argoGrid <- function(argo, p, debug=getOption("oceDebug"), ...)
{
oceDebug(debug, "argoGrid() {\n", sep="", unindent=1)
warningMessages <- NULL
dim <- dim(argo@data$pressure)
## ndepth <- dim[1]
nprofile <- dim[2]
## FIXME: modify sal, temp, and pre. In the end, pre constant along first index
res <- argo
res[["flags"]] <- NULL
warningMessages <- c(warningMessages,
"Data flags are omitted from the gridded argo object. Use handleFlags() first to remove bad data.")
pressure <- argo[["pressure"]]
if (missing(p)) {
pt <- apply(pressure, 1, median, na.rm=TRUE)
} else if (length(p) == 1 && p == "levitus") {
pt <- standardDepths()
pt <- pt[pt < max(pressure, na.rm=TRUE)]
} else if (is.numeric(p)) {
if (length(p) == 1) {
if (p < 1)
stop("'p' must exceed 1")
pt <- seq(0, max(pressure, na.rm=TRUE), length.out=p)
} else {
pt <- p
}
} else {
stop("value of 'p' must be numeric, or \"levitus\"")
}
##message("pt=c(", paste(round(pt), collapse=","), ")")
npt <- length(pt)
res@data$pressure <- matrix(NA, ncol=nprofile, nrow=npt)
for (field in names(res@data)) {
if (!(field %in% c('time', 'longitude', 'latitude'))) {
res@data[[field]] <- matrix(NA, ncol=nprofile, nrow=npt)
for (profile in 1:nprofile) {
ndata <- sum(!is.na(argo@data[[field]][, profile]))
if (ndata > 2 && sum(is.finite(diff(pressure[, profile])))
&& 0 < max(abs(diff(pressure[, profile])), na.rm=TRUE)) {
res@data[[field]][, profile] <- approx(pressure[, profile], argo@data[[field]][, profile], pt, ...)$y
} else {
res@data[[field]][, profile] <- rep(NA, npt)
}
res@data$pressure[, profile] <- pt
}
}
}
res@processingLog <- processingLogAppend(res@processingLog, paste("Grid to regular pressures with: ", deparse(match.call()), sep="", collapse=""))
for (w in warningMessages)
res@processingLog <- processingLogAppend(res@processingLog, w)
res
}
argoDecodeFlags <- function(f) # local function
{
res <- unlist(lapply(seq_along(f), function(i) strsplit(f[i], split="")))
dim(res) <- c(length(res)/length(f), length(f))
mode(res) <- "numeric"
res
}
#' Read an Argo Data File
#'
#' `read.argo` is used to read an Argo file, producing an object of type
#' `argo`. The file must be in the ARGO-style NetCDF format described at
#' in the Argo documentation (see references 2 and 3).
#'
#' @details
#'
#' Items are inferred from the data file in a straightforward way, using
#' [ncdf4::ncvar_get()], converting from one-column matrices
#' to vectors, and trimming leading and trailing blank space in character
#' values, using [trimString()].
#'
#' Items are renamed from the argo ('snake case') forms to oce ('camel
#' case') forms with [argoNames2oceNames()]. For example,
#' `FIRMWARE_VERSION` in the data file is renamed as
#' `firmwareVersion` in the return value.
#' Note that some files use upper-case for items, while other files
#' use lower-case for the same things; `read.argo` attempts
#' to ignore this variation.
#'
#' See the Argo documentation (see references 2 and 3) for some details on what files contain.
#' Many items listed in section 2.2.3 of reference 3 are read from the
#' file and stored in the `metadata` slot, with the exception of
#' `longitude` and `latitude`, which are stored in the
#' `data` slot, alongside hydrographic information.
#'
#' The following global attributes stored within the netcdf file are stored in the
#' `metadata` slot: `title`, `institution`, `source`,
#' `history`, `references`, `userManualVersion`, `conventions`,
#' and `featureType`. These names are identical to those in the netcdf
#' file, except that `userManualVersion` is named
#' `user_manual_version` in the file, and `conventions` is
#' named `Conventions` in the file.
#'
#' It is assumed that the profile data are as listed in the NetCDF variable
#' called `STATION_PARAMETERS`. Each item can have variants, as
#' described in Sections 2.3.4 of reference 3.
#' For example, if `"PRES"` is found in `STATION_PARAMETERS`,
#' then `PRES` (pressure) data are sought in the file, along with
#' `PRES_QC`, `PRES_ADJUSTED`, `PRES_ADJUSTED_QC`, and
#' `PRES_ERROR`. The same pattern works for other profile data. The variables
#' are stored with different names within the resultant [argo-class]
#' object, to match with `oce` conventions. Thus, `PRES` gets renamed
#' `pressure`, while `PRES_ADJUSTED` gets renamed `pressureAdjusted`,
#' and `PRES_ERROR` gets renamed `pressureError`; all of these are
#' stored in the `data` slot. Meanwhile, the quality-control flags
#' `PRES_QC` and `PRES_ADJUSTED_QC` are stored as `pressure`
#' and `pressureAdjusted` in the `metadata$flags` slot.
#'
#' Once a predefined series of items are inferred and stored in either the
#' `metadata` or `data` slot, `read.argo` then reads all
#' non-empty variables in the file, storing them in the `metadata`
#' slot, using with the original ('snake case') name that is used in
#' the data file. (Note that the sample dataset accessed with `data(argo)`
#' lacks metadata items with names starting with `HISTORY_`, because
#' those are zero-length in the source file.)
#'
#' @param file a character string giving the name of the file to load.
#'
#' @param debug a flag that turns on debugging. Set to 1 to get a moderate amount
#' of debugging information, or to 2 to get more.
#'
#' @param processingLog if provided, the action item to be stored in the log.
#' (Typically only provided for internal calls; the default that it provides is
#' better for normal calls by a user.)
#'
#' @param ... additional arguments, passed to called routines.
#'
#' @return
#' An [argo-class] object.
#'
#' @examples
#'\dontrun{
#' ## Example 1: read from a local file
#' library(oce)
#' d <- read.argo("/data/OAR/6900388_prof.nc")
#' summary(d)
#' plot(d)
#'
#' ## Example 2: construct URL for download (brittle)
#' id <- "6900388"
#' url <- "https://www.usgodae.org/ftp/outgoing/argo"
#' if (!length(list.files(pattern="argo_index.txt")))
#' download.file(paste(url, "ar_index_global_meta.txt", sep="/"), "argo_index.txt")
#' index <- readLines("argo_index.txt")
#' line <- grep(id, index)
#' if (0 == length(line)) stop("id ", id, " not found")
#' if (1 < length(line)) stop("id ", id, " found multiple times")
#' dac <- strsplit(index[line], "/")[[1]][1]
#' profile <- paste(id, "_prof.nc", sep="")
#' float <- paste(url, "dac", dac, id, profile, sep="/")
#' download.file(float, profile)
#' argo <- read.argo(profile)
#' summary(argo)
#'}
#'
#'
#' @seealso
#' The documentation for the [argo-class] class explains the structure of argo
#' objects, and also outlines the other functions dealing with them.
#'
#' @references
#' 1. \url{http://www.argo.ucsd.edu/}
#'
#' 2. Argo User's Manual Version 3.3, Nov 89th, 2019, available at
#' \url{https://archimer.ifremer.fr/doc/00187/29825/}
#'
#' 3. User's Manual (ar-um-02-01) 13 July 2010, available at
#' \url{http://www.argodatamgt.org/content/download/4729/34634/file/argo-dm-user-manual-version-2.3.pdf},
#' which is the main document describing argo data.
#'
#' @section Data sources:
#' Argo data are made available at several websites. A bit of detective
#' work can be required to track down the data.
#'
#' Some servers provide data for floats that surfaced in a given ocean
#' on a given day, the anonymous FTP server
#' \code{ftp://usgodae.org/pub/outgoing/argo/geo/} being an example.
#'
#' Other servers provide data on a per-float basis. A complicating
#' factor is that these data tend to be categorized by "dac" (data
#' archiving centre), which makes it difficult to find a particular
#' float. For example,
#' \code{https://www.usgodae.org/ftp/outgoing/argo/} is the top level of
#' a such a repository. If the ID of a float is known but not the
#' "dac", then a first step is to download the text file
#' \code{http://www.usgodae.org/ftp/outgoing/argo/ar_index_global_meta.txt}
#' and search for the ID. The first few lines of that file are header,
#' and after that the format is simple, with columns separated by slash
#' (`/`). The dac is in the first such column and the float ID in the
#' second. A simple search will reveal the dac.
#' For example `data(argo)` is based on float 6900388, and the line