R – Tips and Tricks

Collection of R function I tend to forget about or I am frequnetly asked by people to help them out with…
Functions:

  1. Round value to nearest x
  2. Calculate difference between the n and n-1 element of a variable
  3. Calculate time difference between two strings
  4. Generate Date-Time string from character
  5. Trigonometric functions in degrees x
  6. Check if variable exists
  7. Check if folder exists, and create if False
  8. Parallelise Processes
  9. Load functions from external file
  10. List all files in folder recursively
  11. Read a NetCDF (*.nc) file
  12. Publish directly to WordPress
  13. Select numeric, character or factor columns only from dataframe
  14. Check package installation and install if missing
  15. Read table with missing values
  16. Replace NA with last available value
    Fisheries Acoustic Related
  17. Nearfield of an acoustic transducer at given soundspeed and frequency
  18. Beam overlap of standard Simrad EK60/ES70 Transducers
  19. Speed of Sound under water, Coppens & MacKenzie
  20. Height of deadzone


Round to nearest x

#Round to nearest x
roundx <- function(value,x){
    round(value/x,0)*x
}


Calculate difference between the n and n-1 element of a variable

    variable<-sample(1:1000,100,replace=T)
    diff(variable)


Calculate time difference between two strings

    timestr <- c("10:53:15","10:56:47")
    #transform into numeric seconds
    seconds <- as.numeric(as.difftime(timestr,units="secs"))
    #calculate difference
    diff(seconds)


Generate Time-Date string from character

   dat <- c("2012-07-13", "2012-06-13","2013-07-13","2015-07-14")
   time <- c("15:01:32", "15:02:42", "18:26:31", "18:27:19" )
   strptime( paste(dat,time), "%Y-%m-%d %H:%M:%S")


Trigonometric functions – degrees to radians

    radians <- degrees*pi/180
    cos(radians)
    sin(radians)


Check if variable exists

    ifelse(exists("ff"), paste(x, "- exists already!"), "Variable not existing")


Folder check

    #Check if folder exists, otherwise create it
    dir.create(file.path(getwd(), "check/this/folder"), showWarnings=F)

Parallelise Processes

    library(doParallel)
    #Define cluser with N cores
    N <- 10
    cl <- makeCluster(N)
    registerDoParallel(cl)

    #Run foreach in parallel
    foreach(N) %dopar%{
         myfunctiondoessomestuff(N)
    }

Load functions from external workspace

source("myfunctionfile.R"]

List all files in a folder recursively, given a pattern

allRaw<-list.files(datawd, pattern='\\.raw$', recursive=T, full.names=T)


Read a netCDF (*.nc) file

library(ncdf4)
ncdat <- nc_open(".\\myNetCDFfile.nc")
var <- "selected/dataname"
#Get the entire variable out
cloud <- ncvar_get( ncdat, var)
#Get a subset of the variable, in this case the 10th flag
cloud2 <- ncvar_get( ncdat, var,c(10,1),verbose=TRUE) #Make it a mask cloud2[cloud2&amp;amp;amp;gt;0]<-1
#Simple plot
image(cloud2)
#Ggplot verion
library(ggplot2)
ggplot(melt(cloud2), aes(Var1,Var2, fill=value)) + geom_raster()


Publish to wordpress from within R

#Install the RWordPress library if needed
if (!require('RWordPress'))
  install.packages('RWordPress', repos = 'http://www.omegahat.org/R', type = 'source')
#load the libraries
library(knitr)
library(RWordPress)
#set login options
options(WordPressLogin = c(user = 'password'),
        WordPressURL = 'https://user.wordpress.com/xmlrpc.php')
#knit the R Markdown file to wordpress
knit2wp('MyMarkdownFile.Rmd',action='newPage', title = 'Bottom Detection', shortcode=TRUE)

user and password have to be change in WordPressLogin and WordPressURL
For the knit2wp function shortcode=TRUE is recommended for WordPress.com users only, otherwise a ncier code highlighting is recommended



Select columns based on class

