###script R, version 23/05/2016
###authors: Christine Plumejeaud-Perreau and Anne S. Philippe

###load thworking directory and the benthos dataset
setwd("C:/Users/Desktop/XXXXXXX")
benthos<- read.csv("benthos.csv", sep=";", dec=",", h=T)


							#####################



#####Script 1: Creates a bubble plot of Scrobicularia biomass densities in Aiguillon Bay in 2014(site AIC)#####
	#load packages
library(sp)
library(gstat)
library(rgdal)
library(raster)
library(reshape)

	#Create a dataframe with Scrobicularia biomass for each station at site AIC in 2014
scrAIC<-subset(benthos, benthos$st_sit_id=="aic"&benthos$year==2014&benthos$abr=="scr")
subs<-scrAIC[,c("abr","AFDM","Poskey","year")]
temp<-melt(subs, id=c("abr","Poskey","year"))
tempaic<-cast(temp, year+Poskey~abr,sum) ####now a data frame with summed densities per Poskey per Year per Group
	#Create a data frame "areaaic" with corresponding sampling areas and calculates biomass densities at each station
areaaic<-scrAIC[,c(1,21)]
areaaic<-areaaic[!duplicated(areaaic), ]#only keep the eara value corresponding to aech station
	#Merge the two dataframe and calculate biomass densities at each station
aicscr<-merge(tempaic, areaaic, by="Poskey", all=TRUE)###
aicscr$scr.m²<-aicscr$scr/aicscr$area
	#In order to plot these densities spatially, add coordinates to the dataframe "aicscr"
coordinates<-subset(benthos, benthos$st_sit_id=="aic")
coordinates<-coordinates[,c(1,25,26)]
coordinates<-coordinates[!duplicated(coordinates), ] ##creates a dataframe with coordinates for each station in site AIC
aicscr<-merge(aicscr, coordinates, by="Poskey", all=T)
	#Plot biomass densities of Scrobicularia plana in 2014 in site AIC using bubble function
aicscr$scr.m²[is.na(aicscr$scr.m²)] <- 0 ##Nas are true zeros
coordinates(aicscr) <- ~ st_long_dd + st_lat_dd
proj4string(aicscr) <- CRS("+init=epsg:4326") # def L93
bubble(aicscr, "scr.m²", main="AIC 2014 scr_AFDM", maxsize=3, minsize=0,col="black", key.entries=c(0,0.1,5,10,20,40))

rm(scrAIC,subs,temp,areaaic,coordinates)#remove object that is not necessary anymore



							#####################É



#####Script 2: Creates a plot of Scrobicularia biomass densities at station 1698 in Aiguillon Bay between 2004 and 2014(site AIC)#####
scr1698<-subset(benthos, benthos$abr=="scr"&benthos$st_sit_id=="aic"&benthos$Poskey==1698)
library(reshape)
subs<-scr1698[,c("abr","biomass_dens","Poskey","year")]
temp<-melt(subs, id=c("abr","Poskey","year"))
temp<-cast(temp,,sum)
scr1698<-as.data.frame(temp)
plot(scr1698$year,scr1698$biomass_dens, type="o", pch=20, ylab="Biomass",xlab="year", bty="n")

rm(subs, temp)#remove objects that are not necessary anymore


								#####################É



#####Script 3: Creates pie plots representing the share of principal mollusc genus in the total biomass density for each site in 2014#####

	#Create a dataframe "pies" with only the principal Genus : Abra, Scrobicularia, Macoma, Peringia and Cerastoderma
v<-c("abr","abo","abt","abn")#in the dataset, 4 abreviations refer to the Genus Abra (see Table 1)
benthos$Genus[benthos$abr%in%v]<-"Abra"
benthos$Genus[benthos$abr=="scr"]<-"Scrobicularia"
benthos$Genus[benthos$abr=="mac"]<-"Macoma"
benthos$Genus[benthos$abr=="hyd"]<-"Peringia"
benthos$Genus[benthos$abr=="cer"]<-"Cerestoderma"
pies<-subset(benthos, !is.na(benthos$Genus)) ##only keep the rows where Genus is either Abra, Scrobicularia, Macoma, Peringia or Cerastoderma (62777 rows)

	#work site by site to create dataframes for the future pie plot with total biomass of principal species per site and station, here we start with site AIC, then AIV, MO and OL
piesAIC<-subset(pies, pies$st_sit_id=="aic"&pies$year==2014) #"piesAIC" is a dataframe with only rows from the delected Genus, site AIC and year 2014 (884 rows)
library(reshape)
subs<-piesAIC[,c("Genus","AFDM","Poskey","year")]
temp<-melt(subs, id=c("Genus","Poskey","year"))
piesAIC<-cast(temp, year+Poskey~Genus,sum) #now a data frame with summed densities per Poskey per Year per Group
piesAIC[is.na(piesAIC)] <- 0 #all NAs are true zeros
rm(subs, temp) # remove the uncessary objects from the working environment
	#same thing for AIV
