forked from dankelley/oce
-
Notifications
You must be signed in to change notification settings - Fork 0
/
section.R
2916 lines (2866 loc) · 147 KB
/
section.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 Hydrographic Section Data
#'
#' This class stores data from oceanographic section surveys.
#'
#' Sections can be read with \code{\link{read.section}} or created with
#' \code{\link{read.section}} or created from CTD objects by using
#' \code{\link{as.section}} or by adding a ctd station to an existing section with
#' \code{\link{sectionAddStation}}.
#'
#' Sections may be sorted with \code{\link{sectionSort}}, subsetted with
#' \code{\link{subset,section-method}}, smoothed with \code{\link{sectionSmooth}}, and
#' gridded with \code{\link{sectionGrid}}. Gridded sections may be plotted with
#' \code{\link{plot,section-method}}.
#'
#' Statistical summaries are provided by \code{\link{summary,section-method}}, while
#' overviews are provided by \code{show}.
#'
#' The sample dataset \code{\link{section}} contains data along WOCE line A03.
#'
#' @templateVar class section
#'
#' @templateVar dataExample {}
#'
#' @templateVar metadataExample Examples that are of common interest include \code{stationId}, \code{longitude}, \code{latitude} and \code{time}.
#'
#' @template slot_summary
#'
#' @template slot_put
#'
#' @template slot_get
#'
#' @examples
#' library(oce)
#' data(section)
#' plot(section[['station', 1]])
#' pairs(cbind(z=-section[["pressure"]],T=section[["temperature"]],S=section[["salinity"]]))
#' ## T profiles for first few stations in section, at common scale
#' par(mfrow=c(3,3))
#' Tlim <- range(section[["temperature"]])
#' ylim <- rev(range(section[["pressure"]]))
#' for (stn in section[["station", 1:9]])
#' plotProfile(stn, xtype='potential temperature', ylim=ylim, Tlim=Tlim)
#'
#' @author Dan Kelley
#'
#' @family classes provided by \code{oce}
#' @family things related to \code{section} data
setClass("section", contains="oce")
#' @title Hydrographic section
#'
#' @description
#' This is line A03 (ExpoCode 90CT40_1, with nominal sampling date 1993-09-11).
#' The chief scientist was Tereschenkov of SOI, working aboard the Russian ship
#' Multanovsky, undertaking a westward transect from the Mediterranean outflow
#' region across to North America, with a change of heading in the last few dozen
#' stations to run across the nominal Gulf Stream axis.
#' The data flags follow the WHP CTD convention, i.e. 1 for uncalibrated,
#' 2 for an acceptable measurement, 3 for a questionable measurement, 4
#' for a bad measurement, etc; see \url{https://www.nodc.noaa.gov/woce/woce_v3/wocedata_1/whp/exchange/exchange_format_desc.htm}
#' for further details.
#'
#' @examples
#' \dontrun{
#' library(oce)
#' # Gulf Stream
#' data(section)
#' GS <- subset(section, 109<=stationId&stationId<=129)
#' GSg <- sectionGrid(GS, p=seq(0, 5000, 100))
#' plot(GSg, map.xlim=c(-80,-60))
#' }
#'
#' @name section
#'
#' @docType data
#'
#' @usage data(section)
#'
#' @source This is based on the WOCE file named \code{a03_hy1.csv}, downloaded
#' from \url{https://cchdo.ucsd.edu/cruise/90CT40_1}, 13 April 2015.
#'
#' @family datasets provided with \code{oce}
#' @family things related to \code{section} data
NULL
setMethod(f="initialize",
signature="section",
definition=function(.Object, filename="", sectionId="") {
.Object@metadata$filename <- filename
.Object@metadata$sectionId <- sectionId
.Object@processingLog$time <- as.POSIXct(Sys.time())
.Object@processingLog$value <- "create 'section' object"
return(.Object)
})
## 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 Section Objects
#' @details
#' If \code{flags} and \code{actions} are not provided, the
#' default is to use WHP (World Hydrographic Program) flags [1], in which the
#' value 2 indicates good data, and other values indicate either unchecked,
#' suspicious, or bad data. Any data not flagged as good are set
#' to \code{NA} in the returned value. Since WHP flag codes run
#' from 1 to 9, this default is equivalent to
#' setting \code{flags=list(all=c(1, 3:9))} along with
#' \code{action=list("NA")}.
#' @param object An object of \code{\link{section-class}}.
#' @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)
#' section2 <- handleFlags(section)
#' par(mfrow=c(2, 1))
#' plotTS(section)
#' plotTS(section2)
#'
#' @family things related to \code{section} data
setMethod("handleFlags",
c(object="section", flags="ANY", actions="ANY", debug="ANY"),
function(object, flags=list(), actions=list(), debug=integer()) {
## DEVELOPER 1: alter the next comment to explain your setup
## Default to the World Hydrographic Program system, with
## flags from 1 to 9, with flag=2 for acceptable data.
if (missing(flags))
flags <- list(c(1, 3:9)) # DEVELOPER 2: alter this line to suit a newdata class
if (missing(actions)) {
actions <- list("NA") # DEVELOPER 3: alter this line to suit a new data class
names(actions) <- names(flags)
}
if (missing(debug))
debug <- getOption("oceDebug")
if (any(names(actions)!=names(flags))) {
stop("names of flags and actions must match")
}
res <- object
for (i in seq_along(res@data$station)) {
res@data$station[[i]] <- handleFlags(res@data$station[[i]], flags, actions, debug)
}
res
})
#' @title Summarize a Section Object
#'
#' @description
#' Pertinent summary information is presented, including station locations,
#' distance along track, etc.
#'
#' @param object An object of class \code{"section"}, usually, a result of a call
#' to \code{\link{read.section}}, \code{\link{read.oce}}, or
#' \code{\link{as.section}}.
#'
#' @param ... Further arguments passed to or from other methods.
#'
#' @return \code{NULL}
#'
#' @examples
#' library(oce)
#' data(section)
#' summary(section)
#'
#' @family things related to \code{section} data
#'
#' @author Dan Kelley
setMethod(f="summary",
signature="section",
definition=function(object, ...) {
numStations <- length(object@data$station)
cat("Section Summary\n---------------\n\n")
cat("* Source: \"", object@metadata$filename, "\"\n", sep="")
cat("* ID: \"", object@metadata$sectionId, "\"\n", sep="")
##stn.sum <- matrix(nrow=numStations, ncol=5)
if (numStations > 0) {
cat("Overview of stations\n```\n")
cat(sprintf("%5s %5s %8s %8s %7s %5s\n", "Index", "ID", "Lon", "Lat", "Levels", "Depth"))
for (i in 1:numStations) {
##stn <- object@data$station[[i]]
thisStn <- object@data$station[[i]]
id <- if (!is.null(thisStn@metadata$station) && "" != thisStn@metadata$station)
thisStn@metadata$station else ""
depth <- if (!is.finite(thisStn@metadata$waterDepth) || 0 == thisStn@metadata$waterDepth)
max(thisStn@data$pressure, na.rm=TRUE) else thisStn@metadata$waterDepth
cat(sprintf("%5d %5s %8.3f %8.3f %7.0f %5.0f\n",
i, id,
thisStn@metadata$longitude, thisStn@metadata$latitude,
length(thisStn@data$pressure), depth))
}
cat("```\n")
} else {
cat("* No stations\n")
}
processingLogShow(object)
invisible()
})
#' @title Extract Something From a Section Object
#' @param x A \code{section} object, i.e. one inheriting from \code{\link{section-class}}.
#'
#' @templateVar class section
#'
#' @section Details of the specialized \code{section} method:
#'
#' There are several possibilities, depending on the nature of \code{i}.
#'\itemize{
#' \item If \code{i} is the string \code{"station"}, then the method
#' will return a \code{\link{list}} of
#' \code{\link{ctd-class}} objects holding the station data.
#' If \code{j} is also given, it specifies a station (or set of stations) to be returned.
#' if \code{j} contains just a single value, then that station is returned, but otherwise
#' a list is returned. If \code{j} is an integer, then the stations are specified by index,
#' but if it is character, then stations are specified by the names stored within
#' their metadata. (Missing stations yield \code{NULL} in the return value.)
#'
#' \item If \code{i} is \code{"station ID"}, then the IDs of the stations in the
#' section are returned.
#'
#' \item If \code{i} is \code{"dynamic height"}, then an estimate of dynamic
#' height is returned (as calculated with \code{\link{swDynamicHeight}(x)}).
#'
#' \item If \code{i} is \code{"distance"}, then the distance along the section is
#' returned, using \code{\link{geodDist}}.
#'
#' \item If \code{i} is \code{"depth"}, then a vector containing the depths
#' of the stations is returned.
#'
#' \item If \code{i} is \code{"z"}, then a vector containing the z
#' coordinates is returned.
#'
#' \item If \code{i} is \code{"theta"} or \code{"potential temperature"}, then
#' the potential temperatures of all the stations are returned in one
#' vector. Similarly, \code{"spice"} returns the property known
#' as spice, using \code{\link{swSpice}}.
#'
#' \item If \code{i} is a string ending with \code{"Flag"}, then the characters
#' prior to that ending are taken to be the name of a variable contained
#' within the stations in the section. If this flag is available in
#' the first station of the section, then the flag values are looked
#' up for every station.
#'}
#'
#' If \code{j} is \code{"byStation"}, then a list is returned, with
#' one (unnamed) item per station.
## #' If \code{j} is \code{"grid:distance-pressure"}, then a gridded
## #' representation of \code{i} is returned, as a list with elements
## #' \code{distance} (in km), \code{pressure} (in dbar) and
## #' \code{field} (in whatever unit is used for \code{i}). See Example
## #' for in the documentation for \code{\link{plot,section-method}}.
#'
#' @template sub_subTemplate
#'
#' @examples
#' data(section)
#' length(section[["latitude"]])
#' length(section[["latitude", "byStation"]])
#' # Vector of all salinities, for all stations
#' Sv <- section[["salinity"]]
#' # List of salinities, grouped by station
#' Sl <- section[["salinity", "byStation"]]
#' # First station salinities
#' Sl[[1]]
#'
#' @family things related to \code{section} data
#' @author Dan Kelley
setMethod(f="[[",
signature(x="section", i="ANY", j="ANY"),
definition=function(x, i, j, ...) {
## Data-quality flags are a special case
res <- NULL
if (1 == length(grep(".*Flag$", i))) {
baseName <- gsub("Flag$", "", i)
if (baseName %in% names(x@data$station[[1]]@metadata$flags)) {
res <- unlist(lapply(x@data$station, function(ctd) ctd[[i]]))
return(res)
} else {
stop("the stations within this section do not contain a '", baseName, "' flag")
}
}
nstation <- length(x@data$station)
## some derived things (not all ... be sure to document when adding things!)
##20160809 if (i %in% c("theta", "potential temperature", "sigmaTheta")) {
##20160809 res <- unlist(lapply(x@data$station, function(ctd) ctd[[i]]))
##20160809 return(res)
##20160809 }
if (i %in% names(x@metadata)) {
if (i %in% c("longitude", "latitude")) {
if (!missing(j) && j == "byStation") {
return(x@metadata[[i]])
} else {
res <- NULL
for (stn in seq_along(x@data$station))
res <- c(res, rep(x@data$station[[stn]]@metadata[[i]], length(x@data$station[[stn]][["salinity"]])))
return(res)
}
} else {
return(x@metadata[[i]])
}
} else if (i %in% c("absolute salinity", "CT", "conservative temperature",
"density", "depth", "nitrite", "nitrate",
"potential temperature", "SA", "sigmaTheta",
"spice", "theta", "z",
names(x@data$station[[1]]@data))) {
##message('section.R:311 section[["', i, '"]]')
if (!missing(j) && substr(j, 1, 4) == "grid") {
if (j == "grid:distance-pressure") {
numStations <- length(x@data$station)
p1 <- x[["station", 1]][["pressure"]]
np1 <- length(p1)
field <- matrix(NA, nrow=numStations, ncol=np1)
if (numStations > 1) {
field[1, ] <- x[["station", 1]][[i]]
for (istn in 2:numStations) {
pi <- x[["station", istn]][["pressure"]]
if (length(pi) != np1 || any(pi != p1)) {
warning("returning NULL because this section is not gridded")
return(NULL)
}
field[istn, ] <- x[["station", istn]][[i]]
}
res <- list(distance=x[["distance", "byStation"]], pressure=p1, field=field)
return(res)
} else {
warning("returning NULL because this section contains only 1 station")
return(NULL)
}
} else {
warning("returning NULL because only grid:distance-pressure is permitted")
return(NULL)
}
} else {
## Note that nitrite and nitrate might be computed, not stored
if (!missing(j) && j == "byStation") {
##message("section[['", i, "', 'byStation']] START")
res <- vector("list", nstation)
for (istation in seq_len(nstation))
res[[istation]] <- x@data$station[[istation]][[i]]
##message("section[['", i, "', 'byStation']] END")
} else {
##message("section[['", i, "']] START")
res <- NULL
for (station in x[["station"]])
res <- c(res, station[[i]])
##message("section[['", i, "']] END")
}
return(res)
}
} else if (i == "station") {
if (missing(j)) {
res <- x@data$station
} else {
if (is.character(j)) {
nj <- length(j)
stationNames <- unlist(lapply(x@data$station, function(x) x@metadata$station))
if (nj == 1) {
w <- which(stationNames == j)
res <- if (length(w)) x@data$station[[w[1]]] else NULL
} else {
res <- vector("list", nj)
for (jj in j) {
w <- which(stationNames == j)
res[[jj]] <- if (length(w)) x@data$station[[w[1]]] else NULL
}
}
} else {
nj <- length(j)
if (nj == 1) {
res <- x@data$station[[j]]
} else {
res <- vector("list", nj)
for (jj in j)
res[[jj]] <- x@data$station[[jj]]
}
}
}
} else if ("station ID" == i) {
res <- NULL
for (stn in x[['station']])
res <- c(res, stn[['station']])
} else if ("dynamic height" == i) {
res <- swDynamicHeight(x)
} else if ("distance" == i) {
res <- NULL
for (stn in seq_along(x@data$station)) {
distance <- geodDist(x@data$station[[stn]]@metadata$longitude,
x@data$station[[stn]]@metadata$latitude,
x@data$station[[1]]@metadata$longitude,
x@data$station[[1]]@metadata$latitude)
if (!missing(j) && j == "byStation")
res <- c(res, distance)
else
res <- c(res, rep(distance, length(x@data$station[[stn]]@data$temperature)))
}
} else if ("time" == i) {
## time is not in the overall metadata ... look in the individual objects
res <- unlist(lapply(x@data$station, function(stn) stn[["time"]]))
res <- numberAsPOSIXct(res)
} else {
res <- callNextMethod() # [[
}
res
})
#' @title Replace Parts of a Section Object
#' @param x A \code{section} object, i.e. inheriting from \code{\link{section-class}}
#' @family things related to \code{section} data
#' @template sub_subsetTemplate
#' @examples
#' # 1. Change section ID from a03 to A03
#' data(section)
#' section[["sectionId"]]
#' section[["sectionId"]] <- toupper(section[["sectionId"]])
#' section[["sectionId"]]
#' # 2. Add a millidegree to temperatures at station 10
#' section[["station", 10]][["temperature"]] <-
#' 1e-3 + section[["station", 10]][["temperature"]]
#' @author Dan Kelley
setMethod(f="[[<-",
signature(x="section", i="ANY", j="ANY"),
definition=function(x, i, j, ..., value) {
if (i == "station") {
if (missing(j)) {
x@data$station <- value
} else {
x@data$station[[j]] <- value
}
x
} else {
callNextMethod(x=x, i=i, j=j, ...=..., value=value) # [[<-
}
})
setMethod(f="show",
signature="section",
definition=function(object) {
id <- object@metadata$sectionId
n <- length(object@data$station)
if (n == 0) {
cat("Section has no stations\n")
} else {
if (is.null(id) || id == "")
cat("Unnamed section has ", n, " stations:\n", sep="")
else
cat("Section '", id, "' has ", n, " stations:\n", sep="")
cat(sprintf("%5s %5s %8s %8s %7s %5s\n", "Index", "ID", "Lon", "Lat", "Levels", "Depth"))
##cat(sprintf("%4s %5s %10.2f %10.2f %10.0f\n", "Index", "ID", "Lon", "Lat", "Depth\n"))
for (i in 1:n) {
thisStn <- object@data$station[[i]]
id <- if (!is.null(thisStn@metadata$station) && "" != thisStn@metadata$station)
thisStn@metadata$station else ""
depth <- if (is.null(thisStn@metadata$waterDepth))
max(thisStn@data$pressure, na.rm=TRUE) else thisStn@metadata$waterDepth
cat(sprintf("%5d %5s %8.3f %8.3f %7.0f %5.0f\n",
i, id,
thisStn@metadata$longitude, thisStn@metadata$latitude,
length(thisStn@data$pressure), depth))
}
}
})
#' Subset a Section Object
#'
#' Return a subset of a section object.
#'
#' This function is used to subset data within the
#' stations of a section, or to choose a subset of the stations
#' themselves. The first case is handled with the \code{subset} argument,
#' while the second is handled if \code{...} contains a vector named
#' \code{indices}. Either \code{subset} or \code{indices} must be
#' provided, but not both.
#'
#' \strong{In the "subset" method}, \code{subset} indicates
#' either stations to be kept, or data to be kept within the stations.
#'
#' The first step in processing is to check for the presence of certain
#' key words in the \code{subset} expression. If \code{distance}
#' is present, then stations are selected according to a condition on the
#' distance (in km) from the first station to the given station (Example 1).
#' If either \code{longitude} or \code{latitude} is given, then
#' stations are selected according to the stated condition (Example 2).
#' If \code{stationId} is present, then selection is in terms of the
#' station ID (not the sequence number) is used (Example 3). In all of these
#' cases, stations are either selected in their entirety or dropped
#' in their entirety.
#'
#' If none of these keywords is present, then the \code{subset} expression is
#' evaluated in the context of the \code{data} slot of each of the CTD stations
#' stored within \code{x}. (Note that this slot does not normally
#' contain any of the keywords that are listed in the previous
#' paragraph; it does, then odd results may occur.) Each station is examined
#' in turn, with \code{subset} being evaluated individually in each. The evaluation
#' produces a logical vector. If that vector has length 1 (Examples 4 and 5)
#' then the station is retained or discarded, accordingly. If the vector is longer,
#' then the logical vector is used as a sieve to subsample that individual CTD
#' profile.
#'
#' \strong{In the "indices" method}, stations are selected using
#' \code{indices}, which may be a vector of integers that indicate sequence
#' number, or a logical vector, again indicating which stations to keep.
#'
#' @param x A \code{section} object, i.e. one inheriting from \code{\link{section-class}}.
#' @param subset an optional indication of either the stations to be kept,
#' or the data to be kept within the stations. See \dQuote{Details}.
#'
#' @param ... optional arguments, of which only the first is examined. The only
#' possibility is that this argument be named \code{indices}. See \dQuote{Details}.
#'
#' @return A \code{\link{section-class}} object.
#'
#' @examples
#' library(oce)
#' data(section)
#'
#' # 1. Stations within 500 km of the first station
#' starting <- subset(section, distance < 500)
#'
#' # 2. Stations east of 50W
#' east <- subset(section, longitude > (-50))
#'
#' # 3. Gulf Stream
#' GS <- subset(section, 109 <= stationId & stationId <= 129)
#'
#' # 4. Only stations with more than 5 pressure levels
#' long <- subset(section, length(pressure) > 5)
#'
#' # 5. Only stations that have some data in top 50 dbar
#' surfacing <- subset(section, min(pressure) < 50)
#'
#' # 6. Similar to #4, but done in more detailed way
#' long <- subset(section,
#' indices=unlist(lapply(section[["station"]],
#' function(s)
#' 5 < length(s[["pressure"]]))))
#'
#' @family functions that subset \code{oce} objects
#' @family things related to \code{section} data
#'
#' @author Dan Kelley
setMethod(f="subset",
signature="section",
definition=function(x, subset, ...) {
subsetString <- paste(deparse(substitute(subset)), collapse=" ")
res <- x
dots <- list(...)
dotsNames <- names(dots)
indicesGiven <- length(dots) && ("indices" %in% dotsNames)
debug <- getOption("oceDebug")
if (length(dots) && ("debug" %in% names(dots)))
debug <- dots$debug
if (indicesGiven) {
## select a portion of the stations
if (!missing(subset))
stop("cannot specify both 'subset' and 'indices'")
oceDebug(debug, "subsetting by indices\n")
res <- new("section")
indices <- dots$indices
n <- length(indices)
if (is.logical(indices))
indices <- (1:n)[indices]
if (min(indices) < 1) stop("cannot have negative indices")
if (max(indices) > length(x@data$station)) stop("cannot indices exceeding # stations")
stn <- x@metadata$stationId[indices]
lat <- x@metadata$lat[indices]
lon <- x@metadata$lon[indices]
station <- vector("list", length(indices))
for (i in seq_along(indices)) {
station[[i]] <- x@data$station[[indices[i]]]
}
data <- list(station=station)
res@metadata$stationId <- stn
res@metadata$longitude <- lon
res@metadata$latitude <- lat
res@data <- data
res@processingLog <- x@processingLog
res@processingLog <- processingLogAppend(res@processingLog, paste("subset(x, indices=c(", paste(dots$indices, collapse=","), "))", sep=""))
} else {
if (missing(subset))
stop("must give 'subset' or (in ...) 'indices'")
oceDebug(debug, "subsetting by 'subset'\n")
##subsetString <- deparse(substitute(subset))
##oceDebug(debug, "subsetString='", subsetString, "'\n")
res <- x
if (length(grep("stationId", subsetString))) {
keep <- eval(substitute(subset),
envir=data.frame(stationId=as.numeric(x@metadata$stationId)))
res@metadata$stationId <- x@metadata$stationId[keep]
res@metadata$longitude <- x@metadata$longitude[keep]
res@metadata$latitude <- x@metadata$latitude[keep]
res@metadata$time <- x@metadata$time[keep]
res@data$station <- x@data$station[keep]
res@processingLog <- processingLogAppend(res@processingLog, paste("subset(x, subset=", subsetString, ")", sep=""))
} else if (length(grep("distance", subsetString))) {
l <- list(distance=geodDist(res))
keep <- eval(substitute(subset), l, parent.frame(2))
res@metadata$longitude <- res@metadata$longitude[keep]
res@metadata$latitude <- res@metadata$latitude[keep]
res@metadata$stationId <- res@metadata$stationId[keep]
res@data$station <- res@data$station[keep]
} else if (length(grep("levels", subsetString))) {
levels <- unlist(lapply(x[["station"]], function(stn) length(stn[["pressure"]])))
keep <- eval(substitute(subset), list(levels=levels))
res@metadata$longitude <- res@metadata$longitude[keep]
res@metadata$latitude <- res@metadata$latitude[keep]
res@metadata$stationId <- res@metadata$stationId[keep]
res@data$station <- res@data$station[keep]
} else if (length(grep("latitude", subsetString)) || length(grep("longitude", subsetString))) {
n <- length(x@data$station)
keep <- vector(length=n)
for (i in 1:n)
keep[i] <- eval(substitute(subset), x@data$station[[i]]@metadata, parent.frame(2))
nn <- sum(keep)
station <- vector("list", nn)
stn <- vector("character", nn)
lon <- vector("numeric", nn)
lat <- vector("numeric", nn)
j <- 1
for (i in 1:n) {
if (keep[i]) {
stn[j] <- x@metadata$stationId[i]
lon[j] <- x@metadata$longitude[i]
lat[j] <- x@metadata$latitude[i]
station[[j]] <- x@data$station[[i]]
j <- j + 1
}
}
res <- new("section")
res@data$station <- station
res@metadata$header <- x@metadata$header
res@metadata$sectionId <- x@metadata$sectionId
res@metadata$stationId <- stn
res@metadata$longitude <- lon
res@metadata$latitude <- lat
res@processingLog <- x@processingLog
} else {
res <- new("section")
res@data$station <- list()
res@metadata$header <- x@metadata$header
res@metadata$sectionId <- NULL
res@metadata$stationId <- NULL
res@metadata$longitude <- NULL
res@metadata$latitude <- NULL
res@processingLog <- x@processingLog
n <- length(x@data$station)
j <- 1
for (i in 1:n) {
r <- eval(substitute(subset), x@data$station[[i]]@data, parent.frame(2))
oceDebug(debug, "i=", i, ", j=", j, ", sum(r)=", sum(r), "\n", sep="")
if (sum(r) > 0) {
## copy whole station ...
res@data$station[[j]] <- x@data$station[[i]]
## ... but if we are looking for a subset, go through the data fields and do that
if (length(r) > 1) {
## Select certain levels. This occurs e.g. for
## subset(sec, S > 35)
## but not for station-by-station selection, e.g. as a result of
## subset(sec, length(S) > 3)
## since the length of the latter is 1, which means to copy
## the whole station.
for (field in names(res@data$station[[j]]@data)) {
oceDebug(debug, " field='", field, "', i=", i, ", j=", j, " (case A)\n", sep="")
res@data$station[[j]]@data[[field]] <- res@data$station[[j]]@data[[field]][r]
}
}
## copy section metadata
res@metadata$stationId[j] <- x@metadata$stationId[i]
res@metadata$latitude[j] <- x@metadata$latitude[i]
res@metadata$longitude[j] <- x@metadata$longitude[i]
j <- j + 1
} else {
oceDebug(debug, " skipping this station\n")
}
}
## if (j <= n) {
## for (jj in seq.int(n, j)) {
## oceDebug(debug, "erase item at j =", jj, "\n")
## res@data$station[[jj]] <- NULL
## res@metadata$stationId[jj] <- TRUE
## res@metadata$latitude[jj] <- TRUE
## res@metadata$longitude[jj] <- TRUE
## }
## }
}
res@processingLog <- processingLogAppend(res@processingLog, paste("subset(x, subset=", subsetString, ")", sep=""))
}
res
})
#' @title Sort a Section
#'
#' @description
#' Sections created with \code{\link{as.section}} have "stations" that are in the
#' order of the CTD objects (or filenames for such objects) provided. Sometimes,
#' this is not the desired order, e.g. if file names discovered with
#' \code{\link{dir}} are in an order that makes no sense. (For example, a
#' practioner might name stations \code{"stn1"}, \code{"stn2"}, etc., not
#' realizing that this will yield an unhelpful ordering, by file name, if there
#' are more than 9 stations.) The purpose of \code{sectionSort} is to permit
#' reordering the constituent stations in sensible ways.
#'
#' @param section A \code{section} object containing the section whose stations
#' are to be sorted.
#'
#' @param by An optional string indicating how to reorder. If not provided,
#' \code{"stationID"} will be assumed. Other choices are \code{"distance"}, for
#' distance from the first station, \code{"longitude"}, for longitude,
#' \code{"latitude"} for latitude, and \code{"time"}, for time.
#'
#' @return An object of \code{\link{section-class}} that has less lateral
#' variation than the input section.
#'
#' @examples
#' \dontrun{
#' # Eastern North Atlantic, showing Mediterranean water;
#' # sorting by longitude makes it easier to compare
#' # the map and the section.
#' library(oce)
#' data(section)
#' s <- sectionGrid(subset(section, -30 <= longitude))
#' ss <- sectionSort(ss, by="longitude")
#' plot(ss)
#' }
#'
#' @author Dan Kelley
#'
#' @family things related to \code{section} data
sectionSort <- function(section, by)
{
if (missing(by)) {
by <- "stationId"
} else {
byChoices <- c("stationId", "distance", "longitude", "latitude", "time")
iby <- pmatch(by, byChoices, nomatch=0)
if (0 == iby)
stop('unknown by value "', by, '"; should be one of: ', paste(byChoices, collapse=" "))
by <- byChoices[iby]
}
res <- section
if (by == "stationId") {
o <- order(section@metadata$stationId)
} else if (by == "distance") {
o <- order(section[["distance", "byStation"]])
} else if (by == "longitude") {
o <- order(section[["longitude", "byStation"]])
} else if (by == "latitude") {
o <- order(section[["latitude", "byStation"]])
} else if (by == "time") {
## FIXME: should check to see if startTime exists first?
times <- unlist(lapply(section@data$station, function(x) x@metadata$startTime))
o <- order(times)
} else {
o <- seq_along(section[["station"]]) ## cannot ever get here, actually
}
res@metadata$stationId <- res@metadata$stationId[o]
res@metadata$longitude <- res@metadata$longitude[o]
res@metadata$latitude <- res@metadata$latitude[o]
if ("time" %in% names(res@metadata))
res@metadata$time <- res@metadata$time[o]
res@data$station <- res@data$station[o]
res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep="", collapse=""))
res
}
#' @title Make a Section (DEFUNCT)
#'
#' @description
#' This is a defunct function; use \code{\link{as.section}} instead, and
#' see \link{oce-defunct} for more on the oce procedure for retiring functions.
#' @param item Ignored, since this function is defunct
#' @param ... Ignored, since this function is defunct
makeSection <- function(item, ...)
{
.Defunct("as.section",
msg="makeSection() was marked 'defunct' in March 2016, after having been marked 'deprecated' for a CRAN release cycle. Use as.section() instead. See ?'oce-defunct'.")
## .Deprecated("as.section",
## msg="makeSection() will be removed soon; use as.section() instead. See ?'oce-deprecated'.")
## if (inherits(item, "ctd")) {
## extra.args <- list(...)
## numStations <- 1 + length(extra.args)
## station <- vector("list", numStations)
## stn <- vector("character", numStations)
## lon <- vector("numeric", numStations)
## lat <- vector("numeric", numStations)
## stn[1] <- item@metadata$station
## lon[1] <- item@metadata$longitude
## lat[1] <- item@metadata$latitude
## station[[1]] <- item
## if (numStations > 1) {
## for (i in 2:numStations) {
## ## message("DAN ", i)
## thisStn <- extra.args[[i-1]]
## stn[i] <- thisStn@metadata$station
## lon[i] <- thisStn@metadata$longitude
## lat[i] <- thisStn@metadata$latitude
## station[[i]] <- thisStn
## }
## }
## } else if (inherits(item, "list") && !inherits(item[[1]], "oce")) {
## stop("cannot yet handle a list of non-oce objects")
## } else if (inherits(item, "list") && inherits(item[[1]], "oce")) {
## numStations <- length(item)
## station <- vector("list", numStations)
## stn <- vector("character", numStations)
## lon <- vector("numeric", numStations)
## lat <- vector("numeric", numStations)
## if (numStations < 1)
## stop("need more than 1 item in the list, to create a section")
## ## 2015-12-06 if (inherits(item[[1]], "oce")) {
## for (i in 1:numStations) {
## if (!inherits(item[[i]], "oce"))
## stop("list cannot be a mixture of oce and non-oce items")
## thisItem <- item[[i]]
## stn[i] <- if (is.null(thisItem@metadata$station)) i else thisItem@metadata$station
## lon[i] <- thisItem@metadata$longitude
## lat[i] <- thisItem@metadata$latitude
## station[[i]] <- thisItem
## }
## ## 2015-12-06 } else {
## ## 2015-12-06: this code block could not be run
## ## 2015-12-06 ## demand that items contain @data$pressure
## ## 2015-12-06 if ("pressure" %in% names(item[[1]]) || "pressure" %in% names(item[[1]]@data)) {
## ## 2015-12-06 stop("items must contain pressure")
## ## 2015-12-06 }
## ## 2015-12-06 for (thisItem in item) {
## ## 2015-12-06 names <- names(thisItem)
## ## 2015-12-06 if (!("longitude" %in% names)) stop("each item entry must contain longitude")
## ## 2015-12-06 if (!("latitude" %in% names)) stop("each item must entry contain latitude")
## ## 2015-12-06 ## FIXME: maybe permits 'depth' here
## ## 2015-12-06 if (!("pressure" %in% names)) stop("each item must entry contain pressure")
## ## 2015-12-06 if (!("station" %in% names)) thisItem$station <- seq_along(thisItem$longitude)
## ## 2015-12-06 len <- length(thisItem$pressure)
## ## 2015-12-06 print(names)
## ## 2015-12-06 names <- names[names!="longitude"]
## ## 2015-12-06 names <- names[names!="latitude"]
## ## 2015-12-06 names <- names[names!="station"]
## ## 2015-12-06 print(names)
## ## 2015-12-06 data <- list()
## ## 2015-12-06 for (name in names) {
## ## 2015-12-06 if (length(name) == len) {
## ## 2015-12-06 data[[name]] <- thisItem[[name]]
## ## 2015-12-06 }
## ## 2015-12-06 }
## ## 2015-12-06 str(data)
## ## 2015-12-06 }
## ## 2015-12-06 }
## } else if (class(item) == "character") {
## numStations <- length(item)
## station <- vector("list", numStations)
## stn <- vector("character", numStations)
## lon <- vector("numeric", numStations)
## lat <- vector("numeric", numStations)
## if (numStations <= 1)
## stop("need more than 1 station to make a section")
## if (exists(item[1])) {
## ## ctd objects
## ##oceDebug(1, "ctd objects\n")
## for (i in 1:numStations) {
## thisItem <- get(item[[i]])
## stn[i] <- thisItem@metadata$station
## lon[i] <- thisItem@metadata$longitude
## lat[i] <- thisItem@metadata$latitude
## station[[i]] <- thisItem
## }
## } else {
## ## ctd filenames
## ##oceDebug(1, "ctd files\n")
## for (i in 1:numStations) {
## ##oceDebug(1, "file named", item[i], "\n")
## ctd <- read.ctd(item[i])
## stn[i] <- ctd@metadata$station
## lon[i] <- ctd@metadata$longitude
## lat[i] <- ctd@metadata$latitude
## station[[i]] <- ctd
## }
## }
## } else {
## stop("first argument must be a \"ctd\" object, a \"list\" of ctd objects, or a vector of character strings naming ctd objects")
## }
## res <- new("section")
## res@metadata$sectionId <- ""
## res@metadata$stationId <- stn
## res@metadata$longitude <- lon
## res@metadata$latitude <- lat
## res@data <- list(station=station)
## res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep="", collapse=""))
## res
}
#' @title Add a CTD Station to a Section
#'
#' @description
#' Add a CTD profile to an existing section.
#'
#' @section Historical note:
#' Until March 2015, this operation was carried out with the \code{+} operator,
#' but at that time, the syntax was flagged by the development version of R, so it
#' was changed to the present form.
#'
#' @param section A section to which a station is to be added.
#'
#' @param station A ctd object holding data for the station to be added.
#'
#' @aliases sectionAddCtd
#' @return An object of \code{\link[base]{class}} \code{section}.
#'
#' @examples
#' library(oce)
#' data(ctd)
#' ctd2 <- ctd
#' ctd2[["temperature"]] <- ctd2[["temperature"]] + 0.5
#' ctd2[["latitude"]] <- ctd2[["latitude"]] + 0.1
#' section <- as.section(c("ctd", "ctd2"))
#' ctd3 <- ctd
#' ctd3[["temperature"]] <- ctd[["temperature"]] + 1
#' ctd3[["latitude"]] <- ctd[["latitude"]] + 0.1
#' ctd3[["station"]] <- "Stn 3"
#' sectionAddStation(section, ctd3)
#'
#' @author Dan Kelley
#'
#' @family things related to \code{section} data
sectionAddStation <- function(section, station)
{
if (missing(section)) stop("must provide a section to which the ctd is to be added")
if (!inherits(section, "section")) stop("'section' is not a 'section' object")
if (missing(station)) return(section)
if (!inherits(station, "ctd")) stop("'station' is not a 'ctd' object")
res <- section
n.orig <- length(section@data$station)
s <- vector("list", n.orig + 1)
for (i in 1:n.orig)
s[[i]] <- section@data$station[[i]]
s[[n.orig + 1]] <- station
res@data$station <- s
res@metadata$longitude <- c(res@metadata$longitude, station@metadata$longitude)
res@metadata$latitude <- c(res@metadata$latitude, station@metadata$latitude)
res@metadata$stationId <- c(res@metadata$stationId, station@metadata$station)
res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep="", collapse=""))
res
}
sectionAddCtd <- sectionAddStation
#' @title Plot a Section
#'
#' @description
#' Creates a summary plot for a CTD section, with one panel for each value of
#' \code{which}.
#'
#' @details
#' The type of plot is governed by \code{which}, as follows.
#'
#' \itemize{
#' \item \code{which=0} or \code{"potential temperature"} for temperature contours
#' \item \code{which=1} or \code{"temperature"} for temperature contours (the default)
#' \item \code{which=2} or \code{"salinity"} for salinity contours
#' \item \code{which=3} or \code{"sigmaTheta"} for sigma-theta contours
#' \item \code{which=4} or \code{"nitrate"} for nitrate concentration contours
#' \item \code{which=5} or \code{"nitrite"} for nitrite concentration contours
#' \item \code{which=6} or \code{"oxygen"} for oxygen concentration contours
#' \item \code{which=7} or \code{"phosphate"} for phosphate concentration contours
#' \item \code{which=8} or \code{"silicate"} for silicate concentration contours
#' \item \code{which=9} or \code{"u"} for eastward velocity
#' \item \code{which=10} or \code{"uz"} for vertical derivative of eastward velocity
#' \item \code{which=11} or \code{"v"} for northward velocity
#' \item \code{which=12} or \code{"vz"} for vertical derivative of northward velocity
#' \item \code{which=20} or \code{"data"} for a dot for each data location
#' \item \code{which=99} or \code{"map"} for a location map
#' }
#'
#' The y-axis for the contours is pressure, plotted in the conventional reversed
#' form, so that the water surface appears at the top of the plot. The x-axis is
#' more complicated. If \code{at} is not supplied, then the routine calculates x
#' as the distance between the first station in the section and each of the other
#' stations. (This will produce an error if the stations are not ordered
#' geographically, because the \code{\link{contour}} routine cannot handle
#' non-increasing axis coordinates.) If \code{at} is specified, then it is taken
#' to be the location, in arbitrary units, along the x-axis of labels specified by
#' \code{labels}; the way this works is designed to be the same as for
#' \code{\link{axis}}.
#'
#'
#' @param x a \code{section} object, e.g. as created by \code{\link{as.section}}
#' or \code{\link{read.section}}.
#'
#' @param which a list of desired plot types, as explained in \dQuote{Details}.
#' There may be up to four panels in total, and the desired plots are placed in
#' these panels, in reading order. If only one panel is plotted, \code{par} is
#' not adjusted, which makes it easy to add to the plot with subsequent plotting
#' commands.
#'