#######################################################################################################################################################
############# ICES TRAWL SURVEY AND EVALUATION course #################################################################################################
################# Practical 01 - Exploring the Chrons ###########################################################################################################
##################Doug Beare and Dave Reid ##########################################################################################################
## The aim of this exercise is download ICES exchange format data from www.ices.dk, extract and visualise the chron records.
#######################################################################################################################################################
## 1. Install R (http://cran.r-project.org/), packages and perhaps Tinn-R (http://sourceforge.net/projects/tinn-r/).
#######################################################################################################################################################
install.packages('maps'); # Eg. install the packages. Remember you only need to do this once.
install.packages('maptools');
install.packages('rgdal');
install.packages('sp');
install.packages('spatial');
install.packages('PBSmapping');
library(sp);library(maptools);library(rgdal);library(maps);library(PBSmapping) # attach to search path
#######################################################################################################################################################
## 2. Select some data from the ICES website - http://datras.ices.dk/Data_products/Download/Download_Data_public.aspx
#######################################################################################################################################################
#######################################################################################################################################################
## 3. Save the data to a directory of your choice, eg. 'c:/datras/data' and extract them from the zip files.
#######################################################################################################################################################
#######################################################################################################################################################
## 4. Install the R library datras. Do this by downloading datras*zip from the sharepoint. Open R, go to packages, install package from local zip files.
#######################################################################################################################################################
#######################################################################################################################################################
## 5. Attach datras package written by me to the R search path
#######################################################################################################################################################
library(datrasR)
#######################################################################################################################################################
## 6. Extract the chron records in the relevant directory. In the first instance chose a survey area in which you are interested. The suggestions
#######################################################################################################################################################
## below are simply options. Remember that you will need to adjust the working directory (wd) path so that R can find the data. The ICES website
## only allows relatively small subsets of data to be downloaded at once. For the time-being stick to fairly small subsets of data.
#chrons <- parseExchangeFormatChrons(wd="D:/bearedo/Consultancy/TrawlSurveyCourseCopenhagen2011/Practicals/RockallData")
chrons <- parseExchangeFormatChrons(wd="D:/bearedo/Consultancy/TrawlSurveyCourseCopenhagen2011/Practicals/BTS-TRIDENS-1987-2010")
#chrons <- parseExchangeFormatChrons(wd="D:/bearedo/Consultancy/TrawlSurveyCourseCopenhagen2011/Practicals/NS-IBTS-2010-2011")
#chrons <- parseExchangeFormatChrons(wd="D:/bearedo/Consultancy/TrawlSurveyCourseCopenhagen2011/Practicals/NS-IBTS-2005-2011")
#chrons <- parseExchangeFormatChrons(wd="D:/bearedo/Consultancy/TrawlSurveyCourseCopenhagen2011/Practicals/NS-IBTS-2001-2011-Q1")
#######################################################################################################################################################
## 7. Here only examine valid hauls, non-duplicates and hauls where we have both shoot and haul positions.
#######################################################################################################################################################
chrons <- chrons[chrons$haulval=='V',] # only valid hauls
chrons <- chrons[!duplicated(chrons),] # chuck any duplicates
chrons <- chrons[!is.na(chrons$haullong),] # don't want data without haul longs
chrons <- chrons[chrons$country != 'NOR',] # throw out Norway because shoot and haul positions are the same.
#######################################################################################################################################################
## 8. Print and examine first 5 lines of chron file
#######################################################################################################################################################
head (chrons)
#######################################################################################################################################################
## 9. Explore the data with R function table
#######################################################################################################################################################
table(chrons$haulval) # number of valid hauls (check)
table(chrons$hauldur) # number of hauls at different durations
table(chrons$gear) # number of gears
table(chrons$quarter,chrons$year) # number of hauls by year and quarter
table(chrons$country,chrons$ship) # number of hauls by country and ship
table(chrons$daynight,chrons$depth) # number of hauls by day/night and depth
table(chrons$ship,cut(chrons$depth,breaks=c(seq(0,250,by=25)))) # number of hauls by ship and depth bin.
table(chrons$statrec,as.character(chrons$ship)) # number of hauls by statistical rectangle and ship
table(chrons$year) # number of hauls by year
#######################################################################################################################################################
## 10. Some simply summary statistics
#######################################################################################################################################################
summary(chrons$depth)
summary(chrons$hauldur)
summary(chrons$datim.shot)
#######################################################################################################################################################
## 11. average temperature (if available) in each statrectangle.
#######################################################################################################################################################
tapply(chrons$surtemp,list(chrons$quarter,chrons$statrec),mean,na.rm=T)
#######################################################################################################################################################
## 12. average depth (if available) in each statrectangle.
#######################################################################################################################################################
tapply(chrons$depth,list(chrons$quarter,chrons$statrec),mean,na.rm=T)
#######################################################################################################################################################
## 13. average number, duration and distance of hauls in ICES area
#######################################################################################################################################################
chrons$ices_area <- ICESarea(chrons,roman=F) # Add on ICES area using function datras
tapply(chrons$hauldur,chrons$ices_area,length) # number of hauls in each ICES area
round(tapply(chrons$hauldur,chrons$ices_area,sum,na.rm=T)/60) # total haul duration in each ICES area
round(tapply(chrons$distance,chrons$ices_area,sum,na.rm=T)/1000) # total distance hauled in each ICES area
#######################################################################################################################################################
## 14. You might want to write data out into a csv file for reading with excel or some other software
#######################################################################################################################################################
write.table(chrons,file='chrons.csv',sep=',',row.names=F)
#######################################################################################################################################################
## 15. Explore the data with graphs
#######################################################################################################################################################
# Simple map of shoot and haul locations
par(mfrow=c(1,1))
plot(chrons$shootlong,chrons$shootlat,pch=".",col='blue',ylab='',xlab='',cex=2)
# add points
points(chrons$haullong,chrons$haullat,pch=".",col='red',cex=2)
# add lines for stat rectangles
abline(h=0:80,v=seq(-30,30,by=.5),lty=2,col='grey',lwd=.5)
# add map
library(maps)
map('worldHires',add=T,col='darkseagreen',fill=T)
# Add depth contours (different files for different areas since the data sets are so big).
data(nseaBathy) # north sea bathymetry
data(nwestBathy) # north west atlantic shelf bathymetry
data(balticBathy) # baltic sea bathymetry
contour(nseaBathy$x,nseaBathy$y,nseaBathy$z,levels=c(20,40,60,100),add=T) # add on depth contours for north sea
contour(nwestBathy$x,nwestBathy$y,nwestBathy$z,levels=c(100,200,500),add=T) # add on depth contours for baltic
#######################################################################################################################################################
## 16. Explore data in various ways looking for potential confounding ###
#######################################################################################################################################################
what.quarter <- 3 # remember to add on the correct quarter you are examining
chrons1 <- chrons[chrons$quarter == what.quarter,]
par(mfrow=c(1,2),mar=c(4,4,2,2))
plot(chrons1$shootlong,chrons1$shootlat,pch=16,col='blue',ylab='',xlab='') # just shoot long
# add lines for stat rectangles
abline(h=0:80,v=seq(-30,30,by=.5),lty=2,col='grey')
# add map
library(maps)
map('world',add=T,col='green',fill=T) # low res map. Try 'worldHires' too.
tt <- tapply(chrons1$depth,format(chrons1$datim.shot,"%H"),median,na.rm=T) # calculate mean depths at different times of day
plot(as.numeric(names(tt)),tt,type='l',ylab="average depth",xlab="time of day")
## Plot the survey by year
chrons1 <- chrons[chrons$quarter == what.quarter,] # extract data according to what.quarter
uy <- sort(unique(chrons$year))
print(uy)
print(length(uy)) # how many unique years in the chron file ??
par(mfrow=c(3,5),mar=c(2,2,2,2)) # use the result of print(length(uy)) to set up plot for number of graphs. Remember mfrow=c(3,4) means 'plot 12 graphs'
for (i in uy){
plot(chrons1$shootlong[chrons1$year == i],chrons1$shootlat[chrons1$year==i],pch=16,col='blue',ylab='',xlab='')
# add points
#points(jitter(chrons1$haullong),jitter(chrons1$haullat),pch=16,col='red')
# add lines for stat rectangles
abline(h=0:80,v=seq(-30,30,by=.5),lty=2,col='grey')
# add map
library(maps)
map('world',add=T,col='green',fill=T) # low res map. Try 'worldHires' too.
title(as.character(i))
}
# Comment on the distribution of the survey stations
#######################################################################################################################################################
## 17. Find chrons within a polygon. Note: The 'polygon' can also come from a shapefile.
#######################################################################################################################################################
# First create a data frame with the Round 3 Zone 5 Boundary Co-ordinates (a proposed area for an offshore wind farm).
dat.r3 <- data.frame (
lon.d=c(3,3,2,2,2,1,1,2,2,2),
lon.m=c(4,10,27,6,8,52,55,13,14,2),
lon.s=c(15.035,46.476,11.669,56.469,30.217,6.891,56.591,32.66,29.436,40.5),
lat.d=c(53,52,51,52,52,52,52,52,52,52),
lat.m=c(14,51,59,0,17,15,24,26,36,49),
lat.s=c(46.991,46.076,49.19,33.467,25.312,57.710,7.55,12.405,1.973,33.309))
# Convert to decimal lat and long co-ordinates
x <- ((dat.r3$lon.m*60)+dat.r3$lon.s)/3600
dat.r3$lon <- dat.r3$lon.d+x
y <- ((dat.r3$lat.m*60)+dat.r3$lat.s)/3600
dat.r3$lat <- dat.r3$lat.d+y
# Create polyset data.frame (from PBSmapping but there are many other alternatives).
round.3 <- data.frame(PID=rep(1,11), POS=1:11,X=c(dat.r3$lon,dat.r3$lon[1]),Y=c(dat.r3$lat,dat.r3$lat[1]))
round.3 <- as.PolySet(round.3,projection='UTM')
head(round.3) # examine data (typical method for storage of spatial data).
data(europa) # map of europe (with a detailed Dutch coastline) included in the package datras
par(mfrow=c(1,1))
plotMap(europa,xlim=c(-1,5),ylim=c(50,54),col='darkgreen',xlab='',ylab='')
addLines(round.3,lwd=2)
addPolys(round.3,col='red')
points(chrons$shootlong,chrons$shootlat,pch=".",cex=2)
# Identify the points within the polygon
events <- data.frame(EID=1:length(chrons[,1]),X=chrons$shootlong,Y=chrons$shootlat) # Must first creat an events data frame
events <- as.EventData(events,projection="UTM")
r3<-findPolys(events,round.3,maxRows=1e+06)
print(r3) #examine the output
points(chrons$shootlong[r3$EID],chrons$shootlat[r3$EID],col='blue',pch='*') # plot points INSIDE the polygon
points(chrons$shootlong[-r3$EID],chrons$shootlat[-r3$EID],col='black',pch='+') # plot points OUTSIDE the polygon
#######################################################################################################################################################
## 18. Create and plot points and tow tracks using R library sp. Output the result as kml (Google Earth) and shapefiles.
#######################################################################################################################################################
spDF <- Chrons2SpatialPointsDF(data=chrons) # Convert the data to a spatial points data.frame
# Plot map in R
data(europa) # in built map
plotMap(europa,xlim=c(-15,10),ylim=c(47,62),xlab='',ylab='') # Note: change xlim and ylim to focus in on correct area.
points(spDF,pch=".",col='red',cex=2)
# The magic of Google Earth. Note you can change the driver for other formats.
writeOGR(obj=spDF, dsn="chrons.kml", layer="shoot-locations", driver="KML")
writeOGR(obj=spDF, dsn="chrons.shp", layer="shoot-locations", driver="ESRI Shapefile") # A shapefile.
print(getwd()) # see which directory the file has gone too, find it and open it in Google Earth
# Similarly tow tracks instead of points
uships <- unique(chrons$ship) # how many unique ships in the dataset ?
print(uships)
spLDF <- ChronsTrack2SpatialLineDF(input=chrons,what.ship = c("TRI2")) # convert to a spatial line data frame for a particular ship. Note the function takes too long if you select too many ships/data.
par(mfrow=c(1,1))
plotMap(europa,xlim=c(-5,10),ylim=c(47,62),xlab='',ylab='',col='darkseagreen') # Plot to see
lines(spLDF,col='blue',lwd=1)
writeOGR(spLDF, "towtracks.kml", "tracks", "KML")
print(getwd()) # see which directory the file has gone too, find it and open file towtracks.kml in Google Earth
####################### FURTHER EXERCISES ####################################
# Explore chron files from different datasets downloaded from ICES website.
# Make your own polygons for selecting data.
# Find examples of 'confounding'/non-random sampling in survey data and share them with the group at Plenary, eg. create maps with dots color-coded for different ships, gears, etc.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.