forked from dankelley/oce
-
Notifications
You must be signed in to change notification settings - Fork 0
/
ctd.R
5733 lines (5677 loc) · 271 KB
/
ctd.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 CTD (or general hydrographic) Data
#'
#' This class stores hydrographic data such as measured with a CTD (conductivity,
#' temperature, depth) instrument, or with other systems that produce
#' similar data. Data repositories may store conductivity, temperature
#' and depth, as in the instrument name, but it is also common to store
#' salinity, temperature and pressure instead (or in addition). For this
#' reason, `ctd` objects are required to hold `salinity`,
#' `temperature` and `pressure` in their `data` slot,
#' with other data being optional. Formulae are available for converting
#' between variants of these data triplets, e.g. [swSCTp()]
#' can calculate `salinity` given `conductivity`, `temperature`
#' and `pressure`, and these are used by the main functions that
#' create `ctd` objects. For example, if [read.ctd.sbe()]
#' is used to read a Seabird file that contains only conductivity, temperature
#' and pressure, then that function will automatically append a data
#' item to hold salinity. Since [as.ctd()] does the same with
#' salinity, the result this is that all `ctd` objects hold `salinity`,
#' `temperature` and `pressure`, which are henceforth called
#' the three basic quantities.
#'
#' Different units and scales are permitted for the three basic quantities, and
#' most `oce` functions check those units and scales before
#' doing calculations (e.g. of seawater density), because those calculations
#' demand certain units and scales. The way this is handled is that the
#' accessor function \code{\link{[[,ctd-method}}] returns values in standardized
#' form. For example, a `ctd` object might hold temperature defined on the
#' IPTS-68 scale, but e.g. `ctd[["temperature"]]` returns a value on the ITS-90
#' scale. (The conversion is done with [T90fromT68()].) Similarly,
#' pressure may be stored in either dbars or PSI, but e.g. `ctd[["pressure"]]`
#' returns a value in dbars, after dividing by 0.689476 if the value is
#' stored in PSI. Luckily, there is (as of early 2016) only one salinity scale in
#' common use in data files, namely PSS-78.
#'
#' @templateVar class ctd
#'
#' @templateVar dataExample The key items stored in this slot are: `salinity`, `temperature`, and `pressure`, although in many instances there are quite a few additional items.
#'
#' @templateVar metadataExample An example of the former might be the location at which a `ctd` measurement was made, stored in `longitude` and `latitude`, and of the latter might be `filename`, the name of the data source.
#'
#' @template slot_summary
#'
#' @template slot_put
#'
#' @template slot_get
#'
#' @section Reading/creating `ctd` objects:
#' A file containing CTD profile data may be read with
#' [read.ctd()], and a CTD object can also be created with
#' [as.ctd()]. See [read.ctd()] for references on data
#' formats used in CTD files. Data can also be assembled into
#' `ctd` objects with [as.ctd()].
#'
#' Statistical summaries are provided by [summary,ctd-method()], while
#' [show()] displays an overview.
#'
#' CTD objects may be plotted with [plot,ctd-method()], which does much of its
#' work by calling [plotProfile()] or [plotTS()], both of
#' which can also be called by the user, to get fine control over the plots.
#'
#' A CTD profile can be isolated from a larger record with [ctdTrim()],
#' a task made easier when [plotScan()] is used to examine the results.
#' Towyow data can be split up into sets of profiles (ascending or descending)
#' with [ctdFindProfiles()]. CTD data may be smoothed and/or cast onto
#' specified pressure levels with [ctdDecimate()].
#'
#' As with all oce objects, low-level manipulation may be done with
#' [oceSetData()] and [oceSetMetadata()]. Additionally,
#' many of the contents of CTD objects may be altered with the
#' \code{\link{[[<-,ctd-method}} scheme,
#' and sufficiently skilled users may even manipulate the contents directly.
#'
#' @section Data sources:
#'
#' Archived CTD (and other) data may be found on servers such as
#' 1. \url{https://cchdo.ucsd.edu/}
#' 2. \url{https://www.nodc.noaa.gov/woce/wdiu/}
## 3. \url{https://www.nodc.noaa.gov/cgi-bin/OC5/SELECT/builder.pl}
#'
#' @examples
#'
#' # 1. Create a ctd object with fake data.
#' a <- as.ctd(salinity=35+1:3/10, temperature=10-1:3/10, pressure=1:3)
#' summary(a)
#'
#' # 2. Fix a typo in a station latitude (fake! it's actually okay)
#' data(ctd)
#' ctd <- oceSetMetadata(ctd, "latitude", ctd[["latitude"]]-0.001,
#' "fix latitude typo in log book")
#'
#' @author Dan Kelley
#'
#' @family things related to ctd data
#' @family classes provided by oce
setClass("ctd", contains="oce")
setAs("list", "ctd", function(from) {
as.ctd(from) #salinity=from$salinity, temperature=from$temperature, pressure=from$pressure)
})
#' A CTD profile in Halifax Harbour
#'
#' This is a CTD profile measured in Halifax Harbour in 2003, based
#' on [ctdRaw()], but trimmed to just the downcast with
#' [ctdTrim()], using indices inferred by inspection of the
#' results from [plotScan()].
#'
#' This station was sampled by students enrolled in the Dan Kelley's
#' Physical Oceanography class at Dalhousie University.
#' The data were acquired near the centre of the Bedford Basin of the
#' Halifax Harbour, during an October 2003 field trip of Dalhousie University's
#' Oceanography 4120/5120 class. The original `.cnv` data file had
#' temperature in the IPTS-68 scale, but this was converted to the more modern
#' scale using [T90fromT68()].
#'
#' @name ctd
#' @docType data
#'
#' @usage data(ctd)
#'
#' @examples
#'\dontrun{
#' library(oce)
#' data(ctd)
#' plot(ctd)
#'}
#'
#' @seealso The full profile (not trimmed to the downcast) is available as
#' `data(ctdRaw)`.
#'
#' @family datasets provided with oce
#' @family things related to ctd data
NULL
#' Seawater CTD Profile, Without Trimming of Extraneous Data
#'
#' This is sample CTD profile provided for testing. It includes not just the
#' (useful) portion of the dataset during which the instrument was being lowered,
#' but also data from the upcast and from time spent near the surface. Spikes are
#' also clearly evident in the pressure record. With such real-world wrinkles,
#' this dataset provides a good example of data that need trimming with
#' [ctdTrim()].
#'
#' This station was sampled by students enrolled in the Dan Kelley's
#' Physical Oceanography class at Dalhousie University.
#' The data were acquired near the centre of the Bedford Basin of the
#' Halifax Harbour, during an October 2003 field trip of Dalhousie University's
#' Oceanography 4120/5120 class. The original `.cnv` data file had
#' temperature in the IPTS-68 scale, but this was converted to the more modern
#' scale using [T90fromT68()].
#'
#' @name ctdRaw
#' @docType data
#'
#' @usage data(ctdRaw)
#'
#' @seealso A similar dataset (trimmed to the downcast) is available as
#' `data(ctd)`.
#'
#' @family things related to ctd data
#' @family datasets provided with oce
NULL
## DEVELOPERS: please pattern functions and documentation on this, for uniformity.
## DEVELOPERS: You will need to change the docs, and the 3 spots in the code
## DEVELOPERS: marked '# DEVELOPER 1:', etc.
#' @title Handle Flags in CTD Objects
#'
#' @param object a [ctd-class] object.
#'
#' @template handleFlagsTemplate
#'
#' @references
#' 1. \url{https://www.nodc.noaa.gov/woce/woce_v3/wocedata_1/whp/exchange/exchange_format_desc.htm}
#'
#' @examples
#' library(oce)
#' data(section)
#' stn <- section[["station", 100]]
#' # 1. Default: anything not flagged as 2 is set to NA, to focus
#' # solely on 'good', in the World Hydrographic Program scheme.
#' STN1 <- handleFlags(stn, flags=list(c(1, 3:9)))
#' data.frame(old=stn[["salinity"]], new=STN1[["salinity"]], salinityFlag=stn[["salinityFlag"]])
#'
#' # 2. Use bottle salinity, if it is good and ctd is bad
#' replace <- 2 == stn[["salinityBottleFlag"]] && 2 != stn[["salinityFlag"]]
#' S <- ifelse(replace, stn[["salinityBottle"]], stn[["salinity"]])
#' STN2 <- oceSetData(stn, "salinity", S)
#'
#' # 3. Use smoothed TS relationship to nudge questionable data.
#' f <- function(x) {
#' S <- x[["salinity"]]
#' T <- x[["temperature"]]
#' df <- 0.5 * length(S) # smooths a bit
#' sp <- smooth.spline(T, S, df=df)
#' 0.5 * (S + predict(sp, T)$y)
#' }
#' par(mfrow=c(1,2))
#' STN3 <- handleFlags(stn, flags=list(salinity=c(1,3:9)), action=list(salinity=f))
#' plotProfile(stn, "salinity", mar=c(3, 3, 3, 1))
#' p <- stn[["pressure"]]
#' par(mar=c(3, 3, 3, 1))
#' plot(STN3[["salinity"]] - stn[["salinity"]], p, ylim=rev(range(p)))
#'
#' # 4. Single-variable flags (vector specification)
#' data(section)
#' # Multiple-flag scheme: one per data item
#' A <- section[["station", 100]]
#' deep <- A[["pressure"]] > 1500
#' flag <- ifelse(deep, 7, 2)
#' for (flagName in names(A[["flags"]]))
#' A[[paste(flagName, "Flag", sep="")]] <- flag
#' Af <- handleFlags(A)
#' expect_equal(is.na(Af[["salinity"]]), deep)
#'
#' # 5. Single-variable flags (list specification)
#' B <- section[["station", 100]]
#' B[["flags"]] <- list(flag)
#' Bf <- handleFlags(B)
#' expect_equal(is.na(Bf[["salinity"]]), deep)
#'
#' @family things related to ctd data
setMethod("handleFlags", signature=c(object="ctd", flags="ANY", actions="ANY", where="ANY", debug="ANY"),
definition=function(object, flags=NULL, actions=NULL, where=NULL, debug=getOption("oceDebug")) {
## DEVELOPER 1: alter the next comment to explain your setup
if (is.null(flags)) {
flags <- defaultFlags(object)
if (is.null(flags))
stop("must supply 'flags', or use initializeFlagScheme() on the ctd object first")
}
if (is.null(actions)) {
actions <- list("NA") # DEVELOPER 3: alter this line to suit a new data class
names(actions) <- names(flags)
}
if (any(names(actions)!=names(flags)))
stop("names of flags and actions must match")
handleFlagsInternal(object=object, flags=flags, actions=actions, where=where, debug=debug)
})
#' @templateVar class ctd
#'
#' @templateVar note Since all the entries in the `data` slot of ctd objects are vectors, `i` must be a vector (either logical as in Example 1 or integer as in Example 2), or a function taking a `ctd` object and returning such a vector (see \dQuote{Indexing rules}).
#'
#' @template setFlagsTemplate
#'
#' @examples
#' library(oce)
#' # Example 1: Range-check salinity
#' data(ctdRaw)
#' ## Salinity and temperature range checks
#' qc <- ctdRaw
#' # Initialize flags to 2, meaning good data in the default
#' # scheme for handleFlags(ctd).
#' qc <- initializeFlags(qc, "salinity", 2)
#' qc <- initializeFlags(qc, "temperature", 2)
#' # Flag bad salinities as 4
#' oddS <- with(qc[["data"]], salinity < 25 | 40 < salinity)
#' qc <- setFlags(qc, name="salinity", i=oddS, value=4)
#' # Flag bad temperatures as 4
#' oddT <- with(qc[["data"]], temperature < -2 | 40 < temperature)
#' qc <- setFlags(qc, name="temperature", i=oddT, value=4)
#' # Compare results in TS space
#' par(mfrow=c(2, 1))
#' plotTS(ctdRaw)
#' plotTS(handleFlags(qc, flags=list(1, 3:9)))
#'
#' # Example 2: Interactive flag assignment based on TS plot, using
#' # WHP scheme to define 'acceptable' and 'bad' codes
#'\dontrun{
#' options(eos="gsw")
#' data(ctd)
#' qc <- ctd
#' qc <- initializeFlagScheme(qc, "WHP CTD")
#' qc <- initializeFlags(qc, "salinity", 2)
#' Sspan <- diff(range(qc[["SA"]]))
#' Tspan <- diff(range(qc[["CT"]]))
#' n <- length(qc[["SA"]])
#' par(mfrow=c(1, 1))
#' plotTS(qc, type="o")
#' message("Click on bad points; quit by clicking to right of plot")
#' for (i in seq_len(n)) {
#' xy <- locator(1)
#' if (xy$x > par("usr")[2])
#' break
#' i <- which.min(abs(qc[["SA"]] - xy$x)/Sspan + abs(qc[["CT"]] - xy$y)/Tspan)
#' qc <- setFlags(qc, "salinity", i=i, value=4)
#' qc <- handleFlags(qc, flags=list(salinity=4))
#' plotTS(qc, type="o")
#' }
#'}
#'
#' @family things related to ctd data
setMethod("setFlags",
c(object="ctd", name="ANY", i="ANY", value="ANY", debug="ANY"),
function(object, name=NULL, i=NULL, value=NULL, debug=getOption("oceDebug")) {
oceDebug(debug, "setFlags,ctd-method name=", name, ", i, value=", value, "\n")
if (is.null(i))
stop("must supply i")
if (!is.vector(i) && !is.function(i))
stop("'i' must be a vector or a function returning a vector")
res <- setFlagsInternal(object, name, i, value, debug-1)
res
})
#' @templateVar class ctd
#' @templateVar details {NA}
#' @template initializeFlagSchemeTemplate
setMethod("initializeFlagScheme",
signature=c(object="ctd", name="ANY", mapping="ANY", default="ANY", debug="ANY"),
definition=function(object, name=NULL, mapping=NULL, default=NULL, debug=0) {
if (is.null(name))
stop("must supply 'name'")
invisible(callNextMethod())
})
#' Initialize storage for a ctd object
#'
#' This function creates [ctd-class] objects. It is mainly
#' used by `oce` functions such as [read.ctd()] and [as.ctd()],
#' and it is not intended for novice users, so it may change at any time, without
#' following the usual rules for transitioning to deprecated and defunct status
#' (see [oce-deprecated]).
#'
#' @details
#' To save storage, this function has arguments only for quantities that are often present in data
#' files all cases. For example, not
#' all data files will have oxygen, so that's not present here.
#' Extra data may be added after the object is created, using
#' [oceSetData()].
#' Similarly, [oceSetMetadata()] may be used to add metadata (station ID, etc),
#' while bearing in mind that other functions look for such information
#' in very particular places (e.g. the station ID is a string named `station`
#' within the `metadata` slot). See [ctd-class] for more information
#' on elements stored in `ctd` objects.
#'
#' @param .Object the string `"ctd"`
#' @param pressure optional numerical vector of pressures.
#' @param salinity optional numerical vector of salinities.
#' @param temperature optional numerical vector of temperatures.
#' @param conductivity optional numerical vector of conductivities.
#' @param units optional list indicating units for the quantities specified
#' in the previous arguments. If this
#' is not supplied, a default is set up, based on which of the
#' `pressure` to `conductivity` arguments were specified.
#' If all of those 4 arguments were specified, then `units` is set
#' up as if the call included the following:
#' \code{units=list(temperature=list(unit=expression(degree*C), scale="ITS-90"),
#' salinity=list(unit=expression(), scale="PSS-78"),
#' conductivity=list(unit=expression(), scale=""),
#' pressure=list(unit=expression(dbar), scale=""),
#' depth=list(unit=expression(m), scale=""))}. This list is trimmed
#' of any of the 4 items that were not specified in the previous
#' arguments. Note that if `units` is specified, then it is just
#' copied into the `metadata` slot of the returned object, so the user
#' must be careful to set up values that will make sense to other `oce`
#' functions.
#' @param pressureType optional character string indicating the type of pressure;
#' if not supplied, this defaults to `"sea"`, which indicates the excess of
#' pressure over the atmospheric value, in dbar.
#' @param deploymentType optional character string indicating the type of deployment, which may
#' be `"unknown"`, `"profile"`, `"towyo"`, or `"thermosalinograph"`.
#' If this is not set, the value defaults to `"unknown"`.
#'
#' @family things related to ctd data
#'
#' @examples
#'
#' ## 1. empty
#' new("ctd")
#'
#' ## 2. fake data with no location information, so can only
#' ## plot with the UNESCO equation of state.
#' ## NOTE: always name arguments, in case the default order gets changed
#' ctd <- new("ctd", salinity=35+1:3/10, temperature=10-1:3/10, pressure=1:3)
#' summary(ctd)
#' plot(ctd, eos="unesco")
#'
#' ## 3. as 2, but insert location and plot with GSW equation of state.
#' ctd <- oceSetMetadata(ctd, "latitude", 44)
#' ctd <- oceSetMetadata(ctd, "longitude", -63)
#' plot(ctd, eos="gsw")
#'
#' @aliases initialize,ctd-method
setMethod(f="initialize",
signature="ctd",
definition=function(.Object, pressure, salinity, temperature, conductivity,
units,
pressureType, deploymentType) {
## Assign to some columns so they exist if needed later (even if they are NULL)
.Object@data$pressure <- if (missing(pressure)) NULL else pressure
.Object@data$temperature <- if (missing(temperature)) NULL else temperature
.Object@data$salinity <- if (missing(salinity)) NULL else salinity
.Object@data$conductivity <- if (missing(conductivity)) NULL else conductivity
names <- names(.Object@data)
##.Object@metadata$names <- names
##.Object@metadata$labels <- titleCase(names) # paste(toupper(substring(names,1,1)), substring(names,2),sep="")
##.Object@metadata$filename <- filename
if (missing(units)) {
.Object@metadata$units <- list()
if (!missing(pressure))
.Object@metadata$units$pressure <- list(unit=expression(dbar), scale="")
if (!missing(salinity))
.Object@metadata$units$salinity <- list(unit=expression(), scale="PSS-78")
if (!missing(temperature))
.Object@metadata$units$temperature <- list(unit=expression(degree*C), scale="ITS-90")
if (!missing(conductivity))
.Object@metadata$units$conductivity <- list(unit=expression(), scale="")
} else {
.Object@metadata$units <- units # CAUTION: we are being quite trusting here
}
.Object@metadata$pressureType <- if (!missing(pressureType)) pressureType else "sea" # guess on the unit
.Object@metadata$deploymentType <- if (!missing(deploymentType)) deploymentType
else "unknown" # "profile" "mooring" "towyo" "thermosalinograph"
.Object@metadata$waterDepth <- NA
#.Object@metadata$latitude <- NA
#.Object@metadata$longitude <- NA
#.Object@metadata$waterDepth <- NA
.Object@processingLog$time <- presentTime()
.Object@processingLog$value <- "create 'ctd' object"
return(.Object)
})
#' Summarize a CTD Object
#'
#' Summarizes some of the data in a `ctd` object, presenting such information
#' as the station name, sampling location, data ranges, etc. If the object was read
#' from a `.cnv` file or a `.rsk` file, then the `OriginalName`
#' column for the data summary will contain the original names of data within
#' the source file.
#'
#' @param object a [ctd-class] object.
#'
#' @param ... Further arguments passed to or from other methods.
#'
#' @examples
#' library(oce)
#' data(ctd)
#' summary(ctd)
#'
#' @author Dan Kelley
#'
#' @family things related to ctd data
setMethod(f="summary",
signature="ctd",
definition=function(object, ...) {
##mnames <- names(object@metadata)
cat("CTD Summary\n-----------\n\n")
type <- object@metadata$type
model <- object@metadata$model
mnames <- names(object@metadata)
if (!is.null(type) && nchar(type)) {
if (is.null(model)) {
cat("* Instrument: ", type, "\n")
} else {
cat("* Instrument: ", type, model, "\n")
}
}
##showMetadataItem(object, "type", "Instrument: ")
##showMetadataItem(object, "model", "Instrument model: ")
##showMetadataItem(object, "serialNumber", "Instr. serial no.: ")
showMetadataItem(object, "serialNumberTemperature", "Temp. serial no.: ")
showMetadataItem(object, "serialNumberConductivity", "Cond. serial no.: ")
showMetadataItem(object, "filename", "File: ", quote=TRUE)
showMetadataItem(object, "hexfilename", "Original file: ")
showMetadataItem(object, "institute", "Institute: ")
showMetadataItem(object, "scientist", "Chief scientist: ")
##showMetadataItem(object, "date", "Date: ", isdate=TRUE)
showMetadataItem(object, "startTime", "Start time: ", isdate=TRUE)
##showMetadataItem(object, "systemUploadTime", "System upload time: ", isdate=TRUE)
if ("sampleInterval" %in% mnames && "sampleIntervalUnits" %in% mnames)
cat("* Sample interval: ",
object@metadata$sampleInterval, " ", object@metadata$sampleIntervalUnits, "\n", sep="")
showMetadataItem(object, "cruise", "Cruise: ")
showMetadataItem(object, "ship", "Vessel: ")
showMetadataItem(object, "station", "Station: ")
deploymentType <- object@metadata$deploymentType
if (!is.null(deploymentType) && deploymentType != "unknown")
showMetadataItem(object, "deploymentType", "Deployment type: ")
if ("longitude" %in% names(object@data)) {
cat("* Mean location: ", latlonFormat(mean(object@data$latitude, na.rm=TRUE),
mean(object@data$longitude, na.rm=TRUE),
digits=5), "\n")
} else if ("longitude" %in% names(object@metadata) && !is.na(object@metadata$longitude)) {
cat("* Location: ", latlonFormat(object@metadata$latitude,
object@metadata$longitude,
digits=5), "\n")
}
showMetadataItem(object, "waterDepth", "Water depth: ")
showMetadataItem(object, "levels", "Number of levels: ")
names <- names(object@data)
invisible(callNextMethod()) # summary
})
#' @title Extract Something From a CTD Object
#'
#' @param x a [ctd-class] object.
#'
#' @examples
#' data(ctd)
#' head(ctd[["temperature"]])
#'
#' @template sub_subTemplate
#'
#' @section Details of the specialized `ctd` method:
#'
#' Some uses of \code{\link{[[,ctd-method}} involve direct retrieval of
#' items within the `data` slot of the `ctd` object,
#' while other uses involve calculations based on items in that
#' `data` slot. For an example, all `ctd` objects
#' should hold an item named `temperature` in the `data`
#' slot, so for example `x[["temperature"]]` will retrieve that
#' item. By contrast, `x[["sigmaTheta"]]` is taken to be a
#' request to compute \eqn{\sigma_\theta}{sigma[theta]}, and so
#' it yields a call to [swTheta]`(x)` even if
#' the `data` slot of `x` might happen to contain an item
#' named `theta`. This can be confusing at first, but it tends
#' to lead to fewer surprises in everyday work, for otherwise the
#' user would be forced to check the contents of any `ctd`
#' object under analysis, to determine whether that item will be looked
#' up or computed. Nothing is lost in this scheme, since the data
#' within the object are always accessible with [oceGetData()].
#'
#' It should be noted that the accessor is set up to retrieve quantities
#' in conventional units. For example, [read.ctd.sbe()] is
#' used on a `.cnv` file that stores pressure in psi, it will
#' be stored in the same unit within the `ctd` object, but
#' `x[["pressure"]]` will return a value that has been converted
#' to decibars. (To get pressure in PSI, use `x[["pressurePSI"]]`.)
#' Similarly, temperature is
#' returned in the ITS-90 scale, with a conversion having been performed with
#' [T90fromT68()], if the object holds temperature in
#' IPTS-68. Again, temperature on the IPTS-68
#' scale is returned with `x@@data$temperature`.
#'
#' This preference for computed over stored quantities is accomplished
#' by first checking for computed quantities, and then falling
#' back to the general `[[` method if no match is found.
#'
#' Some quantities are optionally computed. For example, some data files
#' (e.g. the one upon which the [section()] dataset is based)
#' store `nitrite` along with the sum of nitrite and nitrate, the
#' latter with name ``NO2+NO3``. In this case, e.g. `x[["nitrate"]]`
#' will detect the setup, and subtract nitrite from the sum to yield
#' nitrate.
#'
#' The list given below provides notes on some quantities that are,
#' or may be, computed.
#'
#' * `conductivity` without a second argument (e.g. `a[["conductivity"]]`)
#' returns the value stored in the object. However, if a second argument is given,
#' and it is string specifying a unit, then conversion is made to that unit. The
#' permitted units are: either `""` or `"ratio"` (for ratio),
#' `"uS/cm"`, `"mS/cm"` and `"S/m"`. The calculations are based on
#' the definition of conductivity ratio as the ratio between measured conductivity
#' and the standard value 4.2914 S/m.
#'
#' * `CT` or `Conservative Temperature`: Conservative Temperature,
#' computed with [gsw::gsw_CT_from_t()].
#'
#' * `density`: seawater density, computed with [swRho]`(x)`.
#' (Note that it may be better to call that function directly, to gain
#' control of the choice of equation of state, etc.)
#'
#' * `depth`: Depth in metres below the surface, computed
#' with [swDepth]`(x)`.
#'
#' * `N2`: Square of Brunt-Vaisala frequency, computed with [swN2]`(x)`.
#'
#' * `potential temperature`: Potential temperature in the
#' UNESCO formulation, computed with [swTheta]`(x)`.
#' This is a synonym for `theta`.
#'
#' * `Rrho`: Density ratio, computed with [swRrho]`(x)`.
#'
#' * `SA` or `Absolute Salinity`: Absolute Salinity,
#' computed with [gsw::gsw_SA_from_SP()].
#' The calculation involves location as well as measured water properties.
#' If the object `x` does not containing information on the location,
#' then 30N and 60W is used for the calculation, and a warning is generated.
#'
#' * `sigmaTheta`: A form of potential density anomaly, computed with
#' [swSigmaTheta]`(x)`.
#'
#' * `sigma0` Equal to `sigmaTheta`, i.e. potential density anomaly
#' referenced to a pressure of 0dbar, computed with [swSigma0]`(x)`.
#'
#' * `sigma1`: Potential density anomaly
#' referenced to a pressure of 1000dbar, computed with [swSigma1]`(x)`.
#'
#' * `sigma2`: Potential density anomaly
#' referenced to a pressure of 2000dbar, computed with [swSigma2]`(x)`.
#'
#' * `sigma3`: Potential density anomaly
#' referenced to a pressure of 3000dbar, computed with [swSigma3]`(x)`.
#'
#' * `sigma4`: potential density anomaly
#' referenced to a pressure of 4000dbar, computed with [swSigma4]`(x)`.
#'
#' * `SP`: Salinity on the Practical Salinity Scale, which is
#' `salinity` in the `data` slot.
#'
#' * `spice` or `spiciness0`: a variable that is in some sense
#' orthogonal to density, calculated with [swSpice]`(x)`.
#' Note that this is defined differently for `eos="unesco"` and
#' `eos="gsw"`.
#'
#' * `SR`: Reference Salinity computed with [gsw::gsw_SR_from_SP()].
#'
#' * `Sstar`: Preformed Salinity computed with [gsw::gsw_SR_from_SP()].
#' See `SA` for a note on longitude and latitude.
#'
#' * `theta`: potential temperature in the UNESCO formulation,
#' computed with [swTheta]`(x)`. This is a synonym for
#' `potential temperature`.
#'
#' * `time`: returns either a vector of times, a single
#' time, or `NULL`. A vector is returned if `time`
#' is present in the `data` slot, or if a time can be
#' inferred from other entries in the `data` slot (some of which,
#' such as the common `timeS`, also employ
#' `startTime` within the `metadata` slot). A single
#' value is returned if the dataset only has information on the start
#' time (which is stored as `startTime` within the `metadata`
#' slot. If it is impossible to determine the sampling time, then
#' `NULL` is returned. These time variants occur, in the
#' present version of oce, only for data read by [read.ctd.sbe()],
#' the documention of which explains how times are computed.
#'
#' * `z`: Vertical coordinate in metres above the surface, computed with
#' [swZ]`(x)`.
#'
#'
#' @author Dan Kelley
#'
#' @family things related to ctd data
setMethod(f="[[",
signature(x="ctd", i="ANY", j="ANY"),
definition=function(x, i, j, ...) {
data <- x@data
metadata <- x@metadata
dataNames <- names(data)
metadataNames <- names(metadata)
## message("i=\"", i, "\"")
if (i == "conductivity") {
C <- data$conductivity
##message("i=", i, ", j=", if (missing(j)) "(missing)" else j)
if (!is.null(C) && !missing(j)) {
if (!(j %in% c("", "ratio", "uS/cm", "mS/cm", "S/m")))
stop("unknown conductivity unit \"", j, "\"; must be \"\", \"ratio\", \"uS/cm\", \"mS/cm\" or \"S/m\"")
if (j == "")
j <- "ratio" # lets us use switch()
unit <- metadata$units$conductivity$unit
if (is.null(unit) || !length(unit)) {
## FIXME: maybe should look at median value, to make a guess
## warning("ctd object lack conductivity units; assuming \"ratio\"")
unit <- "ratio"
}
##message("A")
unit <- as.character(unit)
##message("next is unit:")
##print(dput(unit))
C <- data$conductivity
##message("B")
## Rather than convert from 3 inputs to 3 outputs, express as ratio, then convert as desired
if (!unit %in% c("ratio", "uS/cm", "mS/cm", "S/m"))
stop("object has unknown conductivity unit \"", unit, "\"; must be \"ratio\", \"uS/cm\", \"mS/cm\" or \"S/m\"")
C <- C / switch(unit, "uS/cm"=42914, "mS/cm"=42.914, "S/m"=4.2914, "ratio"=1)
C <- C * switch(j, "uS/cm"=42914, "mS/cm"=42.914, "S/m"=4.2914, "ratio"=1)
}
C
} else if (i == "salinity" || i == "SP") {
if ("salinity" %in% dataNames) {
S <- data$salinity
} else {
C <- data$conductivity
if (!is.null(C)) {
if (is.null(metadata$units$conductivity)) {
warning("conductivity has no unit, so guessing it is conductivity-ratio. Be cautious on calculated salinity.")
} else {
unit <- as.character(metadata$units$conductivity$unit)
if (0 == length(unit)) {
S <- swSCTp(C, x[["temperature"]], x[["pressure"]])
warning("constructed salinity from temperature, conductivity-ratio and pressure")
} else if (unit == "uS/cm") {
S <- swSCTp(C/42914.0, x[["temperature"]], x[["pressure"]])
warning("constructed salinity from temperature, conductivity and pressure")
} else if (unit == "mS/cm") {
## e.g. RSK
S <- swSCTp(C/42.914, x[["temperature"]], x[["pressure"]])
warning("constructed salinity from temperature, conductivity and pressure")
} else if (unit == "S/m") {
S <- swSCTp(C/4.2914, x[["temperature"]], x[["pressure"]])
warning("constructed salinity from temperature, conductivity and pressure")
} else {
stop("unrecognized conductivity unit '", unit, "'; only uS/cm, mS/cm and S/m are handled")
}
}
} else {
stop("the object's data slot lacks 'salinity', and it cannot be calculated since 'conductivity' is also missing")
}
}
S
} else if (i == "SR") {
gsw::gsw_SR_from_SP(SP=x[["salinity"]])
} else if (i == "Sstar") {
if (!any(is.finite(x[["longitude"]])) || !any(is.finite(x[["latitude"]])))
stop("object lacks location information, so Sstar cannot be computed")
n <- length(data$salinity)
## Lengthen lon and lat if necessary, by repeating.
lon <- metadata$longitude
if (n != length(lon))
lon <- rep(metadata$longitude, length.out=n)
lat <- metadata$latitude
if (n != length(lat))
lat <- rep(metadata$latitude, length.out=n)
lon <- ifelse(lon < 0, lon + 360, lon) # not required because gsw_saar() does this ... but UNDOCUMENTED
## Do the calculation in two steps
SA <- gsw::gsw_SA_from_SP(SP=x[["salinity"]], p=x[["pressure"]], longitude=lon, latitude=lat)
gsw::gsw_Sstar_from_SA(SA=SA, p=x[["pressure"]], longitude=lon, latitude=lat)
} else if (i == "temperature") {
scale <- metadata$units[["temperature"]]$scale
if (!is.null(scale) && "IPTS-68" == scale)
T90fromT68(data$temperature)
else data$temperature
} else if (i == "pressure") {
if ("pressure" %in% dataNames) {
pressure <- data$pressure
unit <- metadata$units[["pressure"]]$unit
## NOTE: 2019-04-29: The next will always return pressure, from the
## else part of the conditional. This is because oce
## stores pressure as dbar, and copies any original PSI data
## into data$pressurePSI.
if (!is.null(unit) && "psi" == as.character(unit))
pressure * 0.6894757 # 1 psi=6894.757 Pa
else pressure
} else {
if ("depth" %in% dataNames)
swPressure(data$depth)
else stop("object's data slot does not contain 'pressure' or 'depth'")
}
## } else if (i == "longitude") {
## if ("longitude" %in% metadataNames) metadata$longitude else data$longitude
## } else if (i == "latitude") {
## if ("latitude" %in% metadataNames) metadata$latitude else data$latitude
} else if (i == "time") {
## After checking for 'time' literally in the metadata
## and data slots, we turn to the 10 time variants
## listed in the SBE Seasoft data processing manual,
## revision 7.26.8, Appendix VI.
if ("time" %in% dataNames) {
data$time
} else if ("timeS" %in% dataNames) {
if (!("startTime" %in% metadataNames)) {
warning("have timeS in data slot, but no startTime in metadata slot, so cannot compute [[\"time\"]]\n")
NULL
} else {
metadata$startTime + data$timeS
}
} else if ("timeM" %in% dataNames) {
if (!("startTime" %in% metadataNames)) {
warning("have timeM in data slot, but no startTime in metadata slot, so cannot compute [[\"time\"]]\n")
NULL
} else {
metadata$startTime + 60 * data$timeM
}
} else if ("timeH" %in% dataNames) {
if (!("startTime" %in% metadataNames)) {
warning("have timeH in data slot, but no startTime in metadata slot, so cannot compute [[\"time\"]]\n")
NULL
} else {
metadata$startTime + 3600 * data$timeH
}
} else if ("timeJ" %in% dataNames) {
if (!("startTime" %in% metadataNames)) {
warning("have timeJ in data slot, but no startTime in metadata slot, so cannot compute [[\"time\"]]\n")
NULL
} else {
t0 <- ISOdatetime(1900 + as.POSIXlt(metadata$startTime)$year, 1, 1, 0, 0, 0, tz="UTC")
timeJ <- data$timeJ
## Handle wraparound (FIXME: this code is tricky and not well-tested, esp wrt leap year)
wraps <- sum(diff(timeJ) < 0)
n <- length(timeJ)
for (w in rev(which(diff(timeJ)<0))) {
timeJ[seq(w+1, n)] <- round(timeJ[w]) + timeJ[seq(w+1, n)]
}
t0 + 86400 * (timeJ - 1)
}
} else if ("timeN" %in% dataNames) {
as.POSIXct(data$timeN, origin="1970-01-01",tz="UTC")
} else if ("timeQ" %in% dataNames) {
as.POSIXct(data$timeQ, origin="2000-01-01",tz="UTC")
} else if ("timeK" %in% dataNames) {
as.POSIXct(data$timeK, origin="2000-01-01",tz="UTC")
} else if ("timeJV2" %in% dataNames) {
if (!("startTime" %in% metadataNames)) {
warning("have timeJV2 in data slot, but no startTime in metadata slot, so cannot compute [[\"time\"]]\n")
NULL
} else {
t0 <- ISOdatetime(1900 + as.POSIXlt(metadata$startTime)$year, 1, 1, 0, 0, 0, tz="UTC")
timeJV2 <- data$timeJV2
## Handle wraparound (FIXME: this code is tricky and not well-tested, esp wrt leap year)
wraps <- sum(diff(timeJV2) < 0)
n <- length(timeJV2)
for (w in rev(which(diff(timeJV2)<0))) {
timeJV2[seq(w+1, n)] <- round(timeJV2[w]) + timeJV2[seq(w+1, n)]
}
t0 + 86400 * (timeJV2 - 1)
}
} else if ("timeSCP" %in% dataNames) {
if (!("startTime" %in% metadataNames)) {
warning("have timeSCP in data slot, but no startTime in metadata slot, so cannot compute [[\"time\"]]\n")
NULL
} else {
t0 <- ISOdatetime(1900 + as.POSIXlt(metadata$startTime)$year, 1, 1, 0, 0, 0, tz="UTC")
timeSCP <- data$timeSCP
## Handle wraparound (FIXME: this code is tricky and not well-tested, esp wrt leap year)
wraps <- sum(diff(timeSCP) < 0)
n <- length(timeSCP)
for (w in rev(which(diff(timeSCP)<0))) {
timeSCP[seq(w+1, n)] <- round(timeSCP[w]) + timeSCP[seq(w+1, n)]
}
t0 + 86400 * (timeSCP - 1)
}
} else if ("timeY" %in% dataNames) {
as.POSIXct(data$timeY, origin="1970-01-01",tz="UTC")
} else if ("time" %in% metadataNames) {
metadata$time
} else if ("startTime" %in% metadataNames) {
metadata$startTime
} else {
NULL
}
## end of time decoding (whew!)
} else if (i == "N2") {
swN2(x)
} else if (i == "density") {
swRho(x)
} else if (i == "sigmaTheta") {
swSigmaTheta(x)
} else if (i == "sigma0") {
swSigma0(x)
} else if (i == "sigma1") {
swSigma1(x)
} else if (i == "sigma2") {
swSigma2(x)
} else if (i == "sigma3") {
swSigma3(x)
} else if (i == "sigma4") {
swSigma4(x)
} else if (i %in% c("theta", "potential temperature")) {
swTheta(x)
} else if (i == "Rrho") {
swRrho(x)
} else if (i == "spice" || i == "spiciness") {
swSpice(x)
} else if (i %in% c("absolute salinity", "SA")) {
if (!any(is.finite(x[["longitude"]])) || !any(is.finite(x[["latitude"]])))
stop("object lacks location information, so SA cannot be computed")
SP <- x[["salinity"]]
p <- x[["pressure"]]
n <- length(SP)
## Lengthen lon and lat if necessary, by repeating.
lon <- x[["longitude"]]
if (n != length(lon))
lon <- rep(lon, length.out=n)
lat <- x[["latitude"]]
if (n != length(lat))
lat <- rep(lat, length.out=n)
lon <- ifelse(lon < 0, lon + 360, lon) # not required because gsw_saar() does this ... but UNDOCUMENTED
##: Change e.g. NaN to NA ... FIXME: tests show that this is not required:
##: > a<-as.ctd(10:11, c(35, asin(3)), 1:2, lon=-60, lat=50)
##: a[["SA"]]
##: [1] 10.0472934071578 11.0520223880037
##: > a<-as.ctd(10:11, c(35, NA), 1:2, lon=-60, lat=50)
##: > a[["SA"]]
##: [1] 10.0472934071578 11.0520223880037
##: SP[!is.finite(SP)] <- NA
##: p[!is.finite(p)] <- NA
##: lon[!is.finite(lon)] <- NA
##: lat[!is.finite(lat)] <- NA
gsw::gsw_SA_from_SP(SP, p, lon, lat)
} else if (i %in% c("conservative temperature", "CT")) {
if (!any(is.finite(x[["longitude"]])) || !any(is.finite(x[["latitude"]])))
stop("object lacks location information, so CT cannot be computed")
gsw::gsw_CT_from_t(SA=x[["SA"]], t=x[["temperature"]], p=x[["pressure"]])
} else if (i == "nitrate") {
if ("nitrate" %in% dataNames) {
data$nitrate
} else {
if ("nitrite" %in% dataNames && "NO2+NO3" %in% dataNames)
data[["NO2+NO3"]] - data$nitrite
else NULL
}
} else if (i == "nitrite") {
if ("nitrite" %in% dataNames) {
data$nitrite
} else {
if ("nitrate" %in% dataNames && "NO2+NO3" %in% dataNames)
data[["NO2+NO3"]] - data$nitrate
else NULL
}
} else if (i == "z") {
swZ(x) # FIXME-gsw: permit gsw version here
} else if (i == "depth") {
if ("depth" %in% names(data)) data$depth else swDepth(x) # FIXME-gsw: permit gsw version here
} else if (i == "N2") {
swN2(x)
} else {
callNextMethod() # [[
}
})
#' @title Replace Parts of a CTD Object
#'
#' @param x a [ctd-class] object.
#'
#' @template sub_subsetTemplate
#'
#' @examples
#' data(ctd)
#' summary(ctd)
#' # Move the CTD profile a nautical mile north.
#' ctd[["latitude"]] <- 1/60 + ctd[["latitude"]] # acts in metadata
#' # Increase the salinity by 0.01.
#' ctd[["salinity"]] <- 0.01 + ctd[["salinity"]] # acts in data
#' summary(ctd)
#'
#' @family things related to ctd data
setMethod(f="[[<-",
signature(x="ctd", i="ANY", j="ANY"),
definition=function(x, i, j, ..., value) {
callNextMethod(x=x, i=i, j=j, ..., value=value) # [[<-
})
#' Coerce data into CTD object
#'
#' Assemble data into a [ctd-class] object.
#'
#' @param salinity There are several distinct choices for `salinity`.
#' * It can be an [rsk-class] object (see \dQuote{Converting rsk objects} for details).
#'
#' * It can be a
#' vector indicating the practical salinity through the water column. In that case,
#' `as.ctd` employs the other arguments listed below.
#' * It can be an [rsk-class] object (see \dQuote{Converting rsk objects} for details).
#'
#' * It can be something (a data frame, a list or an `oce` object)
#' from which practical salinity, temperature, pressure, and conductivity
#' can be inferred. In this case, the relevant information
#' is extracted and the other arguments to `as.ctd` are ignored, except for
#' `pressureAtmospheric`. If the first argument has salinity, etc., in
#' matrix form (as can happen with some objects of [argo-class]),
#' then only the first column is used, and a warning to that effect is given,
#' unless the `profile` argument is specified and then that specific
#' profile is extracted.
#' * It can be an [rsk-class] object (see \dQuote{Converting rsk objects} for details).
#'
#' * It can be an [rsk-class] object (see \dQuote{Converting rsk objects} for details).
#'
#' * It can be unspecified, in which case `conductivity` becomes a mandatory
#' argument, because it will be needed for computing actual salinity,
#' using [swSCTp()].
#'
#' @param temperature *in-situ* temperature in \eqn{^\circ deg}C on
#' the ITS-90 scale; see \dQuote{Temperature units} in the documentation for
#' [swRho()].
#'
#' @param pressure Vector of pressure values, one for each `salinity` and
#' `temperature` pair, or just a single pressure, which is repeated to match
#' the length of `salinity`.
#'
#' @param conductivity electrical conductivity ratio through the water column
#' (optional). To convert from raw conductivity in milliSeimens per centimeter
#' divide by 42.914 to get conductivity ratio (see Culkin and Smith, 1980).
#'
##1108 @param SA absolute salinity (as in TEOS-10). If given, the supplied absolute
##1108 salinity is converted internally to UNESCO-defined practical salinity.
##1108
##1108 @param CT conservative temperature (as in TEOS-10). If given, the supplied
##1108 conservative temperature is converted internally to UNESCO-defined in-situ
##1108 temperature.
##1108
##1108 @param oxygen optional oxygen concentration
##1108
##1108 @param nitrate optional nitrate concentration
##1108
##1108 @param nitrite optional nitrite concentration
##1108
##1108 @param phosphate optional phosphate concentration
##1108
##1108 @param silicate optional silicate concentration
##1108
#' @param scan optional scan number. If not provided, this will be set to