piesAIV<-subset(pies, pies$st_sit_id=="aiv"&pies$year==2014)  #"piesAIV" is a dataframe with only rows from the delected Genus, site AIV and year 2014 (1256 rows)
subs<-piesAIV[,c("Genus","AFDM","Poskey","year")]
temp<-melt(subs, id=c("Genus","Poskey","year"))
piesAIV<-cast(temp, year+Poskey~Genus,sum) ####now a data frame with summed densities per Poskey per Year per Group
piesAIV[is.na(piesAIV)] <- 0
rm(subs, temp) # remove the uncessary objects from the working environment
	#same thing for MO
piesMO<-subset(pies, pies$st_sit_id=="mo"&pies$year==2014) #"piesMO" is a dataframe with only rows from the delected Genus, site MO and year 2014 (986 rows)
subs<-piesMO[,c("Genus","AFDM","Poskey","year")]
temp<-melt(subs, id=c("Genus","Poskey","year"))
piesMO<-cast(temp, year+Poskey~Genus,sum) ####now a data frame with summed densities per Poskey per Year per Group
piesMO[is.na(piesMO)] <- 0
rm(subs, temp) # remove the uncessary objects from the working environment
	#same thing for OL
piesOL<-subset(pies, pies$st_sit_id=="ol"&pies$year==2014) #"piesOL" is a dataframe with only rows from the delected Genus, site OL and year 2014 (2186 rows)
subs<-piesOL[,c("Genus","AFDM","Poskey","year")]
temp<-melt(subs, id=c("Genus","Poskey","year"))
piesOL<-cast(temp, year+Poskey~Genus,sum) ####now a data frame with summed densities per Poskey per Year per Group
piesOL[is.na(piesOL)] <- 0
rm(subs, temp) # remove the uncessary objects from the working environment


	###calculate the total area sampled for all Genus (except Peringia ) for all stations, starting with site AIC 
areaaic<-subset(pies, pies$st_sit_id=="aic"&pies$year==2014)
areaaic<-areaaic[,c(1,21)] #selects the column "area"
areaaic<-areaaic[!duplicated(areaaic), ]
areaaic<-sum(areaaic$area)# sum the areas sampled at each station to have the total area sampled in site AIC (for all species ecept Peringia) [areaaic=1.00531 m²]
###calculate total area sampled for Peringia (another sampling protocol is used for Peringia, see section Sampling Method)
areaaichyd<-subset(pies, pies$st_sit_id=="aic"&pies$year==2014)
areaaichyd<-areaaichyd[,c(1,22)]
areaaichyd<-areaaichyd[!duplicated(areaaichyd), ]
areaaichyd<-sum(areaaichyd$area_hyd) # sum the areas sampled at each station to have the total area sampled in site AIC for Peringia [areaaichyd=0.5026548 m²]

piesAIC$Abra.m²<-sum(piesAIC$Abra)/areaaic
piesAIC$Cerastoderma.m²<-sum(piesAIC$Cerastoderma)/areaaic
piesAIC$Macoma.m²<-sum(piesAIC$Macoma)/areaaic
piesAIC$Peringia.m²<-sum(piesAIC$Peringia)/areaaichyd  ##!!! the sampled surface is different for Peringia
piesAIC$Scrobicularia.m²<-sum(piesAIC$Scrobicularia)/areaaic
aic<-piesAIC[1,c(1,8,9,10,11,12)]

	#same for sites AIV, MO and OL
     ###calculate total area sampled for molluscs and Peringia
areaaiv<-subset(pies, pies$st_sit_id=="aiv"&pies$year==2014)
areaaiv<-areaaiv[,c(1,21)]
areaaiv<-areaaiv[!duplicated(areaaiv), ]
areaaiv<-sum(areaaiv$area)
areaaivhyd<-subset(pies, pies$st_sit_id=="aiv"&pies$year==2014)
areaaivhyd<-areaaivhyd[,c(1,22)]
areaaivhyd<-areaaivhyd[!duplicated(areaaivhyd), ]
areaaivhyd<-sum(areaaivhyd$area_hyd)
areamo<-subset(pies, pies$st_sit_id=="mo"&pies$year==2014)
areamo<-areamo[,c(1,21)]
areamo<-areamo[!duplicated(areamo), ]
areamo<-sum(areamo$area)
areamohyd<-subset(pies, pies$st_sit_id=="mo"&pies$year==2014)
areamohyd<-areamohyd[,c(1,22)]
areamohyd<-areamohyd[!duplicated(areamohyd), ]
areamohyd<-sum(areamohyd$area_hyd)
areaol<-subset(pies, pies$st_sit_id=="ol"&pies$year==2014)
areaol<-areaol[,c(1,21)]
areaol<-areaol[!duplicated(areaol), ]
areaol<-sum(areaol$area)
areaolhyd<-subset(pies, pies$st_sit_id=="ol"&pies$year==2014)
areaolhyd<-areaolhyd[,c(1,22)]
areaolhyd<-areaolhyd[!duplicated(areaolhyd), ]
areaolhyd<-sum(areaolhyd$area_hyd)
    ###create mo, ol and aiv files with mean biomass densities in 2014 (all Poskeys together)
