ISO week

I am working with a model that produces estimates of snow water equivalent through time. Because I deal with large spatial extents, I decided to have the model produce weekly averages. The problem with this is knowing which file to access for a given date. The snow files are saved using an ISO-8061 week number (I had no idea how many date "standards" there were, but this is a common one). I thought that getting the week number out of a date format in R would be easy, but I was wrong. It turns out to be platform specific and regardless, as far as I could tell, there is no built-in way to calculate the ISO week. Thus the following function.

ISOweek<-function(date,format="%Y-%m-%d",return.val="weekofyear"){
  ##converts dates into "dayofyear" or "weekofyear", the latter providing the ISO-8601 week
  ##date should be a vector of class Date or a vector of formatted character strings
  ##format refers to the date form used if a vector of
  ##  character strings  is supplied
  
  ##convert date to POSIXt format 
  if(class(date)[1]%in%c("Date","character")){
    date=as.POSIXlt(date,format=format)
  }
  if(class(date)[1]!="POSIXt"){
    print("Date is of wrong format.")
    break
  }else if(class(date)[2]=="POSIXct"){
  date=as.POSIXlt(date)
}


  if(return.val=="dayofyear"){
    ##add 1 because POSIXt is base zero
    return(date$yday+1)
  }else if(return.val=="weekofyear"){
    ##Based on the ISO8601 weekdate system,
    ## Monday is the first day of the week
    ## W01 is the week with 4 Jan in it.
    year=1900+date$year
    jan4=strptime(paste(year,1,4,sep="-"),format="%Y-%m-%d")
    wday=jan4$wday
    
    wday[wday==0]=7  ##convert to base 1, where Monday == 1, Sunday==7
    
    ##calculate the date of the first week of the year
    weekstart=jan4-(wday-1)*86400  
    weeknum=ceiling(as.numeric((difftime(date,weekstart,units="days")+0.1)/7))

    #########################################################################
    ##calculate week for days of the year occuring in the next year's week 1.
    #########################################################################
    mday=date$mday
    wday=date$wday
    wday[wday==0]=7
    year=ifelse(weeknum==53 & mday-wday>=28,year+1,year)
    weeknum=ifelse(weeknum==53 & mday-wday>=28,1,weeknum)
  
    ################################################################
    ##calculate week for days of the year occuring prior to week 1.
    ################################################################
    
    ##first calculate the numbe of weeks in the previous year
    year.shift=year-1
    jan4.shift=strptime(paste(year.shift,1,4,sep="-"),format="%Y-%m-%d")
    wday=jan4.shift$wday
    wday[wday==0]=7  ##convert to base 1, where Monday == 1, Sunday==7
    weekstart=jan4.shift-(wday-1)*86400
    weeknum.shift=ceiling(as.numeric((difftime(date,weekstart)+0.1)/7))

    ##update year and week
    year=ifelse(weeknum==0,year.shift,year)
    weeknum=ifelse(weeknum==0,weeknum.shift,weeknum)
    
    return(list("year"=year,"weeknum"=weeknum))
  }else{
    print("Unknown return.val")
    break
  }
}

Of course after I wrote this function, I found this thread in R-help. Gustaf Rydevik provided a function that is similar to mine (his function is also about twice as fast); however, on my computer it was giving incorrect years and week numbers for days in January that occur in the previous year's final week (e.g., 1 January 2010 should be week 53 of 2009). Reading through the rest of the thread indicates that the confusion involving week numbers is widespread. I guess the moral of the story is pick a standard and stick to it.

Please let me know in the comments if you encounter any problems running this function (or if you can suggest ways to make it more efficient).

Leap years

A quick function that when provided a numeric vector of years returns a boolean vector where TRUE == Leap year.

is.leapyear=function(year){
  #http://en.wikipedia.org/wiki/Leap_year
  return(((year %% 4 == 0) & (year %% 100 != 0)) | (year %% 400 == 0))
}

