forked from dankelley/oce
-
Notifications
You must be signed in to change notification settings - Fork 0
/
misc.R
5291 lines (5132 loc) · 201 KB
/
misc.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=80:expandtab:shiftwidth=4:softtabstop=4
abbreviateVector <- function(x)
{
if (1 >= length(x)) {
return(x)
} else {
ud <- unique(diff(x))
if (1 == length(ud) && 1 == ud) return(paste(x[1], ":", tail(x, 1), sep="")) else return(x)
}
}
#' Convert angles from 0:360 to -180:180
#'
#' This is mostly used for instrument heading angles, in cases where the
#' instrument is aligned nearly northward, so that small variations in heading
#' (e.g. due to mooring motion) can yield values that swing from small angles
#' to large angles, because of the modulo-360 cut point.
#' The method is to use the cosine and sine of the angle in order to find "x"
#' and "y" values on a unit circle, and then to use [atan2()] to
#' infer the angles.
#'
#' @param theta an angle (in degrees) that is in the range from 0 to 360
#' degrees
#' @return A vector of angles, in the range -180 to 180.
#' @author Dan Kelley
#' @examples
#'
#' library(oce)
#' ## fake some heading data that lie near due-north (0 degrees)
#' n <- 20
#' heading <- 360 + rnorm(n, sd=10)
#' heading <- ifelse(heading > 360, heading - 360, heading)
#' x <- 1:n
#' plot(x, heading, ylim=c(-10, 360), type='l', col='lightgray', lwd=10)
#' lines(x, angleRemap(heading))
angleRemap <- function(theta)
{
toRad <- atan2(1, 1) / 45
atan2(sin(toRad * theta), cos(toRad * theta)) / toRad
}
#' Earth magnetic declination
#'
#' Instruments that use magnetic compasses to determine current direction need
#' to have corrections applied for magnetic declination, to get currents with
#' the y component oriented to geographic, not magnetic, north. Sometimes, and
#' for some instruments, the declination is specified when the instrument is
#' set up, so that the velocities as recorded are already. Other times, the
#' data need to be adjusted. This function is for the latter case.
#'
#' @param x an [oce-class] object.
#'
#' @param declination magnetic declination (to be added to the heading)
#'
#' @param debug a debugging flag, set to a positive value to get debugging.
#'
#' @return Object, with velocity components adjusted to be aligned with
#' geographic north and east.
#'
#' @author Dan Kelley
#'
#' @seealso Use [magneticField()] to determine the declination,
#' inclination and intensity at a given spot on the world, at a given time.
#'
#' @references
#' 1. \samp{https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html}
#'
#' @family things related to magnetism
applyMagneticDeclination <- function(x, declination=0, debug=getOption("oceDebug"))
{
oceDebug(debug, "applyMagneticDeclination(x,declination=", declination, ") {\n", sep="", unindent=1)
if (inherits(x, "cm")) {
oceDebug(debug, "object is of type 'cm'\n")
res <- x
S <- sin(-declination * pi / 180)
C <- cos(-declination * pi / 180)
r <- matrix(c(C, S, -S, C), nrow=2)
uvr <- r %*% rbind(x@data$u, x@data$v)
res@data$u <- uvr[1, ]
res@data$v <- uvr[2, ]
oceDebug(debug, "originally, first u:", x@data$u[1:3], "\n")
oceDebug(debug, "originally, first v:", x@data$v[1:3], "\n")
oceDebug(debug, "after application, first u:", res@data$u[1:3], "\n")
oceDebug(debug, "after application, first v:", res@data$v[1:3], "\n")
} else {
stop("cannot apply declination to object of class ", paste(class(x), collapse=", "), "\n")
}
res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep="", collapse=""))
oceDebug(debug, "} # applyMagneticDeclination\n", unindent=1)
res
}
#' Trilinear interpolation in a 3D array
#'
#' Interpolate within a 3D array, using the trilinear approximation.
#'
#' Trilinear interpolation is used to interpolate within the `f` array,
#' for those (`xout`, `yout` and `zout`) triplets that are
#' inside the region specified by `x`, `y` and `z`. Triplets
#' that lie outside the range of `x`, `y` or `z` result in
#' `NA` values.
#'
#' @param x vector of x values for grid (must be equi-spaced)
#' @param y vector of y values for grid (must be equi-spaced)
#' @param z vector of z values for grid (must be equi-spaced)
#' @param f matrix of rank 3, with the gridded values mapping to the `x`
#' values (first index of `f`), etc.
#' @param xout vector of x values for output.
#' @param yout vector of y values for output (length must match that of
#' `xout`).
#' @param zout vector of z values for output (length must match that of
#' `xout`).
#' @return A vector of interpolated values (or `NA` values), with length
#' matching that of `xout`.
#' @author Dan Kelley and Clark Richards
#'
#' @examples
#' ## set up a grid
#' library(oce)
#' n <- 5
#' x <- seq(0, 1, length.out=n)
#' y <- seq(0, 1, length.out=n)
#' z <- seq(0, 1, length.out=n)
#' f <- array(1:n^3, dim=c(length(x), length(y), length(z)))
#' ## interpolate along a diagonal line
#' m <- 100
#' xout <- seq(0, 1, length.out=m)
#' yout <- seq(0, 1, length.out=m)
#' zout <- seq(0, 1, length.out=m)
#' approx <- approx3d(x, y, z, f, xout, yout, zout)
#' ## graph the results
#' plot(xout, approx, type='l')
#' points(xout[1], f[1, 1, 1])
#' points(xout[m], f[n,n,n])
approx3d <- function(x, y, z, f, xout, yout, zout)
{
## Were all arguments given?
if (missing(x))
stop("must provide x")
if (missing(y))
stop("must provide y")
if (missing(z))
stop("must provide z")
if (missing(f))
stop("must provide f")
if (missing(xout))
stop("must provide xout")
if (missing(yout))
stop("must provide yout")
if (missing(zout))
stop("must provide zout")
## Are there enough data to interpolate?
if (length(x) < 2)
stop("must have more than one x value")
if (length(y) < 2)
stop("must have more than one y value")
if (length(z) < 2)
stop("must have more than one z value")
## Are the array dimensions consistent with x, y, and z?
if (3 != length(dim(f)))
stop("f must be an array with 3 dimensions")
if (length(x) != dim(f)[1])
stop("length of x and dim(f)[1] must agree, but they are ", length(x), " and ", dim(f)[1])
if (length(y) != dim(f)[2])
stop("length of y and dim(f)[2] must agree, but they are ", length(y), " and ", dim(f)[2])
if (length(z) != dim(f)[3])
stop("length of z and dim(f)[3] must agree, but they are ", length(z), " and ", dim(f)[3])
## Are x, y and z equi-spaced?
if (length(x) > 2) {
equispaced <- function(a) sd(diff(a)) / mean(diff(a)) < 1e-5
if (!equispaced(x))
stop("x values must be equi-spaced")
if (!equispaced(y))
stop("y values must be equi-spaced")
if (!equispaced(z))
stop("z values must be equi-spaced")
}
do_approx3d(x, y, z, f, xout, yout, zout)
}
#' Show an argument to a function, e.g. for debugging
#'
#' @param x the argument
#'
#' @param nshow number of values to show at first (if length(x)> 1)
#'
#' @param last indicates whether this is the final argument to the function
#'
#' @param sep the separator between name and value
argShow <- function(x, nshow=4, last=FALSE, sep="=")
{
if (missing(x))
return("")
name <- paste(substitute(expr=x, env=environment()))
res <- ""
nx <- length(x)
if (missing(x)) {
res <- "(missing)"
} else {
if (is.null(x)) {
res <- NULL
} else {
if (is.function(x)) {
res <- "(provided)"
} else if (nx==1) {
res <- if (is.character(x)) paste('"', x[1], '"', sep="") else x[1]
} else {
look <- seq.int(1L, min(nshow, nx)-1)
res <- paste(format(x[look], digits=4), collapse=",")
res <- paste(res, if (nx > nshow) ", ..., " else ",", x[nx], sep="")
}
}
}
res <- if (nx == 1) paste(name, "=", res, sep="") else paste(name, "=c(", res, ")", sep="")
if (!last)
res <- paste(res, ",", sep="")
res
}
#' Get first finite value in a vector or array, or NULL if none
#' @param v A numerical vector or array.
firstFinite <- function(v)
{
if (!is.vector(v))
v <- as.vector(v)
first <- which(is.finite(v))
if (length(first) > 0) v[first[1]] else NULL
}
#' Read a World Ocean Atlas NetCDF File
#'
#' @param file character string naming the file
#' @param name of variable to extract. If not provided, an
#' error message is issued that lists the names of data in the file.
#' @param positive logical value indicating whether `longitude` should
#' be converted to be in the range from 0 to 360, with `name`
#' being shuffled accordingly. This is set to `FALSE` by default,
#' because the usual oce convention is for longitude to range between -180
#' to +180.
#'
#' @return A list containing vectors `longitude`, `latitude`,
#' `depth`, and an array with the specified name. If `positive`
#' is true, then `longitude` will be converted to range from 0
#' to 360, and the array will be shuffled accordingly.
#'
#' @examples
#'\dontrun{
#' ## Mean SST at 5-degree spatial resolution
#' tmn <- read.woa("/data/woa13/woa13_decav_t00_5dv2.nc", "t_mn")
#' imagep(tmn$longitude, tmn$latitude, tmn$t_mn[,,1], zlab="SST")
#'}
read.woa <- function(file, name, positive=FALSE)
{
if (!missing(file) && is.character(file) && 0 == file.info(file)$size)
stop("empty file")
if (!is.character(file))
stop("'file' must be a character string")
con <- ncdf4::nc_open(file)
if (missing(name)) {
varnames <- names(con$var)
stop("must supply a name from the list: ", paste(varnames, collapse=", "))
return(NULL)
}
longitude <- as.vector(ncdf4::ncvar_get(con, "lon"))
latitude <- as.vector(ncdf4::ncvar_get(con, "lat"))
depth <- as.vector(ncdf4::ncvar_get(con, "depth"))
field <- ncdf4::ncvar_get(con, name)
if (positive) {
lon2 <- ifelse(longitude < 0, longitude + 360, longitude)
i <- order(lon2)
longitude <- longitude[i]
## Crude method to reorder field on first index, whether it is 2D, 3D or 4D,
## although I'm not sure that any 4D items occur in the World Ocean Atlas.
if (is.array(field)) {
ndim <- length(dim(field))
if (ndim == 2)
field <- field[i,]
else if (ndim == 3)
field <- field[i,,]
else if (ndim == 4)
field <- field[i,,,]
}
}
rval <- list(longitude=longitude, latitude=latitude, depth=depth, field=field)
names(rval) <- c(head(names(rval), -1), name)
rval
}
## alphabetized functions END
## unalphabetized functions START
shortenTimeString <- function(t, debug=getOption("oceDebug"))
{
tc <- as.character(t)
oceDebug(debug, "shortenTimeString() {\n", sep="", unindent=1)
oceDebug(debug, "A: '", paste(t, collapse="' '"), "'\n")
tc <- gsub(" [A-Z]{3}$", "", tc) # remove timezone
if (all(grepl("^[0-9]{4}", tc))) {
## leading years
years <- substr(tc, 1, 4)
if (1 == length(unique(years))) {
tc <- gsub("^[0-9]{4}", "", tc)
tc <- gsub("^-", "", tc) # works for ISO dates
oceDebug(debug, "B: '", paste(tc, collapse="' '"), "'\n", sep='')
}
} else if (any(grepl("[a-zA-Z]", tc))) {
## Change e.g. 'Jul 01' to 'Jul' if all labels end in 01
if (all(grepl("01\\s*$", tc))) {
tc <- gsub(" 01\\s*$", "", tc)
oceDebug(debug, "B: '", paste(tc, collapse="' '"), "'\n", sep='')
}
}
oceDebug(debug, "C: '", paste(tc, collapse="' '"), "'\n", sep='')
tc <- gsub("^\\s*", "", tc)
tc <- gsub("\\s*$", "", tc)
oceDebug(debug, "D: '", paste(tc, collapse="' '"), "'\n", sep='')
oceDebug(debug, "}\n", unindent=1)
tc
}
#' Convert each of a vector of strings from SNAKE_CASE to camelCase
#'
#' `snakeToCamel` converts "snake-case" characters such as `"NOVA_SCOTIA"`
#' to "camel-case" values, such as `"NovaScotia"`. It was written for
#' use by [read.argo()], but it also may prove helpful in other contexts.
#'
#' The basic procedure is to chop the string up into substrings separated by
#' the underline character, then to upper-case the first letter of
#' all substrings except the first, and then to paste the substrings
#' together.
#'
#' However, there are exceptions. First, any upper-case string that contains no
#' underlines is converted to lower case, but any mixed-case string with no
#' underlines is returned as-is (see the second example). Second, if
#' the `specialCases` argument contains `"QC"`, then the
#' `QC` is passed through directly (since it is an acronym) and
#' if the first letter of remaining text is upper-cased (contrast
#' see the four examples).
#'
#' @param s A vector of character values.
#'
#' @param specialCases A vector of character values that tell which
#' special-cases to apply, or `NULL` (the default) to turn off special
#' cases. The only permitted special case at the moment is `"QC"` (see
#' \dQuote{Details}) but the idea of this argument is that other cases
#' can be added later, if needed.
#'
#' @return A vector of character values
#'
#' @examples
#' library(oce)
#' snakeToCamel("PARAMETER_DATA_MODE") # "parameterDataMode"
#' snakeToCamel("PARAMETER") # "parameter"
#' snakeToCamel("HISTORY_QCTEST") # "historyQctest"
#' snakeToCamel("HISTORY_QCTEST", "QC") # "historyQCTest"
#' snakeToCamel("PROFILE_DOXY_QC") # "profileDoxyQc"
#' snakeToCamel("PROFILE_DOXY_QC", "QC") # "profileDoxyQC"
#' @author Dan Kelley
snakeToCamel <- function(s, specialCases=NULL)
{
ns <- length(s)
if ("QC" %in% specialCases) {
s <- gsub("QCTEST", "Q_C_TEST", s) # for e.g. HISTORY_QCTEST
s <- gsub("QC$", "Q_C", s) # for e.g. PROFILE_DOXY_QC
s <- gsub("Qc$", "QC", s) # for e.g. positionQc (converted previously)
}
if (ns < 1)
stop("'s' must be a vector of character values")
res <- vector("character", length=ns)
for (is in seq(1L, length(s))) {
if (!grepl("_", s[is])) {
## Handle the single-word case. If all upper-case, convert to lower,
## but otherwise, leave it as it is.
res[is] <- if (s[is] == toupper(s[is])) tolower(s[is]) else s[is]
} else {
## Handle the multi-word case. Start by making it lower case.
s[is] <- tolower(s[is])
## Now, split and then work through the words
w <- strsplit(s[is], "_")[[1]]
nw <- length(w)
res[is] <- w[1]
for (iw in 2:nw) {
wl <- tolower(w[iw])
res[is] <- paste0(res[is], toupper(substring(wl,1,1)), substring(wl,2))
}
}
}
res
}
#' Decode units, from strings
#'
#' This is mainly intended for internal use within the package, e.g. by
#' [read.odf()], and so the list of string-to-unit mappings is not
#' documented, since developers can learn it from simple examination
#' of the code. The focus of `unitFromString()` is on strings that are
#' found in oceanographic files available to the author, *not* on all
#' possible units.
#'
#' @param unit a character value indicating the unit. These
#' are matched according to rules developed to work with actual
#' data files, and so the list is not by any means exhaustive.
#'
#' @param scale a character value indicating the scale. The default value
#' of `NULL` dictates that the scale is to be inferred from the unit. If
#' a non-`NULL` value is supplied, it will be used, even if it makes no sense
#' in relation to value of `unit`.
#'
#' @return A [list()] of two items: `unit` which is an
#' [expression()], and `scale`, which is a string.
#'
#' @examples
#' unitFromString("dbar") # dbar (no scale)
#' unitFromString("deg c") # modern temperature (ITS-90 scale)
#' @family functions that interpret variable names and units from headers
unitFromString <- function(unit, scale=NULL)
{
if (length(unit) > 1L)
stop("cannot work with a vector of strings")
u <- trimws(unit) # remove any leading/trailing whitespace
U <- toupper(u) # simplify some match tests
#> message("unit=\"", unit, "\", scale=\"", scale, "\"")
if (U == "" || U == "NONE" || U == "(NONE)")
return(list(unit=expression(), scale=if (is.null(scale)) "" else scale))
if (U == "10**3CELLS/L")
return(list(unit=expression(10^3*cells/l), scale=if (is.null(scale)) "" else scale))
if (U == "CODE")
return(list(unit=expression(), scale=if (is.null(scale)) "" else scale))
if (U == "COUNTS")
return(list(unit=expression(), scale=if (is.null(scale)) "" else scale))
if (U == "DB") # NOTE: this really should be decibel, but ODF files use it for decibar, sometimes
return(list(unit=expression(dbar), scale=if (is.null(scale)) "" else scale))
if (U == "DBAR" || U == "DECIBAR" || U == "DECIBARS")
return(list(unit=expression(dbar), scale=if (is.null(scale)) "" else scale))
if (U == "DEG C" || U == "DEGREES C")
return(list(unit=expression(degree*C), scale=if (is.null(scale)) "ITS-90" else scale))
if (U == "FMOL/KG")
return(list(unit=expression(fmol/kg), scale=if (is.null(scale)) "" else scale))
if (U == "G")
return(list(unit=expression(g), scale=if (is.null(scale)) "" else scale))
if (U == "GMT")
return(list(unit=expression(), scale=if (is.null(scale)) "" else scale))
if (U == "HPA")
return(list(unit=expression(hPa), scale=if (is.null(scale)) "" else scale))
if (U == "HZ")
return(list(unit=expression(Hz), scale=if (is.null(scale)) "" else scale))
if (U == "ITS-90" || U == "ITS-90 DEGC")
return(list(unit=expression(degree*C), scale=if (is.null(scale)) "ITS-90" else scale))
if (U == "IPTS-68" || U == "IPTS-68 DEGC" || U == "IPTS-68, DEG C")
return(list(unit=expression(degree*C), scale=if (is.null(scale)) "IPTS-68" else scale))
if (U == "ITS-68" || U == "ITS-68 DEGC" || U == "ITS-68, DEG C") # not an accepted unit, but seen in ODF files
return(list(unit=expression(degree*C), scale=if (is.null(scale)) "IPTS-68" else scale))
if (U == "KG/M^3" || U == "KG/M**3")
return(list(unit=expression(kg/m^3), scale=if (is.null(scale)) "" else scale))
if (U == "M" || U == "METER" || U == "METRE" || U == "METERS" || U == "METRES")
return(list(unit=expression(m), scale=if (is.null(scale)) "" else scale))
if (U == "M**3/KG" || U == "M^3/KG")
return(list(unit=expression(m^3/kg), scale=if (is.null(scale)) "" else scale))
if (U == "MA")
return(list(unit=expression(ma), scale=if (is.null(scale)) "" else scale))
if (U == "M/S" || U == "METER/SEC" || U == "M/S" || U == "METRE/SEC")
return(list(unit=expression(m/s), scale=if (is.null(scale)) "" else scale))
if (U == "MG/M^3" || U == "MG/M**3")
return(list(unit=expression(mg/m^3), scale=if (is.null(scale)) "" else scale))
if (U == "MICRON" || U == "MICRONS")
return(list(unit=expression(mu*m), scale=if (is.null(scale)) "" else scale))
if (grepl("^\\s*MICRO[ ]?MOL[E]?S/M(\\*){0,2}2/S(EC)?\\s*$", U))
return(list(unit=expression(mu*mol/m/s), scale=if (is.null(scale)) "" else scale))
if (U == "ML/L")
return(list(unit=expression(ml/l), scale=if (is.null(scale)) "" else scale))
if (U == "M/S" || U == "M/SEC")
return(list(unit=expression(m/s), scale=if (is.null(scale)) "" else scale))
if (U == "M^-1/SR")
return(list(unit=expression(1/m/sr), scale=if (is.null(scale)) "" else scale))
if (U == "MHO/CM" || U == "MHOS/CM") # 1 mho (archaic) = 1 Siemen (modern)
return(list(unit=expression(S/cm), scale=if (is.null(scale)) "" else scale))
if (U == "MHO/M" || U == "MHOS/M") # 1 mho (archaic) = 1 Siemen (modern)
return(list(unit=expression(S/m), scale=if (is.null(scale)) "" else scale))
if (U == "MHO/CM" || U == "MHOS/CM") # 1 mho (archaic) = 1 Siemen (modern)
return(list(unit=expression(S/cm), scale=if (is.null(scale)) "" else scale))
if (U == "MMHO") # 1 mho (archaic) = 1 Siemen (modern)
return(list(unit=expression(mS), scale=if (is.null(scale)) "" else scale))
if (U == "NBS SCALE")
return(list(unit=expression(), scale=if (is.null(scale)) "NBS" else scale))
if (U == "NTU")
return(list(unit=expression(NTU), scale=if (is.null(scale)) "" else scale))
if (U == "PPM")
return(list(unit=expression(ppm), scale=if (is.null(scale)) "" else scale))
if (U == "PSS-78" || U == "PSU")
return(list(unit=expression(), scale=if (is.null(scale)) "PSS-78" else scale))
if (U == "PMOL/KG")
return(list(unit=expression(pmol/kg), scale=if (is.null(scale)) "" else scale))
if (U == "PSU")
return(list(unit=expression(), scale=if (is.null(scale)) "PSS-78" else scale))
if (U == "ML/L")
return(list(unit=expression(ml/l), scale=if (is.null(scale)) "" else scale))
if (U == "S" || U == "SEC" || U == "SECOND")
return(list(unit=expression(s), scale=if (is.null(scale)) "" else scale))
if (u == "s/m")
return(list(unit=expression(s/m), scale=if (is.null(scale)) "" else scale))
if (u == "S/m")
return(list(unit=expression(S/m), scale=if (is.null(scale)) "" else scale))
if (u == "Total scale")
return(list(unit=expression(), scale=if (is.null(scale)) "" else scale))
if (u == "True degrees")
return(list(unit=expression(degree), scale=if (is.null(scale)) "" else scale))
if (u == "uA")
return(list(unit=expression(mu*a), scale=if (is.null(scale)) "" else scale))
# > stringi::stri_escape_unicode() indicates that Greek mu is "\u00b5"
if (u == "\u00b5einsteins/s/m^2" || U == "UEINSTEINS/S/M**2" || U == "UEINSTEINS/S/M^2")
return(list(unit=expression(mu*mol/m^2/s), scale=if (is.null(scale)) "" else scale))
# > stringi::stri_escape_unicode() indicates that Greek mu is "\u00b5"
if (u == "\u00b5M")
return(list(unit=expression(mu*M), scale=if (is.null(scale)) "" else scale))
if (U == "UMOL/KG")
return(list(unit=expression(mu*mol/kg), scale=if (is.null(scale)) "" else scale))
if (u == "ug/l")
return(list(unit=expression(mu*g/l), scale=if (is.null(scale)) "" else scale))
if (grepl("^mmol/m\\*\\*3$", unit, ignore.case=TRUE))
return(list(unit=expression(mmol/m^3), scale=if (is.null(scale)) "" else scale))
if (grepl("^mmol/m\\^3$", unit, ignore.case=TRUE))
return(list(unit=expression(mmol/m^3), scale=if (is.null(scale)) "" else scale))
if (U == "UMOL/KG")
return(list(unit=expression(mmol/kg), scale=if (is.null(scale)) "" else scale))
if (U == "UMOL/M**3")
return(list(unit=expression(mu*mol/m^3), scale=if (is.null(scale)) "" else scale))
if (U == "UMOL/M**2/S")
return(list(unit=expression(mu*mol/m^2/s), scale=if (is.null(scale)) "" else scale))
if (U == "UMOL PHOTONS/M2/S")
return(list(unit=expression(mu*mol/m^2/s), scale=if (is.null(scale)) "" else scale))
if (U == "UTC")
return(list(unit=expression(), scale=if (is.null(scale)) "" else scale))
if (U == "V")
return(list(unit=expression(V), scale=if (is.null(scale)) "" else scale))
if (U == "1/CM")
return(list(unit=expression(1/cm), scale=if (is.null(scale)) "" else scale))
if (U == "1/M")
return(list(unit=expression(1/m), scale=if (is.null(scale)) "" else scale))
if (U == "VOLT" || U == "VOLTS")
return(list(unit=expression(V), scale=if (is.null(scale)) "" else scale))
if (U == "%")
return(list(unit=expression(percent), scale=if (is.null(scale)) "" else scale))
# If none of the above worked, just try our best.
return(list(unit=as.expression(unit), scale=if (is.null(scale)) "" else scale))
}
## #' Rename a duplicated item (used in reading CTD files)
## #'
## #' Determine a new name for an item that is already in a list of names. This is
## #' done by e.g. appending a `2` to the second occurrence of a name, etc.
## #' The purpose is to create distinct variable names for
## #' [read.ctd.sbe()].
## #'
## #' @param existingNames Vector of strings with names already processed.
## #' @param name String with a candidate name.
## #' @return names String with an unduplicated name.
## #' @seealso [unduplicateNames()] is similar, but considers
## #' a vector of names.
## #'
## #' @examples
## #' unduplicateName("a", c("a", "b", "a")) # returns "a3"
## unduplicateName <- function(name, existingNames)
## {
## counter <- 0
## for (i in seq_along(existingNames)) {
## if (name == existingNames[i])
## counter <- counter + 1
## }
## res <- if (counter > 0) paste(name, counter+1, sep="") else name
## ## message("unduplicateName() name: '", name, "'")
## ## message(" existingNames '", paste(existingNames, collapse="' '"))
## ## message(" returning '", res, "'")
## res
## }
#' Rename duplicated character strings
#'
#' Append numeric suffices to character strings, to avoid repeats.
#' This is used by various data
#' input functions, to handle the fact that several oceanographic data
#' formats permit the reuse of variable names within a given file.
#'
#' @param strings Vector of character strings.
#'
#' @param style An integer giving the style. If `style`
#' is `1`, then e.g. a triplicate of `"a"` yields
#' `"a"`, `"a1"`, and `"a2"`.
#' If `style` is `2`, then the same input yields
#' `"a_001"`, `"a_002"`, and `"a_003"`.
#'
#' @return Vector of strings with repeats distinguished by suffix.
#'
#' @seealso Used by [read.ctd.sbe()] with `style=1` to
#' rename repeated data elements (e.g. for multiple temperature sensors)
#' in CTD data, and by [read.odf()] with `style=2` on
#' key-value pairs within ODF metadata.
#'
#' @examples
#' unduplicateNames(c("a", "b", "a", "c", "b"))
#' unduplicateNames(c("a", "b", "a", "c", "b"), style=2)
unduplicateNames <- function(strings, style=1)
{
## Handle duplicated names
if (style == 1) {
for (i in seq_along(strings)) {
w <- which(strings == strings[i])
lw <- length(w)
if (lw > 1) {
w <- w[-1]
strings[w] <- paste(strings[i], 1+seq.int(1, length(w)), sep="")
}
}
} else if (style == 2) {
for (i in seq_along(strings)) {
w <- which(strings == strings[i])
lw <- length(w)
if (lw > 1) {
suffices <- seq_len(lw)
strings[w] <- sprintf("%s_%03d", strings[i], suffices)
}
}
} else {
stop("unknown style=", style, "; it must be 1 or 2")
}
strings
}
#' Rename items in the data slot of an oce object (**deprecated**)
#'
#' This was deprecated in December 2019, because [oceRenameData()] does
#' a better job and is more consistent with other functions that work
#' with items in the `data` and `metadata` slots.
## This function may be used to rename elements within the
## `data` slot of `oce` objects. It also updates
## the processing log of the returned object, indicating
## the changes.
#'
#' @param x an [oce-class] object.
#'
#' @param old Vector of strings, containing old names.
#'
#' @param new Vector of strings, containing old names.
#'
## @examples
## data(ctd)
## new <- renameData(ctd, "temperature", "temperature68")
## new <- oceSetData(new, name="temperature",
## value=T90fromT68(new[["temperature68"]]),
## unit=list(unit=expression(degree*C),scale="ITS=90"))
renameData <- function(x, old=NULL, new=NULL)
{
.Deprecated("oceRenameData", msg="Superceded by oceRenameData(), as of December 2019")
if (is.null(old)) stop("need to supply old")
if (is.null(new)) stop("need to supply new")
n <- length(old)
if (n != length(new)) stop("lengths of old and new must match")
Old <- names(x@data)
New <- Old
for (i in 1:n) {
w <- which(Old == old[i])
if (length(w) == 0) stop("'", old[i], "' is not in the data slot of x")
if (length(w) > 1) stop("multiple matches are not permitted")
## message("i: ", i, ", old[i]: ", old[i], ", w:", w)
New[w] <- new[i]
}
## ensure unique ... this is a common user error
if (length(New) != length(unique(New))) stop("cannot have two columns of same name")
names(x@data) <- New
x@processingLog <- processingLogAppend(x@processingLog, paste(deparse(match.call()), sep="", collapse=""))
x
}
#' Calculate a rounded bound, rounded up to mantissa 1, 2, or 5
#'
#' @param x a single positive number
#'
#' @return for positive x, a value exceeding x that has mantissa 1, 2, or 5; otherwise, x
bound125 <- function(x)
{
x <- x[1] # ignore all but first element
if (x <= 0) {
res <- x
} else {
exp10 <- 10^floor(log10(x))
xx <- x / exp10
m <- if (xx <= 1) 1 else if (xx <=2) 2 else if (xx <= 5) 5 else 10
res <- m * exp10
##> r <- 10^rnorm(1e4)
##> R <- unlist(lapply(1:1e4, function(i) bound125(r[i]))
##> range(r/R)
##> message("x: ", x, ", exp10: ", exp10, ", m: ", m, ", xx: ", xx, ", res: ", res)
}
res
}
#' Put longitude in the range from -180 to 180
#'
#' @param longitude in degrees East, possibly exceeding 180
#'
#' @return longitude in signed degrees East
#'
#' @seealso
#' [matrixShiftLongitude()] and [shiftLongitude()] are more
#' powerful relatives to `standardizeLongitude`.
standardizeLongitude <- function(longitude) ifelse(longitude > 180, longitude-360, longitude)
#' Try to associate data names with units, for use by summary()
#'
#' Note that the whole object is not being given as an argument;
#' possibly this will reduce copying and thus storage impact.
#'
#' @param names the names of data within an object
#'
#' @param units the units from metadata
#'
#' @return a vector of strings, with blank entries for data with unknown units
#'
#' @examples
#' library(oce)
#' data(ctd)
#' dataLabel(names(ctd@@data), ctd@@metadata$units)
dataLabel <- function(names, units)
{
res <- names
## message("in dataLabel()")
if (!missing(units)) {
## message(" dataLabel(); next line is names")
## print(names)
## message(" dataLabel(); next line is units")
## print(units)
unitsNames <- names(units)
##message(" dataLabel(); next line is unitsNames")
##print(unitsNames)
for (i in seq_along(names)) {
##message(" i: ", i, ", name: ", names[i])
w <- which(unitsNames == names[i])
if (length(w)) {
##message(" we match a unit at index w=", paste(w, collapse=" "))
u <- units[w]
if (!is.null(u)) {
if (is.character(u)) {
res[i] <- paste(res[i], " [", u, "]", sep="")
} else if (is.list(u)) {
res[i] <- paste(res[i], " [", u$unit[[1]], u$scale, "]", sep="")
}
}
}
}
}
##> message("names:", paste(names, collapse=" | "))
##> message("units:", paste(units, collapse=" | "))
##> message("res:", paste(res, collapse=" | "))
res <- gsub(" *\\[\\]", "", res)
##message("dataLabel() returning:")
##print(res)
res
}
#' Capitalize first letter of each of a vector of words
#'
#' This is used in making labels for data names in some ctd functions
#'
#' @param w vector of character strings
#'
#' @return vector of strings patterned on `w` but with first letter
#' in upper case and others in lower case
titleCase <- function(w)
{
unlist(lapply(seq_along(w),
function(i) paste(toupper(substr(w[i], 1, 1)),
tolower(substr(w[i], 2, nchar(w[i]))), sep="")))
}
#' Curl of 2D vector field
#'
#' Calculate the z component of the curl of an x-y vector field.
#'
#' The computed component of the curl is defined by \eqn{\partial }{dv/dx -
#' du/dy}\eqn{ v/\partial x - \partial u/\partial y}{dv/dx - du/dy} and the
#' estimate is made using first-difference approximations to the derivatives.
#' Two methods are provided, selected by the value of `method`.
#'
#' * For `method=1`, a centred-difference, 5-point stencil is used in
#' the interior of the domain. For example, \eqn{\partial v/\partial x}{dv/dx}
#' is given by the ratio of \eqn{v_{i+1,j}-v_{i-1,j}}{v[i+1,j]-v[i-1,j]} to the
#' x extent of the grid cell at index \eqn{j}{j}. (The cell extents depend on
#' the value of `geographical`.) Then, the edges are filled in with
#' nearest-neighbour values. Finally, the corners are filled in with the
#' adjacent value along a diagonal. If `geographical=TRUE`, then `x`
#' and `y` are taken to be longitude and latitude in degrees, and the
#' earth shape is approximated as a sphere with radius 6371km. The resultant
#' `x` and `y` are identical to the provided values, and the
#' resultant `curl` is a matrix with dimension identical to that of
#' `u`.
#'
#' * For `method=2`, each interior cell in the grid is considered
#' individually, with derivatives calculated at the cell center. For example,
#' \eqn{\partial v/\partial x}{dv/dx} is given by the ratio of
#' \eqn{0.5*(v_{i+1,j}+v_{i+1,j+1}) -
#' 0.5*(v_{i,j}+v_{i,j+1})}{0.5*(v[i+1,j]+v[i+1,j+1]) - 0.5*(v[i,j]+v[i,j+1])}
#' to the average of the x extent of the grid cell at indices \eqn{j}{j} and
#' \eqn{j+1}{j+1}. (The cell extents depend on the value of
#' `geographical`.) The returned `x` and `y` values are the
#' mid-points of the supplied values. Thus, the returned `x` and `y`
#' are shorter than the supplied values by 1 item, and the returned `curl`
#' matrix dimensions are similarly reduced compared with the dimensions of
#' `u` and `v`.
#'
#' @param u matrix containing the 'x' component of a vector field
#'
#' @param v matrix containing the 'y' component of a vector field
#'
#' @param x the x values for the matrices, a vector of length equal to the
#' number of rows in `u` and `v`.
#'
#' @param y the y values for the matrices, a vector of length equal to the
#' number of cols in `u` and `v`.
#'
#' @param geographical logical value indicating whether `x` and `y`
#' are longitude and latitude, in which case spherical trigonometry is used.
#'
#' @param method A number indicating the method to be used to calculate the
#' first-difference approximations to the derivatives. See \dQuote{Details}.
#'
#' @return A list containing vectors `x` and `y`, along with matrix
#' `curl`. See \dQuote{Details} for the lengths and dimensions, for
#' various values of `method`.
#'
#' @section Development status.: This function is under active development as
#' of December 2014 and is unlikely to be stabilized until February 2015.
#'
#' @author Dan Kelley and Chantelle Layton
#'
#' @examples
#' library(oce)
#' ## 1. Shear flow with uniform curl.
#' x <- 1:4
#' y <- 1:10
#' u <- outer(x, y, function(x, y) y/2)
#' v <- outer(x, y, function(x, y) -x/2)
#' C <- curl(u, v, x, y, FALSE)
#'
#' ## 2. Rankine vortex: constant curl inside circle, zero outside
#' rankine <- function(x, y)
#' {
#' r <- sqrt(x^2 + y^2)
#' theta <- atan2(y, x)
#' speed <- ifelse(r < 1, 0.5*r, 0.5/r)
#' list(u=-speed*sin(theta), v=speed*cos(theta))
#' }
#' x <- seq(-2, 2, length.out=100)
#' y <- seq(-2, 2, length.out=50)
#' u <- outer(x, y, function(x, y) rankine(x, y)$u)
#' v <- outer(x, y, function(x, y) rankine(x, y)$v)
#' C <- curl(u, v, x, y, FALSE)
#' ## plot results
#' par(mfrow=c(2, 2))
#' imagep(x, y, u, zlab="u", asp=1)
#' imagep(x, y, v, zlab="v", asp=1)
#' imagep(x, y, C$curl, zlab="curl", asp=1)
#' hist(C$curl, breaks=100)
#' @family things relating to vector calculus
curl <- function(u, v, x, y, geographical=FALSE, method=1)
{
if (missing(u)) stop("must supply u")
if (missing(v)) stop("must supply v")
if (missing(x)) stop("must supply x")
if (missing(y)) stop("must supply y")
if (length(x) <= 1) stop("length(x) must exceed 1 but it is ", length(x))
if (length(y) <= 1) stop("length(y) must exceed 1 but it is ", length(y))
if (length(x) != nrow(u)) stop("length(x) must equal nrow(u)")
if (length(y) != ncol(u)) stop("length(x) must equal ncol(u)")
if (nrow(u) != nrow(v)) stop("nrow(u) and nrow(v) must match")
if (ncol(u) != ncol(v)) stop("ncol(u) and ncol(v) must match")
if (!is.logical(geographical)) stop("geographical must be a logical quantity")
method <- as.integer(round(method))
if (1 == method)
res <- do_curl1(u, v, x, y, geographical)
else if (2 == method)
res <- do_curl2(u, v, x, y, geographical)
else
stop("method must be 1 or 2")
res
}
#' Calculate Range, Extended a Little, as is Done for Axes
#'
#' This is analogous to what is done as part of the R axis range calculation,
#' in the case where `xaxs="r"`.
#'
#' @param x a numeric vector.
#'
#' @param extend fraction to extend on either end
#'
#' @return A two-element vector with the extended range of `x`.
#'
#' @author Dan Kelley
rangeExtended <- function(x, extend=0.04) # extend by 4% on each end, like axes
{
if (length(x) == 1) {
x * c(1 - extend, 1 + extend)
} else {
r <- range(x, na.rm=TRUE)
d <- diff(r)
c(r[1] - d * extend, r[2] + d * extend)
}
}
#' Apply a function to vector data
#'
#' The function `FUN` is applied to `f` in bins specified by
#' `xbreaks`. (If `FUN` is [mean()],
#' consider using [binMean2D()] instead, since it should be faster.)
#'
#' @param x a vector of numerical values.
#'
#' @param f a vector of data to which the elements of `FUN` may be
#' supplied
#'
#' @param xbreaks values of x at the boundaries between bins; calculated using
#' [pretty()] if not supplied.
#'
#' @param FUN function to apply to the data
#'
#' @param \dots arguments to pass to the function `FUN`
#'
#' @return A list with the following elements: the breaks in x and y
#' (`xbreaks` and `ybreaks`), the break mid-points (`xmids` and
#' `ymids`), and a matrix containing the result of applying function
#' `FUN` to `f` subsetted by these breaks.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' ## salinity profile with median and quartile 1 and 3
#' data(ctd)
#' p <- ctd[["pressure"]]
#' S <- ctd[["salinity"]]
#' q1 <- binApply1D(p, S, pretty(p, 30), function(x) quantile(x, 1/4))
#' q3 <- binApply1D(p, S, pretty(p, 30), function(x) quantile(x, 3/4))
#' plotProfile(ctd, "salinity", col='gray', type='n')
#' polygon(c(q1$result, rev(q3$result)),
#' c(q1$xmids, rev(q1$xmids)), col='gray')
#' points(S, p, pch=20)
#'
#' @family bin-related functions
binApply1D <- function(x, f, xbreaks, FUN, ...)
{
if (missing(x)) stop("must supply 'x'")
if (missing(f)) stop("must supply 'f'")
if (missing(xbreaks)) xbreaks <- pretty(x, 20)
if (missing(FUN)) stop("must supply 'FUN'")
if (!is.function(FUN)) stop("'FUN' must be a function")
## FIXME: maybe employ the code below to get data from oce objects
##if ("data" %in% slotNames(x)) # oce objects have this
## x <- x@data
##t <- try(x <- data.frame(x), silent=TRUE)
##if (class(t) == "try-error")
## stop("cannot coerce 'data' into a data.frame")
fSplit <- split(f, cut(x, xbreaks, include.lowest=TRUE, labels=FALSE))
##message("length(xbreaks)=", length(xbreaks))
##message("length(fSplit)=", length(fSplit))
result <- unlist(lapply(fSplit, FUN, ...))
result[!is.finite(result)] <- NA
names(result) <- NULL
## Put some NAs at start and end of 'result', if required because of
## 'xbreaks' bins that have no 'x' data.
xmin <- min(x, na.rm=TRUE)
xmax <- max(x, na.rm=TRUE)
nxbreaks <- length(xbreaks)
if (xmin > xbreaks[2]) {
m <- which(xbreaks < xmin)[1]
##message("condition 1; m=", m)
result <- c(rep(NA, m), result)
}
if (xmax < xbreaks[nxbreaks]) {
m <- which(xbreaks > xmax)[1]
##message("condition 2; m=", m)
result <- c(result, rep(NA, nxbreaks-m))
}
list(xbreaks=xbreaks, xmids=xbreaks[-1]-0.5*diff(xbreaks), result=result)
}
#' Apply a function to matrix data
#'
#' The function `FUN` is applied to `f` in bins specified by