#create a data frame
data <- as.data.frame(cbind(Numbers=sample(1:100,10),
                             Fact=sample(c("Bla","blub"),10,replace=T),
                             Char=sample(c("character","string","text"),10,replace=T)))
data[,1]<-as.numeric(data[,1])
data[,2]<-as.character(data[,2])#for numeric values
#numeric
data[sapply(data,is.numeric)]
#character
data[sapply(data, is.character)]
#factor
data[sapply(data, is.factor)]


Check for package and install if missing

#list of packages to be installed
list.of.packages <- c("lunar")
#check if pakage if installed
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
#if a package is missing
if(length(new.packages)) install.packages(new.packages)


Read table with missing data

#read table with space as separator 
#read.csv for csv files or sep="\t" for tab separated data
#skip for skipping first lines
#fill to fill missing values with empty cells
#na.strings - replace empty by NA
read.table("myfile.txt",sep=" ",skip=3,fill=TRUE, na.strings=c("","NA"))


Replace NA by last available value

#zoo library is required
# na.locf where locf = last one carried forward
data<- cbind(A=c(1,sample(c(1:5,NA),10,replace=T)),
             B=c(2,sample(c(1:5,NA),10,replace=T)),
             C=c(3,sample(c(1:5,NA),10,replace=T)))
transform(data, A = na.locf(A), B = na.locf(B), C = na.locf(C))

 

Nearfield of an acoustic transducer

#near field
#R=2*L^2*l^-1
#L=transducer diameter
#l=wavelength
Rnf <- function(c,f,beamangle){
  #first wavelength
  lambda <- c/(f*1000)
  #wave number
  k=2*pi/lambda
  #active radius
  a=1.6*100/(k*sin((beamangle/2)*pi/180))
  #near-field
  (2*a/100)^2/lambda
}
#################################Example################################
c=1500
f=120
beamangle=7
Rnf(c,f,beamangle)
f=c(18,38,50,70,120,333)
b=c(11,7,7,7,7,7)
plot(f,Rnf(c,f,b),type="ol",xlab="Frequency [kHz]",ylab="Near-field [m]")

Beam overlap

#Compute overlay of typical SIMRAD EK60 Transducers, 
#f1 = First Frequencyy
#f2 = Second Frequency
#d = distance between the transducers [m]

over<-function(f1,f2,d){
if(f1==18){
  bwr0=11
  
}else{
  bwr0=7
}
if(f2==18){
  bwr1=11
}else{
  bwr1=7
}

range <-c(1,seq(2,3.8,by=.1),seq(3.81,3.9,by=.01),seq(4,5,by=.1),seq(6,17,by=1),seq(17,18,by=.1),seq(19,234))
r0 <- tan((pi/180)*bwr0/2)*range
r1 <- tan((pi/180)*bwr1/2)*range

CBA=range
for(i in 1:length(range)){
if(r1[i]+r0[i]-d<0){
  CBA[i]=0
}else{
  if(d+r1[i]-r0[i]<0.1){
    CBA[i]=1
  }else{
    CBA[i]=acos(((r1[i]^2+d^2-r0[i]^2)/(2*r1[i]*d)))
  }
}
}

CBD=2*CBA

CAB=acos((r0^2+d^2-r1^2)/(2*r0*d))
CAD=2*CAB

a_over<-(CBD*r1^2/2)-((sin(CBD)*r1^2)/2)+(CAD*r0^2/2)-((sin(CAD)*r0^2)/2)
a_TXDCR <- pi*r0^2

perc_o<-range
for (i in 1:length(range)){
if(CBA[i]==0){
  perc_o[i]<-0
}else{
  if(CBA[i]==1){
    perc_o[i]<-(pi*r1^2)/(pi*r0^2)*100
  }else{
    perc_o[i]<-a_over[i]/a_TXDCR[i]*100
  }
}

}
as.data.frame(cbind(range=range,
                    r0=r0,
                    r1=r1,
                    CBA=CBA,
                    CBD=CBD,
                    CAB=CAB,
                    CAD=CAD,
                    AreaOver=a_over,
                    AreaTXDCR=a_TXDCR,
                    perOver=perc_o))
}

