detections.r

#Written: R Branton, Mar 2012

# modified R Branton, Apr 9, 2012 to include modifications made at Orono workshop

# Purpose: for use at Orono Maine to demonstrate use OTN detection data with RStudio, RCurl, Lattice, RgoogleMaps etc ...
# Notes: 1st time users must edit and execute the following command line to create .OTNUserID file ...
# cat('USERNAME:PASSWORD@',file=".OTNUserID")

# Set these variables to specify desired server
my.serviceName<-'kil-otn-2w1.its.dal.ca'  # use this when running from RStudio on OTN server
my.serviceName<-'kil-otn-2w1-x.its.dal.ca' # use either this or the next when running from RStudio on desktop

# Remaining code is fixed ...

options(
useFancyQuotes = FALSE,
stringsAsFactors=FALSE,
as.is=TRUE)

library(RCurl)
library(RgoogleMaps)
library(manipulate)
library(geosphere)
library(lattice)

# define fuction to get data from wfs server

get.data<-function(typeName=''){
my.url<-paste('http://',scan(".OTNUserID","character",quiet=TRUE),my.serviceName,":8080/geoserver/wfs?request=getfeature&service=wfs&version=1.1.0&outputformat=CSV&typename=",typeName,sep="")
cat(getURL(my.url),file="temp.csv")
read.csv(file="temp.csv",sep=",",header=TRUE)}

a<-rbind(get.data('pbnUnit:mv_anm_released'),
get.data('dasUnit:mv_anm_released'),
get.data('pbuUnit:mv_anm_released'))

table(a$collectioncode,a$scien)

table(a$sci,a$collectioncode)
a<-a[a$sci=='Salmo salar',]
table(a$sci,a$life)
a<-a[a$life!='ADULT',]

a2<-aggregate(rep(1,length(a$jul)),by=list(year=a$yearc,day=a$jul,coll=a$collectioncode),FUN=sum)

xyplot(x~day|factor(year)*coll,data=a2,type='h',main="numbers released by jday, year and collection")

x<-rbind(get.data('pbnUnit:mv_anm_compressed'),get.data('pbuUnit:mv_anm_compressed'))

# create line segment series by merging detection point series onto itself

rm(x2)

x2<-merge(
x=data.frame(collectioncode=x$col,catalognumber=x$cat,year=substr(x$dates,1,4),scientificname=x$sci,station.start=x$sta,seq_num=x$seq,total_count=x$tot,datestart=x$dates,longitude.start=x$lon,latitude.start=x$lat),
y=data.frame(catalognumber=x$cat,station.end=x$sta,seq_num=x$seq-1,total_count=x$tot,dateend=x$dates,longitude.end=x$lon,latitude.end=x$lat),
by=c("catalognumber","seq_num")
)

# map line segments as arrows

MyMap<-GetMap.bbox(range(c(x2$longitude.s,x2$longitude.e)),range(c(x2$latitude.s,x2$latitude.e)),destfile="temp.png",maptype='hybrid')
PlotArrowsOnStaticMap(MyMap,lon0=x2$longitude.s,lat0=x2$latitude.s,lon1=x2$longitude.e,lat1=x2$latitude.e,FUN=arrows,col="green",verbose=0,length=.1,lwd=1)

# merge animal length to line segment series and examine using hist

x2<-merge(
x=x2,
y=data.frame(catalognumber=a$cat,length=a$length),
by="catalognumber"
)

hist(a$length)
table(a$sci,a$collectioncode)

# calculate duration, distance and speed (blps) for each line segment

x2$duration<-as.numeric(strptime(x2$datee,'%Y-%m-%dT%H:%M:%S')-strptime(x2$dates,'%Y-%m-%dT%H:%M:%S'))
x2$distance<-distCosine(matrix(c(x2$longitude.start,x2$latitude.start),ncol=2),matrix(c(x2$longitude.end,x2$latitude.end),ncol=2))
x2$speed<-(x2$dist/x2$length)/x2$dur

print(histogram(~speed|collectioncode,data=x2,layout=c(1,2)))

print(densityplot(~length|collectioncode,data=x2,group=year,layout=c(1,2),
auto.key=list(columns=length(table(x2$year))),
main="size composition"))

# select subset of release ending on halifax line

x3<-x2[x2$seq_num==1&substr(x2$station.end,1,3)=='HFX',]

# compare size at release with size at detection by year

x4=rbind(
data.frame(g='released',l=a$length,y=a$yearcollected,c=a$collectioncode)
,data.frame(g='detected',l=x3$length,y=x3$y,c=x3$coll))

print(densityplot(~l|y*c,group=g,data=x4,subset=(x4$y>=2008),xlab="length",auto.key=list(column=2)))

# xyplot days v. km by year grouped by collectioncode

print(xyplot(duration/(60*60*24)~distance/1000|year,group=collectioncode,data=x3,auto.key=TRUE,
xlab="kilometres",ylab="days",main="Halifax Line Transits"))

# examine speed

hist(x3$speed)

print(bwplot(speed~year|collectioncode,data=x3,main="body lengths / second by year and release "))

# calculate julian dates and examine using bwplot

x3$js<-as.numeric(format(strptime(x3$dates,'%Y-%m-%d'),format='%j'))
x3$je<-as.numeric(format(strptime(x3$datee,'%Y-%m-%d'),format='%j'))
x3$jm<-(x3$je+x3$js)/2
print(bwplot(jm~year|collectioncode,data=x3,main="median date of transit"))

# use manipulate to map by by catalog number

cat(paste(
'manipulate(my.plot(catalogNumber)',',\n',
'catalogNumber=picker("all",',paste(dQuote(names(table(x2$catalognumber))),collapse=","),')','\n',
')',
collapse=''),file='temp.r')

my.plot<-function(catalogNumber="all"){
if (catalogNumber!='all'){
x3<-x2[x2$catalognumber==catalogNumber,]
} else {
x3<-x2
}
MyMap<-GetMap.bbox(range(c(x3$longitude.s,x3$longitude.e)),range(c(x3$latitude.s,x3$latitude.e)),destfile="temp.png",maptype='hybrid')
PlotArrowsOnStaticMap(MyMap,lon0=x3$longitude.s,lat0=x3$latitude.s,lon1=x3$longitude.e,lat1=x3$latitude.e,FUN=arrows,col="green",verbose=0,length=.1,lwd=1)
legend("top",legend=paste("catalognumber =",catalogNumber))
}
source('temp.r')