applyMagneticDeclination

[EOGrady21]

Short summary of problem

Attempting to use the applyMagneticDeclination() function to correct magnetic variation in velocity components on an ADCP timeseries.

Details (optional)

applyMagneticDeclination() seems to reject all objects (tried array, matrix and list of velocity timeseries) except current meters, may be helpful to expand to apply to adp-class objects / arrays or clarify documentation.

What you did

####magnetic variation correction####
#using average over start and end points of timeseries

#start and end times
startTime <- adp[['time']][1]
endTime <- max(adp[['time']])

#start and end magnetic declinations
mStart <- magneticField(lon, lat, startTime)
mEnd <- magneticField(lon, lat, endTime)
mAvg <- round(mean(c(mStart$declination, mEnd$declination)), digits = 2)

#apply magnetic declination to velocity components
applyMagneticDeclination(adp[['v']], declination = mAvg, debug = getOption("oceDebug"))

What happened

Error in applyMagneticDeclination(adp[['v']], declination = mAvg, debug = getOption("oceDebug")) : 
  cannot apply declination to object of class array

How urgent is this?

Not urgent

Output from sessionInfo()

R version 3.4.4 (2018-03-15)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] oce_0.9-23     gsw_1.0-5      testthat_2.0.0

loaded via a namespace (and not attached):
[1] compiler_3.4.4 magrittr_1.5   R6_2.2.2       tools_3.4.4    rlang_0.2.0 

[dankelley]

Thanks. CR and I discussed this today, and I think we agree that we ought to have another function that does arbitrary rotations of (somewhat) arbitrary objects.

The reason you're getting the error is that applyMagneticDeclination only works for objects of class cm (for current meters), and you are supplying it with an array.

Thanks!!

PS. I'm going to edit your comment, putting in what are called "code fences" so that the code segments will be colourized.

[dankelley]

If @Echisholm21 and @richardsc want to suggest a function name, that would help. I suppose the function should be able to take more than one angle (e.g. a person might want to rotate about a vertical axis and then a horizontal axis) but that might be confusing to the user so maybe it would be best to have a function to rotate about the z axis, another about the x axis, etc. Hm, C, we ought to have written something on paper when we were chatting...

[dankelley]

Below is some R code that shows how to tackle this problem. The first part is a trial function that may go into oce (hence the "roxygen" comments ... these would get turned into documentation when "oce" is built) and the second part is a test (in fact, it's the test that will be accomplished by using example(rotateAboutZ) if this makes it into "oce".

I'd be interested in comments from @Echisholm21 and the oce developers on this potential function (and, particularly, its name and the phrasing of the documentation I've sketched).

##
## 1: code that may go into oce
##

#' Rotate velocity components within an oce object
#' @param x An oce object of class \code{adp}, \code{adv} or \code{cm}.
#' @param angle The rotation angle about the z axis, in degrees counterclockwise.
#' @author Dan Kelley
#' @examples
#' library(oce)
#' par(mfcol=c(2, 3))
#' # adp (acoustic Doppler profiler)
#' data(adp)
#' plot(adp, which="uv")
#' mtext("adp", side=3, line=0, adj=1, cex=0.7)
#' adpRotated <- rotateAboutZ(adp, 30)
#' plot(adpRotated, which="uv")
#' mtext("adp rotated 30 deg", side=3, line=0, adj=1, cex=0.7)
#' # adv (acoustic Doppler velocimeter)
#' data(adv)
#' plot(adv, which="uv")
#' mtext("adv", side=3, line=0, adj=1, cex=0.7)
#' advRotated <- rotateAboutZ(adv, 125)
#' plot(advRotated, which="uv")
#' mtext("adv rotated 125 deg", side=3, line=0, adj=1, cex=0.7)
#' # cm (current meter)
#' data(cm)
#' plot(cm, which="uv")
#' mtext("cm", side=3, line=0, adj=1, cex=0.7)
#' cmRotated <- rotateAboutZ(cm, 30)
#' plot(cmRotated, which="uv")
#' mtext("cm rotated 30 deg", side=3, line=0, adj=1, cex=0.7)
#'
#' @family things related to \code{adp} data
#' @family things related to \code{adv} data
#' @family things related to \code{cm} data
rotateAboutZ <- function(x, angle)
{
    S <- sin(angle * pi / 180)
    C <- cos(angle * pi / 180)
    rotation <- matrix(c(C, S, -S, C), nrow=2)
    res <- x
    if (inherits(x, "adp")) {
        if (x[["oceCoordinate"]] != "enu")
            stop("cannot rotate adp unless coordinate system is 'enu'; see ?toEnu or ?xyzToEnu")
        V <- x[["v"]]
        ## Work through the bins, transforming a 3D array operation to a
        ## sequence of 2D matrix operations.
        for (j in seq_len(dim(V)[2])) {
            uvr <- rotation %*% t(V[, j, 1:2])
            V[, j, 1] <- uvr[1, ]
            V[, j, 2] <- uvr[2, ]
        }
        res@data$v <- V
    } else if (inherits(x, "adv")) {
        if (x[["oceCoordinate"]] != "enu")
            stop("cannot rotate adv unless coordinate system is 'enu'; see ?toEnu or ?xyzToEnu")
        V <- x[["v"]]
        uvr <- rotation %*% t(V[, 1:2])
        V[, 1] <- uvr[1, ]
        V[, 2] <- uvr[2, ]
        res@data$v <- V
    } else if (inherits(x, "cm")) {
        uvr <- rotation %*% rbind(x@data$u, x@data$v)
        res@data$u <- uvr[1, ]
        res@data$v <- uvr[2, ]
    } else {
        stop("cannot handle objects of type '", class(x), "'")
    }
    ## Update processing log
    res@processingLog <- processingLogAppend(res@processingLog,
                                             paste("rotateAboutZ(x, angle=", angle, ")", sep=""))
    res
}

##
## 2: test (bottom panels should have data aligned vertically)
##
library(oce)
par(mfcol=c(2, 3))
#
data(adp)
plot(adp, which="uv")
mtext("adp", side=3, line=0, adj=1, cex=0.7)
adpRotated <- rotateAboutZ(adp, 30)
plot(adpRotated, which="uv")
mtext("adp rotated 30 deg", side=3, line=0, adj=1, cex=0.7)
#
data(adv)
plot(adv, which="uv")
mtext("adv", side=3, line=0, adj=1, cex=0.7)
advRotated <- rotateAboutZ(adv, 125)
plot(advRotated, which="uv")
mtext("adv rotated 125 deg", side=3, line=0, adj=1, cex=0.7)
#
data(cm)
plot(cm, which="uv")
mtext("cm", side=3, line=0, adj=1, cex=0.7)
cmRotated <- rotateAboutZ(cm, 30)
plot(cmRotated, which="uv")
mtext("cm rotated 30 deg", side=3, line=0, adj=1, cex=0.7)

[EOGrady21]

Thanks Dan!

The name makes sense and I appreciate the more general nature of the function. Very helpful!

[dankelley]

Done in "develop" branch, commit aeb9163d1ffc0fd7f4a29b2879156d1dccd013e6

[dankelley]

Dear reporter,

If you think the issue (as defined by its title) has been addressed, please close it. This comment is added to an issue when it seems that an issue has been addressed and that discussion has ceased. Please bear in mind that we work on issues as defined by their titles. Narrow-scoped titles help us to focus our development efforts, and benefit other oce users.

If you think the issue has not been addressed, please look through all the comments and write a new one after this, stating what still needs to be done.

If you see that the issue has meandered from the topic defined by its title, please consider closing this and opening a new issue.

Thanks!

PS. This is a standardized reply.