Truncated Normal Distribution

Many distributions may be used to describe patterns that are non-negative; however, there are not as many choices when an upper bound is also needed (although the beta distribution is very flexible). For various reasons, truncated distributions are sometimes preferred, and the truncated normal is particularly popular. While R has a package that includes the standard functions for this distribution (see rtnorm, dtnorm, etc. in the msm pacakge), the true expectation and variance of the distribution may be of interest. It turns out that the first two moments of the truncated normal are not too hard to calculate (but worth writing functions for):

mean.tnorm<-function(mu,sd,lower,upper){
  ##return the expectation of a truncated normal distribution
  lower.std=(lower-mu)/sd
  upper.std=(upper-mu)/sd
  mean=mu+sd*(dnorm(lower.std)-dnorm(upper.std))/
      (pnorm(upper.std)-pnorm(lower.std))
  return(mean)
}

var.tnorm<-function(mu,sd,lower,upper){
  ##return the variance of a truncated normal distribution
  lower.std=(lower-mu)/sd
  upper.std=(upper-mu)/sd
  variance=sd^2*(1+(lower.std*dnorm(lower.std)-upper.std*dnorm(upper.std))/
    (pnorm(upper.std)-pnorm(lower.std))-((dnorm(lower.std)-dnorm(upper.std))/
    (pnorm(upper.std)-pnorm(lower.std)))^2)
  return(variance)
}

###Testing
> library(msm)
> a=rtnorm(1000000,-5,2,1,3)
> paste(mean(a),var(a))
[1] "1.52135857341077 0.197281057170982"
> paste(mean.tnorm(-5,2,1,3),var.tnorm(-5,2,1,3))
[1] "1.52090857118 0.197111175109889"

Plotting contours

Plenty of packages allow you to plot contours of a "z" value; however, I wanted to be able to plot a specific density contour of a sample from a bivariate distribution over a plot that was a function of the x and y parameters. The example only plots the density of x and y; however, if you set plot.dens=FALSE, the contours will be added to the active display device (i.e., over an image(x,y,function(x,y){...})).

This function also lets you specify the line widths and types for each of the contours.

###########################################
## drawcontour.R
## Written by J.D. Forester, 17 March 2008
###########################################

##This function draws an approximate density contour based on
##empirical, bivariate data.


##change testit to FALSE if sourcing the file
testit=TRUE


draw.contour<-function(a,alpha=0.95,plot.dens=FALSE, line.width=2, line.type=1, limits=NULL, density.res=300,spline.smooth=-1,...){
  ##a is a list or matrix of x and y coordinates (e.g., a=list("x"=rnorm(100),"y"=rnorm(100)))
  ## if a is a list or dataframe, the components must be labeled "x" and "y"
  ## if a is a matrix, the first column is assumed to be x, the second y
  ##alpha is the contour level desired
  ##if plot.dens==TRUE, then the joint density of x and y are plotted,
  ##   otherwise the contour is added to the current plot.
  ##density.res controls the resolution of the density plot

  ##A key assumption of this function is that very little probability mass lies outside the limits of
  ## the x and y values in "a". This is likely reasonable if the number of observations in a is large.
  
  require(MASS)
  require(ks)
  if(length(line.width)!=length(alpha)){
    line.width <- rep(line.width[1],length(alpha))
  }

  if(length(line.type)!=length(alpha)){
    line.type <- rep(line.type[1],length(alpha))
  }
  
  if(is.matrix(a)){
    a=list("x"=a[,1],"y"=a[,2])
  }
  ##generate approximate density values
  if(is.null(limits)){
    limits=c(range(a$x),range(a$y))
  }
  f1<-kde2d(a$x,a$y,n=density.res,lims=limits)

  ##plot empirical density
  if(plot.dens) image(f1,...)
  
  if(is.null(dev.list())){
    ##ensure that there is a window in which to draw the contour
    plot(a,type="n",xlim=limits[1:2],ylim=limits[3:4],...)
  }
    
  ##estimate critical contour value
  ## assume that density outside of plot is very small

  zdens <- rev(sort(f1$z))
  Czdens <- cumsum(zdens)
  Czdens <- (Czdens/Czdens[length(zdens)])
  for(cont.level in 1:length(alpha)){
    ##This loop allows for multiple contour levels
    crit.val <- zdens[max(which(Czdens<=alpha[cont.level]))]
   
    ##determine coordinates of critical contour
    b.full=contourLines(f1,levels=crit.val)
    for(c in 1:length(b.full)){
      ##This loop is used in case the density is multimodal or if the desired contour
      ##  extends outside the plotting region
      b=list("x"=as.vector(unlist(b.full[[c]][2])),"y"=as.vector(unlist(b.full[[c]][3])))
      
      ##plot desired contour
      line.dat<-xspline(b,shape=spline.smooth,open=TRUE,draw=FALSE)
      lines(line.dat,lty=line.type[cont.level],lwd=line.width[cont.level])
    }
  }
}

