plotrivers2=function(data,startyear=1,bis=NA,filterday=21,filteryear=11,colorss=c("#a40606","orange","#c3f1f8","#7be3f3","#0073de","#0000de"),nbcol=100,bgcolor='#c4c4c4',phi=50,r= 50,d=0.1,expand=0.5,ltheta=90,lphi=180,shade=0.35,titlep=NA,filterw='gauss')
{

#Routine to build an animated gif picture showing a 3D representation of daily river flows at a given cross river section
#Package "animation" is needed
#Data is the time series to plot. It must start at Jan 1st of a given year, and be composed of an integer number of years
#Staryear is the first year of the observation period
#Select bis=*, where * is the position of the first leap year in the given record, if you want to eliminate the observations collected
#on Feb 29.
#Filterday is the span (in days) of the filter in the intra-annual direction
#Filteryear is the span (in years) of the filter in the interannual direction
#colorss is the combination of the colors used to draw the river flow surface.
#nbcol is the number of colors used to draw the river flow surface
#bgcolor is the background color. Color "white" can be used.
#phi, r, expand, ltheta, lphi, shade: see the help page of the function persp
#titlep is the title of the picture
#filterw defines the type of filter to be used: movavg, gausss or both.

if(filterw!='gauss' & filterw!='movavg' & filterw!="both")
{
cat('Error! Filter can be gauss or movavg or both only!',fill=T)
return()
}

prova=data
if(is.na(bis)==F)
prova=delbis(prova,bis)
nanni=length(prova)/365
if((nanni-floor(nanni))>0)
{
cat("Error: the series refers to a fractional number of years",fill=T)
return()
}
prova=array(prova,dim=c(365,nanni))

if(filterw=='gauss')
{
if(filterday/2==floor(filterday/2)) filterday=filterday+1
if(filteryear/2==floor(filteryear/2)) filteryear=filteryear+1

if(filterday>=5)
{
pr100=1.96/floor(filterday/2)
pr101=-1.96
for(iiii in 1:(floor(filterday/2)-1))
{
pr101=c(pr101,-1.96+pr100*iiii)
}
pr101=c(pr101,0)
for(iiii in 1:(floor(filterday/2)))
{
pr101=c(pr101,pr100*iiii)
}
pr102=dnorm(pr101)
pr102=pr102/sum(pr102)
}

if(filteryear>=5)
{
pr100=1.96/floor(filteryear/2)
pr121=-1.96
for(iiii in 1:(floor(filteryear/2)-1))
{
pr121=c(pr121,-1.96+pr100*iiii)
}
pr121=c(pr121,0)
for(iiii in 1:(floor(filteryear/2)))
{
pr121=c(pr121,pr100*iiii)
}
pr122=dnorm(pr121)
pr122=pr122/sum(pr122)
}

if(filterday==3)
{
pr102=dnorm(c(-1.96,0,1.96))
pr102=pr102/sum(pr102)
}

if(filteryear==3)
{
pr122=dnorm(c(-1.96,0,1.96))
pr122=pr122/sum(pr122)
}

if(filteryear>=3)
{
for(i in 1:365)
{
prova[i,]=filter(prova[i,],pr122,circular=T)
}
}

if(filterday>=3)
{
for(i in 1:nanni)
{
prova[,i]=filter(prova[,i],pr102,circular=T)
}
}
}

if(filterw=='movavg')
{
for(i in 1:365)
{
prova[i,]=filter(prova[i,],rep(1/filteryear,filteryear),circular=T)
}
for(i in 1:nanni)
{
prova[,i]=filter(prova[,i],rep(1/filterday,filterday),circular=T)
}
}

if(filterw=='both')
{
if(filterday/2==floor(filterday/2)) filterday=filterday+1
if(filteryear/2==floor(filteryear/2)) filteryear=filteryear+1

if(filterday>=5)
{
pr100=1.96/floor(filterday/2)
pr101=-1.96
for(iiii in 1:(floor(filterday/2)-1))
{
pr101=c(pr101,-1.96+pr100*iiii)
}
pr101=c(pr101,0)
for(iiii in 1:(floor(filterday/2)))
{
pr101=c(pr101,pr100*iiii)
}
pr102=dnorm(pr101)
pr102=pr102/sum(pr102)
}

if(filteryear>=5)
{
pr100=1.96/floor(filteryear/2)
pr121=-1.96
for(iiii in 1:(floor(filteryear/2)-1))
{
pr121=c(pr121,-1.96+pr100*iiii)
}
pr121=c(pr121,0)
for(iiii in 1:(floor(filteryear/2)))
{
pr121=c(pr121,pr100*iiii)
}
pr122=dnorm(pr121)
pr122=pr122/sum(pr122)
}

if(filterday==3)
{
pr102=dnorm(c(-1.96,0,1.96))
pr102=pr102/sum(pr102)
}

if(filteryear==3)
{
pr122=dnorm(c(-1.96,0,1.96))
pr122=pr122/sum(pr122)
}

if(filteryear>=3)
{
for(i in 1:365)
{
prova[i,]=filter(prova[i,],pr122,circular=T)
}
}

if(filterday>=3)
{
for(i in 1:nanni)
{
prova[,i]=filter(prova[,i],pr102,circular=T)
}
}

for(i in 1:365)
{
prova[i,]=filter(prova[i,],rep(1/filteryear,filteryear),circular=T)
}
for(i in 1:nanni)
{
prova[,i]=filter(prova[,i],rep(1/filterday,filterday),circular=T)
}
}

#COLORS
jet.colors <-colorRampPalette(colorss)
# Generate the desired number of colors from this palette
color <- jet.colors(nbcol)
# Compute the z-value at the facet centres
ncz=ncol(prova)
nrz=nrow(prova)
zfacet <- prova[-1, -1] + prova[-1, -ncz] + prova[-nrz, -1] + prova[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)

library('animation')
ani.options(interval=0.1)
saveMovie(figure.ani(prova,startyear=startyear,colorss=colorss,phi=phi,r=r,d=d,expand=expand,ltheta=ltheta,lphi=lphi,shade=shade,titlep=titlep,bgcolor=bgcolor,nanni=nanni,color=color,facetcol=facetcol),outdir=getwd())
listfile=0
for(i in 1:360)
listfile[i]=paste("images/",i,".png",sep="")
outdirr=paste("../..",getwd(),"../../../animation.html",sep="")
im.convert(listfile,output=outdirr)
}