############################### Example ##############################
ov <- over(18,38,0.6053509838)
ov2 <- over(38,120,0.358050276)
ov3 <- over(38,200,0.45607017)
ov4 <- over(120,200,0.570087713)
plot(data=ov,perOver~range,type="l",ylim=c(0,100),col="green",ylab="Percent Overlay",xlab="Range [m]")
lines(data=ov2,perOver~range,col="blue")
lines(data=ov3,perOver~range,col="red")
lines(data=ov4,perOver~range,col="orange")
legend(190, 25,c("18 - 38 kHz","38 - 120 kHz","38 - 200 kHz","120 - 200 kHz"), 
       lty=c(1,1,1,1),
       col=c("green","blue","red","orange"),cex=0.7,lwd=2,bty="n")

Speed of sound underwater, Coppens & MacKenzie

####Range of validity: temperature 2 to 30 Celsius, salinity 25 to 40 parts per thousand, depth 0 to 8000 m
#  * K.V. Mackenzie, Nine-term equation for the sound speed in the oceans (1981) J. Acoust. Soc. Am. 70(3), pp 807-812 *
mac <- function(S,T,D) {
  c_mackenzie &amp;amp;amp;amp;lt;-  1448.96 +4.591*T-5.304*10^(-2)*T^2+2.374*10^(-4)*T^3+1.340*(S-35)+1.630*10^(-2)*D+1.675*10^(-7)*D^2-1.025*10^(-2)*T*(S - 35)-7.139*10^(-13)*T*D^3
}
#Coppens (1981)

cop <- function(S,T,D) {
  t <- T/10 # T in 10/ degrees Celsius
  D <- D/1000 #D in km
  c0 <-1449.05+45.7*t-5.21*t^2+0.23*t^3+(1.333-0.126*t+0.009*t^2)*(S - 35) 
  
  c0+(16.23+0.253*t)*D+(0.213-0.1*t)*D^2+(0.016+0.0002*(S-35))*(S-35)*t*D
}

############################## Example ################################################################
  s <- seq(34,36,by=0.1)  
  t <- seq(10,30,by=1)
  d <- seq(0,1000,by=50)
  ca<-cop(T=t,S=s,D=d)
  ma<-mac(T=t,S=s,D=100)
  
  par(mfrow=c(1,3))
  plot(ca~t,type="l",ylab="Speed of sound",xlab="Temperature",cex.lab=2,cex.axis=1.5)
  lines(ma~t,lty=2)
  points(ca~t,col="red",pch=16,cex=1.2)
  points(ma~t,col="blue",pch=16,cex=.8)
  legend("bottomright", legend = c("Coppens", "MacKenzie"),
         pch = 16,col=c("red","blue"),cex=1.5)
  
  
  plot(ca~s,type="l",xlab="Salinity",cex.lab=2,cex.axis=1.5)
  lines(ma~s,lty=2)
  points(ca~s,col="red",pch=16,cex=1.2)
  points(ma~s,col="blue",pch=16,cex=.8)
  legend("bottomright", legend = c("Coppens", "MacKenzie"),
         pch = 16,col=c("red","blue"),cex=1.5)
  
  plot(ca~d,type="l",xlab="Depth [m]",cex.lab=2,cex.axis=1.5)
  lines(ma~d,lty=2)
  points(ca~d,col="red",pch=16,cex=1.2)
  points(ma~d,col="blue",pch=16,cex=.8)
  legend("bottomright", legend = c("Coppens", "MacKenzie"),
         pch = 16,col=c("red","blue"),cex=1.5)
  

Height of acoustic deadzone


#d = depth[m],
#theta = beam angle (typical 7 degrees,
#c=sound speed [m/s],
#tau = pulse duration [sec])
hdz <- function(d,theta,c,tau){
d*(1-cos(theta*pi/180/2))+c*tau/2
}
############ Example ################
hdz(100,7,1450,0.0004)