`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)) }`

## Thursday, October 18, 2007

### 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

## Sunday, October 14, 2007

### 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")

## Wednesday, October 3, 2007

### Reorder factor levels

Very often, especially when plotting data, I need to reorder the levels of a factor because the default order is alphabetical. There must be many ways of reordering the levels; however, I always forget which package to look in. A direct way of reordering, using standard syntax is as follows:

In this case, the level order could be set in the first line; however, if there are many levels (and you don't want to type out all of the levels explicitly), the above method is preferred. Again, if there is a better way to do this, please let me know in the comments.

## generate datax = factor(sample(letters[1:5],100, replace=TRUE))

print(levels(x)) ## This will show the levels of x are "Levels: a b c d e"

## To reorder the levels:

## note, if x is not a factor use levels(factor(x))

x = factor(x,levels(x)[c(4,5,1:3)])

print(levels(x)) ## Now "Levels: d e a b c"

Subscribe to:
Posts (Atom)