sensors.r

# Written: R Branton, Jan 2012
# Revised: R Branton, Jan 26, 2012 - .OTNUserID format changed
# Purpose: Demonstrate RCurl, RgoogleMaps, Lattice and other R functions for visualizing OTN sensor data ...
# Note: 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 and collectioncode
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
my.collectionCode='wdg'

library(RCurl)
library(RgoogleMaps)
library(lattice)

# get data from geoserver using RCurl

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,as.is=TRUE)}

x<-get.data(paste(my.collectionCode,'unit:mv_anm_detections',sep=''))

# briefly inspect data

names(x)
summary(x)
x[1:3,]
x[1:100,9:10]
ccode<-names(table(x$collectioncode))

# map with RgoogleMaps

my.lon<-aggregate(x$longitude,list(x$station),FUN=mean)[,'x']
my.lat<-aggregate(x$latitude,list(x$station),FUN=mean)[,'x']
my.num<-aggregate(rep(1,dim(x)[1]),list(x$station),FUN=sum)[,'x']
MyMap<-GetMap.bbox(range(my.lon),range(my.lat),destfile="temp.png",maptype='hybrid')
PlotOnStaticMap(MyMap,my.lat,my.lon,col='red',cex=log10(my.num)+1,FUN=points,verbose=0)

# note: RStudio uses must click the 'Clear All' button in the 'Plots' window to continue ...

# summarize and subset with basic R functions

table(x$detectedby,x$sensortype,x$scientificname)
x2<-x[!is.na(x$sensorvalue),]
table(x2$sensortype)
x2<-x[!is.na(x$sensorvalue),]
boxplot(split(x2$sensorvalue,x2$sensortype))

#x<-x[x$sci=='Morone saxatilis',]

# examine datecollected for detections for selected line ...

x2<-x[x$detectedby=='GMG',]
table(x2$cat,x2$sta)
x2[1:10,c('catalognumber','station','datecollected')]

# visualize sensorvalues with Lattice histogram function

print(
histogram(~sensorvalue|sensortype*detectedby,
data=x,
subset=!is.na(sensorvalue),
#    scales='free',
main=paste('Distribution of Sensorvalues by Receiver Array for ',ccode,' Tag Releases',sep='')
)
)

# summarize all depth data into 5 metre bins by detectedby

x2<-x[x$sensortype=='depth',c('detectedby','sensorvalue')]
x2<-aggregate(rep(1,dim(x2)[1]),list(depth=floor(x2$sensorvalue),detectedby=x2$detectedby),FUN=sum)
print(tapply(x2$x,list(depth=floor(x2$depth/5)*5,detectedby=x2$detectedby),FUN=sum),na.print='.')

# summarize MPS subset into 5 metre bins by station

x2<-x[x$sensortype=='depth'&x$detectedby=='GMG',c('station','sensorvalue')]
x3<-aggregate(rep(1,dim(x2)[1]),list(depth=floor(x2$sensorvalue),station=x2$station),FUN=sum)
print(tapply(x3$x,list(depth=floor(x3$depth/5)*5,station=x3$station),FUN=sum),na.print='.')

# visualize MPS Subset with boxplot, qnorm  

boxplot(split(-1*x2$sensorvalue,x2$station),main=paste('Distribution of Sensorvalues by HFX Station for ',ccode,' Tag Releases',sep=''))
qqnorm(x2$sensorvalue); qqline(x2$sensorvalue)  

# visualize as percent of total by depth and series

x2<-x[x$sensortype=='depth',c('detectedby','sensorvalue')]
x3<-aggregate(rep(1,dim(x2)[1]),list(detectedby=x2$detectedby,depth=floor(x2$sensorvalue)),FUN=sum)
x4<-sapply(split(x3$x,x3$detectedby),cumsum)
x5<-split(x3$depth,x3$detectedby)
plot(x3$depth,x3$x,ylim=c(0,100),type='n',xlab='depth (m)',ylab='percent (%)')
for (i in 1:length(x4)) {lines(x5[[i]],(x4[[i]]/max(x4[[i]])*100),col=i+1,lwd=2)}
legend("bottomright",names(table(x2$detectedby)),lty=1,lwd=2,col=(1:length(x4))+1)
title(main=paste('Percent by Depth (m) and DetectedBy for ',names(table(x$coll)),':',names(table(x$sci)),sep=''))

# visualize sensorvalue distribtion by detectedb and station using Lattice bwplot function

x2<-x[x$sensortype=='depth',c('station','sensorvalue')]   
print(
bwplot(-1*sensorvalue~station|detectedby*scientificname,
data=x2,
subset=sensortype=='depth',
scales='free',
layout=c(1,3),
ylab='depth (m)'
)
)

print(
bwplot(sensorvalue~station|detectedby*scientificname,
data=x,
subset=sensortype=='temperature',
scales='free',
ylab='temperature (c)'
)
)