piesAIV$Abra.m²<-sum(piesAIV$Abra)/areaaiv
piesAIV$Cerastoderma.m²<-sum(piesAIV$Cerastoderma)/areaaiv
piesAIV$Macoma.m²<-sum(piesAIV$Macoma)/areaaiv
piesAIV$Peringia.m²<-sum(piesAIV$Peringia)/areaaivhyd## Area for Peringia
piesAIV$Scrobicularia.m²<-sum(piesAIV$Scrobicularia)/areaaiv
aiv<-piesAIV[1,c(1,8,9,10,11,12)]
piesMO$Abra.m²<-sum(piesMO$Abra)/areamo
piesMO$Cerastoderma.m²<-sum(piesMO$Cerastoderma)/areamo
piesMO$Macoma.m²<-sum(piesMO$Macoma)/areamo
piesMO$Peringia.m²<-sum(piesMO$Peringia)/areamohyd## Area for Peringia
piesMO$Scrobicularia.m²<-sum(piesMO$Scrobicularia)/areamo
mo<-piesMO[1,c(1,8,9,10,11,12)]
piesOL$Abra.m²<-sum(piesOL$Abra)/areaol
piesOL$Cerastoderma.m²<-sum(piesOL$Cerastoderma)/areaol
piesOL$Macoma.m²<-sum(piesOL$Macoma)/areaol
piesOL$Peringia.m²<-sum(piesOL$Peringia)/areaolhyd### Area for Peringia
piesOL$Scrobicularia.m²<-sum(piesOL$Scrobicularia)/areaol
ol<-piesOL[1,c(1,8,9,10,11,12)]

	#create the column "All", representing the sum of all five densities, and vectors of biomass densities for the 5 taxa at each site
aic$All<-aic$Abra.m²+aic$Cerastoderma.m²+aic$Macoma.m²+aic$Peringia.m²+aic$Scrobicularia.m²
aicpie<-c(aic$Abra.m²,aic$Cerastoderma.m²,aic$Macoma.m²,aic$Peringia.m²,aic$Scrobicularia.m²) 
aiv$All<-aiv$Abra.m²+aiv$Cerastoderma.m²+aiv$Macoma.m²+aiv$Peringia.m²+aiv$Scrobicularia.m²
aivpie<-c(aiv$Abra.m²,aiv$Macoma.m²,aiv$Cerastoderma.m²,aiv$Peringia.m²,aiv$Scrobicularia.m²) 
mo$All<-mo$Abra.m²+mo$Cerastoderma.m²+mo$Macoma.m²+mo$Peringia.m²+mo$Scrobicularia.m²
mopie<-c(mo$Abra.m²,mo$Macoma.m²,mo$Cerastoderma.m²,mo$Peringia.m²,mo$Scrobicularia.m²) 
ol$All<-ol$Abra.m²+ol$Cerastoderma.m²+ol$Macoma.m²+ol$Peringia.m²+ol$Scrobicularia.m²
olpie<-c(ol$Abra.m²,ol$Macoma.m²,ol$Cerastoderma.m²,ol$Peringia.m²,ol$Scrobicularia.m²) 


	#Creates pie plots, with diameter proportionnal to the value "All", i.e. the total biomass density of all five species considered
col<-c(gray(0),gray(1),gray(0.6),gray(0.3),gray(0.9))#sets colours representing each Genus, in grey scale
labs<-c("Abra", "Cerastoderma", "Macoma", "Peringia", "Scrobicularia")###sets labels for the legend
pie(aicpie, labels=NA,col=col, cex=0.6, init.angle=180, main="",radius=0.08*aic$All)
pie(aivpie, labels=NA,col=col, cex=0.6, init.angle=180, main="",radius=0.08*aiv$All)
pie(mopie, labels=NA,col=col, cex=0.6, init.angle=180, main="",radius=0.08*mo$All)
pie(olpie, labels=NA,col=col, cex=0.6, init.angle=180, main="",radius=0.08*ol$All)
legend(-1,-1, labs,cex=0.6,fill=col,ncol=2, bty="n")

##removes the unecessary objects for the environment
rm(aic,aiv,mo,ol,v,pies,areaol, areaaic, areaaiv, areamo,piesAIC, piesAIV,piesMO, piesOL, areaaichyd,areaaivhyd,areamohyd,areaolhyd)