figure.ani=function(prova,startyear,colorss,phi,r,d,expand,ltheta,lphi,shade,titlep,bgcolor,nanni,color,facetcol)
{
ani.start()
par(bg=bgcolor)
for(i in 1:360)
{
persp(seq(1,365),seq(startyear,startyear+nanni-1),prova,col=color[facetcol],theta=i,phi=phi,r=r,d=d,expand=expand,ltheta=ltheta,lphi=lphi,shade=shade,ticktype="detailed" , nticks=5,xlab=" \n\nGiorni",ylab="\n\nAnni",zlab="\n\nPortata (mc/s)",cex.axis=1,cex.lab=1.5,border=NA)
if(is.na(titlep)==F)
title(titlep,cex.main=1.8)
}
ani.stop()
}

delbis=function(serie,bis)
{
        n=length(serie)
        nanni=trunc(n/365)
        ser=c()
        nb=floor(((nanni-bis)/4)+1)
        cat("bisestili",nb,fill=T)
        serie[(bis-1)*365+60]=-100000
        for(i in 1:(nb-1)) {
         serie[(bis-1)*365+60+i*1461]=-100000
        }
        j=0
        for(i in 1:(n-nb)) {
         if(serie[i+j]==-100000) {
          j=j+1
         }
         ser[i]=serie[i+j]
        }
        return(ser)
}