##############################
##Test the function
##############################

##generate data
if(testit){
  n=10000
  a<-list("x"=rnorm(n,400,100),"y"=rweibull(n,2,100))
  draw.contour(a=a,alpha=c(0.95,0.5,0.05),line.width=c(2,1,2),line.type=c(1,2,1),plot.dens=TRUE, xlab="X", ylab="Y")
}

Approximate sunrise and sunset times

This function is not perfect, but it does a reasonable job estimating sunrise and sunset times for my field site. If more accurate data are required, try here. Note, the command to calculate day of year is: strptime(x, "%m/%d/%Y")$yday+1

suncalc<-function(d,Lat=48.1442,Long=-122.7551){
  ## d is the day of year
  ## Lat is latitude in decimal degrees
  ## Long is longitude in decimal degrees (negative == West)
  
  ##This method is copied from:
  ##Teets, D.A. 2003. Predicting sunrise and sunset times.
  ##  The College Mathematics Journal 34(4):317-321.
 
  ## At the default location the estimates of sunrise and sunset are within
  ## seven minutes of the correct times (http://aa.usno.navy.mil/data/docs/RS_OneYear.php)
  ## with a mean of 2.4 minutes error.

  ## Function to convert degrees to radians
  rad<-function(x)pi*x/180
  
  ##Radius of the earth (km)
  R=6378
  
  ##Radians between the xy-plane and the ecliptic plane
  epsilon=rad(23.45)

  ##Convert observer's latitude to radians
  L=rad(Lat)

  ## Calculate offset of sunrise based on longitude (min)
  ## If Long is negative, then the mod represents degrees West of
  ## a standard time meridian, so timing of sunrise and sunset should
  ## be made later.
  timezone = -4*(abs(Long)%%15)*sign(Long)

  ## The earth's mean distance from the sun (km)
  r = 149598000

  theta = 2*pi/365.25*(d-80)

  z.s = r*sin(theta)*sin(epsilon)
  r.p = sqrt(r^2-z.s^2)

  t0 = 1440/(2*pi)*acos((R-z.s*sin(L))/(r.p*cos(L)))
  
  ##a kludge adjustment for the radius of the sun
  that = t0+5 

  ## Adjust "noon" for the fact that the earth's orbit is not circular:
  n = 720-10*sin(4*pi*(d-80)/365.25)+8*sin(2*pi*d/365.25)

  ## now sunrise and sunset are:
  sunrise = (n-that+timezone)/60
  sunset = (n+that+timezone)/60

  return(list("sunrise" = sunrise,"sunset" = sunset))
}

Convert polar coordinates to Cartesian

When I want to calculate the coordinates of a location (e.g., a nest or burrow) based on distance and bearing from a grid point, this function helps me avoid writing down SOH-CAH-TOA every time. Just note that the bearing in this case is from the grid point (known location) to the unknown location.

