# 1.1. Connection to a PostgreSQL database (on a server)
#---------------------------------------------------
library(RPostgreSQL)

drv <- dbDriver("PostgreSQL")

con <- dbConnect(drv,host='acoustica',
                 dbname="Echobase-test",user='echobase',
                 password='echobase',port=5432)

# extract map data (longitude and latitude positions)
#------------------
sqlMapData="SELECT MISSION.NAME, VOYAGE.NAME, CELL.NAME, DATA.DATAVALUE, DATAMETADATA.NAME FROM PUBLIC.VOYAGE AS VOYAGE, PUBLIC.MISSION AS MISSION, PUBLIC.CELL AS CELL, PUBLIC.DATA AS DATA, PUBLIC.DATAMETADATA AS DATAMETADATA WHERE VOYAGE.MISSION = MISSION.TOPIAID AND CELL.VOYAGE = VOYAGE.TOPIAID AND DATA.CELL = CELL.TOPIAID AND DATA.DATAMETADATA = DATAMETADATA.TOPIAID"
MapData <- dbGetQuery(con,sqlMapData)

# exercice 1: 
# - find the number of data is Echobase and check that it is identical 
# to the number of rows in MapData
# - how many MapData per voyage?

# format map data
#------------------
head(MapData)
names(MapData)=c('mission','voyage','cell','value','data')
MapDatas=MapData[MapData$data%in%c('GridCellLongitude','GridCellLatitude'),]
head(MapDatas)
xy.df=unique(merge(MapDatas[MapDatas$data=='GridCellLongitude',c('cell','value')],
                   MapDatas[MapDatas$data=='GridCellLatitude',c('cell','value')],
                   by.x='cell',by.y='cell'))
head(xy.df)
xy.df[,2:3]=apply(xy.df[,2:3],2,as.character)
xy.df[,2:3]=apply(xy.df[,2:3],2,as.numeric)

#xy.df[,2]
names(xy.df)[2:3]=c('x','y')

# plot map cell positions
#------------------
plot(as.numeric(as.character(MapDatas[MapDatas$data=='GridCellLongitude','value'])),
     as.numeric(as.character(MapDatas[MapDatas$data=='GridCellLatitude','value'])))

# extract fish map results: 
#------------------
sqlMapResults="SELECT mission.name AS mission, voyage.name AS voyage, cell.name AS cellname, celltype.id AS celltype, result.resultvalue, datametadata.name AS datametadata, species.baracoudacode FROM public.voyage AS voyage, public.mission AS mission, public.transit AS transit, public.transect AS transect, public.dataacquisition AS dataacquisition, public.cell AS cell, public.celltype AS celltype, public.result AS result, public.datametadata AS datametadata, public.category AS category, public.speciescategory AS speciescategory, public.species AS species WHERE voyage.mission = mission.topiaid AND transit.voyage = voyage.topiaid AND transect.transit = transit.topiaid AND dataacquisition.transect = transect.topiaid AND cell.voyage = voyage.topiaid AND cell.celltype = celltype.topiaid AND result.cell = cell.topiaid AND result.datametadata = datametadata.topiaid AND result.category = category.topiaid AND category.speciescategory = speciescategory.topiaid AND speciescategory.species = species.topiaid AND celltype.id = 'Map'"

MapResults <- dbGetQuery(con,sqlMapResults)
dim(MapResults)
head(MapResults)
unique(MapResults$datametadata)

# format map results
#------------------
# merge map cell data and results
names(MapResults)=c('mission','voyage','cell','celltype','value','parameter','sp')
MapResults=merge(MapResults,xy.df,by.x='cell',by.y='cell')
# convert values to numeric
MapResults$value=as.numeric(as.character(MapResults$value))
head(MapResults)
# reshape dataframe
wMapResults=reshape(MapResults,v.names='value',idvar = c("cell","mission","voyage","celltype","sp","x","y"),
                    timevar = "parameter", direction = "wide")
head(wMapResults)
wMapResults$year=substr(wMapResults$voyage,7,11)

# 3. Plot average maps over the series
#-------------------------------
# Define output path
path.maps='~/.gvfs/donnees2 sur nantes/Campagnes/bd/Echobase/ConnectionExterne/R/'
# set path.maps to NULL if you do not want to save the plots

# species list
lsp=unique(wMapResults$sp)

for (i in 1:length(lsp)){
  # select 1 species
  wMapResultsi=wMapResults[wMapResults$sp==lsp[i],]
  
  # Compute mean and stdev maps per file over the series
  #-----------------------------
  # convert Echobase data to gridFile  format
  mat=Echobase2grid(mat=wMapResultsi)
  #Compute mean and stdev maps over the series
  mmap=aggregate(mat[,c('Zvalue','Nsample')],list(Xgd=mat$Xgd,Ygd=mat$Ygd,I=mat$I,J=mat$J),mean,na.rm=TRUE)
  names(mmap)[5:6]=c("Zvalue",'Nsample')
  sdmap=aggregate(mat$Zvalue,list(Xgd=mat$Xgd,Ygd=mat$Ygd,I=mat$I,J=mat$J),sd,na.rm=TRUE)
  names(sdmap)[5]=c("Zstdev")
  msdmap=merge(mmap,sdmap,by.x=c('Xgd','Ygd','I','J'),by.y=c('Xgd','Ygd','I','J'))
  head(msdmap)
  n1=grep(FALSE,is.letters(unique(mat$Survey)))[1]
  msdmap$Survey=substr(unique(mat$Survey),1,(n1-1))
  msdmap$Year=paste(min(mat$Year),max(mat$Year),sep='-')
  msdmap=msdmap[order(msdmap$J,msdmap$I),]
  msdmap$sp=lsp[i]
  msdmap$Zstdev=1
  head(msdmap)
  # save mean and sd file
  write.table(msdmap,paste(path.maps,unique(msdmap[,c('Survey')]),unique(msdmap[,c('Year')]),lsp[i],'-MeanSD.csv',sep=''),
              sep=';',row.names=FALSE)
  
  #Plot mean maps over the series
  pat=msdmap
  #pat=paste(path1,lfs2$pat[i],'-MeanSD.csv',sep='')
  pat2=path.maps
  namz='average'
  nams='standard deviation'
  namn='average'
  
  # define titles
  namz=paste(lsp[i],'mean biomass')
  nams=paste(lsp[i],'stedv biomass')
  namn=paste(lsp[i],'No. of samples')
  # plot maps
  mat=grid.plot(pat=msdmap,input='gridDataframe',pat2=NULL,a=0,ux11=TRUE,namz=NULL,namn=NULL,nams=NULL,i=NA,
                lon1=NULL,lon2=NULL,lat1=NULL,lat2=NULL,resolution=4,newbathy=FALSE,bathy.plot=TRUE,
                deep=-450,shallow=-50,bstep=450,bcol='grey50',drawlabels=FALSE)
  
  if (i==1){
    msdmap.db=msdmap
  }else{
    msdmap.db=rbind(msdmap.db,msdmap)
  }
}

graphics.off()