polar2cart<-function(x,y,dist,bearing,as.deg=FALSE){
  ## Translate Polar coordinates into Cartesian coordinates
  ## based on starting location, distance, and bearing
  ## as.deg indicates if the bearing is in degrees (T) or radians (F)
  
  if(as.deg){
    ##if bearing is in degrees, convert to radians
    bearing=bearing*pi/180
  }
  
  newx<-x+dist*sin(bearing)  ##X
  newy<-y+dist*cos(bearing)  ##Y
  return(list("x"=newx,"y"=newy))
}


##Example

oldloc=c(0,5) 
bearing=200 #degrees
dist=5

newloc<-polar2cart(oldloc[1],oldloc[2],dist,bearing,TRUE)
plot(oldloc[1],oldloc[2],xlim=c(-10,10),ylim=c(-10,10))
points(newloc$x,newloc$y,col="red")

Variance-Covariance Matrix in glm

This is a small function Venables and Ripley provide in their MASS book. You don't need it anymore because vcov() has a method for the glm class. However, it is useful to see how to extract bits from a fitted model object.

vcov.glm<-function(obj){
  #return the variance-covariance matrix of a glm object
  #from p. 188 in Venables and Ripley. 2002. 
  #Modern Applied Statistics With S. Springer. New York.
  so <- summary(obj,corr=F)
  so$dispersion * so$cov.unscaled
}

Calculate turning angles and step lengths from location data

anglefun <- function(xx,yy,bearing=TRUE,as.deg=FALSE){
  ## calculates the compass bearing of the line between two points
  ## xx and yy are the differences in x and y coordinates between two points
  ## Options:
  ## bearing = FALSE returns +/- pi instead of 0:2*pi
  ## as.deg = TRUE returns degrees instead of radians
  c = 1
  if (as.deg){
    c = 180/pi
  }
  
  b<-sign(xx)
  b[b==0]<-1  #corrects for the fact that sign(0) == 0
  tempangle = b*(yy<0)*pi+atan(xx/yy)
  if(bearing){
    #return a compass bearing 0 to 2pi
    #if bearing==FALSE then a heading (+/- pi) is returned
    tempangle[tempangle<0]<-tempangle[tempangle<0]+2*pi
  }
  return(tempangle*c)
}

bearing.ta <- function(loc1,loc2,loc3,as.deg=FALSE){
  ## calculates the bearing and length of the two lines
  ##    formed by three points
  ## the turning angle from the first bearing to the
  ##    second bearing is also calculated
  ## locations are assumed to be in (X,Y) format.
  ## Options:
  ## as.deg = TRUE returns degrees instead of radians
  if (length(loc1) != 2 | length(loc2) != 2 | length(loc3) !=2){
    print("Locations must consist of either three vectors, length == 2,
or three two-column dataframes")
    return(NaN)
  }
  c = 1
  if (as.deg){
    c = 180/pi
  }
  
  locdiff1<-loc2-loc1
  locdiff2<-loc3-loc2
  bearing1<-anglefun(locdiff1[1],locdiff1[2],bearing=F)
  bearing2<-anglefun(locdiff2[1],locdiff2[2],bearing=F)

  if(is.data.frame(locdiff1)){
    dist1<-sqrt(rowSums(locdiff1^2))
    dist2<-sqrt(rowSums(locdiff2^2))
  }else{
    dist1<-sqrt(sum(locdiff1^2))
    dist2<-sqrt(sum(locdiff2^2))
  }
  
  ta=(bearing2-bearing1)
  
  ta[ta < -pi] = ta[ta < -pi] + 2*pi
  ta[ta > pi] = ta[ta > pi] - 2*pi
  return(list(bearing1=unlist(bearing1*c),bearing2=unlist(bearing2*c),
  ta=unlist(ta*c),dist1=unlist(dist1),dist2=unlist(dist2)))
}