#########
# Datas #
#########
spaceObject="mars"
datas=read.table("../data/3-mars/megt90n000cb_xyz.txt") #chosen map
value0=1436 #chosen value for the sea level. Here with about 70% of the surface covered by water
valuesInKilometre=FALSE #unit for the datas file. False = meter, true = kilometer

datas=datas[order(datas$V2),] #we put in order to have increasing longitude and latitude vectors
datas=as.matrix(datas)
x=unique(datas[,1]) #longitude
y=unique(datas[,2]) #latitude
z=matrix(datas[,3],ncol=length(y)) #height for each ordered pair (x,y)

##########################################
# filled.contour function without legend #
##########################################
#The filled.contour function already exists. We take it and comment out the legend.
#filled.contour(x, y, z, col=colors, levels=intervals, axes=FALSE, ylab="",xlab="", key.title=FALSE)
filled.contour.nolegend = function (x = seq(0, 1, length.out = nrow(z)), y = seq(0, 1, 
    length.out = ncol(z)), z, xlim = range(x, finite = TRUE), 
    ylim = range(y, finite = TRUE), zlim = range(z, finite = TRUE), 
    levels = pretty(zlim, nlevels), nlevels = 20, color.palette = cm.colors, 
    col = color.palette(length(levels) - 1), plot.title, plot.axes, 
    key.title, key.axes, asp = NA, xaxs = "i", yaxs = "i", las = 1, 
    axes = TRUE, frame.plot = axes, ...) 
{
    if (missing(z)) {
        if (!missing(x)) {
            if (is.list(x)) {
                z <- x$z
                y <- x$y
                x <- x$x
            }
            else {
                z <- x
                x <- seq.int(0, 1, length.out = nrow(z))
            }
        }
        else stop("no 'z' matrix specified")
    }
    else if (is.list(x)) {
        y <- x$y
        x <- x$x
    }
    if (any(diff(x) <= 0) || any(diff(y) <= 0)) 
        stop("increasing 'x' and 'y' values expected")
    #mar.orig <- (par.orig <- par(c("mar", "las", "mfrow")))$mar
    #on.exit(par(par.orig))
    #w <- (3 + mar.orig[2L]) * par("csi") * 2.54
    #layout(matrix(c(2, 1), ncol = 2L), widths = c(1, lcm(w)))
    #par(las = las)
    #mar <- mar.orig
    #mar[4L] <- mar[2L]
    #mar[2L] <- 1
    #par(mar = mar)
    #plot.new()
    #plot.window(xlim = c(0, 1), ylim = range(levels), xaxs = "i", 
    #    yaxs = "i")
    #rect(0, levels[-length(levels)], 1, levels[-1L], col = col)
    #if (missing(key.axes)) {
    #    if (axes) 
    #        axis(4)
    #}
    #else key.axes
    #box()
    #if (!missing(key.title)) 
    #    key.title
    #mar <- mar.orig
    #mar[4L] <- 1
    #par(mar = mar)
    plot.new()
    plot.window(xlim, ylim, "", xaxs = xaxs, yaxs = yaxs, asp = asp)
    if (!is.matrix(z) || nrow(z) <= 1L || ncol(z) <= 1L) 
        stop("no proper 'z' matrix specified")
    if (!is.double(z)) 
        storage.mode(z) <- "double"
    .filled.contour(as.double(x), as.double(y), z, as.double(levels), 
        col = col)
    if (missing(plot.axes)) {
        if (axes) {
            title(main = "", xlab = "", ylab = "")
            Axis(x, side = 1)
            Axis(y, side = 2)
        }
    }
    else plot.axes
    if (frame.plot) 
        box()
    if (missing(plot.title)) 
        title(...)
    else plot.title
    invisible()
}

##############################################
# Plot function with shifts of the sea level #
##############################################
coeffMultiplying=1
if(valuesInKilometre==TRUE)
{
 coeffMultiplying=1000
}

topography.plot <- function(x, y, z, shift)
{
 colorsSea=c("#71ABD8", "#79B2DE", "#84B9E3", "#8DC1EA", "#96C9F0", "#A1D2F7", "#ACDBFB", "#B9E3FF", "#C6ECFF", "#D8F2FE")
 colorsEarth=c("#ACD0A5", "#94BF8B", "#A8C68F", "#BDCC96", "#D1D7AB", "#E1E4B5", "#EFEBC0", "#E8E1B6", "#DED6A3", "#D3CA9D", "#CAB982", "#C3A76B", "#B9985A", "#AA8753", "#AC9A7C", "#BAAE9A", "#CAC3B8", "#E0DED8", "#F5F4F2")
 colors=c(colorsSea,colorsEarth)
 intervals=c(-100000, -3500, -3000, -2500, -2000, -1500, -1000, -750, -500, -200, 0, 50, 100, 250, 500, 750, 1000, 1250, 1500, 1750, 2000, 2250, 2500, 2750, 3000, 3250, 3500, 3750, 4000, 100000)
 intervals=intervals/coeffMultiplying
 par(mar=c(0,0,0,0))
 filled.contour.nolegend(x, y, z+shift, col=colors, levels=intervals, axes=FALSE)
 contour(x, y, z+shift, add=T,levels=c(0), col="#0978AB", drawlabels=FALSE)
}

###########
# Outputs #
###########
min(z) #for etopo10_bedrock.xyz : -10421
max(z) #for etopo10_bedrock.xyz : 6527
minChoice=-floor(-coeffMultiplying*min(z)/1000)*1000
maxChoice=floor(coeffMultiplying*max(z)/1000)*1000

if(spaceObject=="monde") #if space object is the world
{
 shiftVector=sort(unique(c(seq(from=minChoice, to=maxChoice, by=1000),seq(from=-1000, to=1000, by=100),seq(from=-100, to=100, by=10),value0)))
} else
{
 shiftVector=sort(unique(c(seq(from=minChoice, to=maxChoice, by=1000),value0)))  
}
shiftVector=shiftVector/coeffMultiplying

dir.create("./sorties480/")
dir.create("./sorties1024/")
dir.create("./sorties1280/")
dir.create("./sorties1920/")

for(i in 1:length(shiftVector))
{
 print(paste(i, "on", length(shiftVector)))
 shift=shiftVector[i] #increasing value of the sea level

 nameOutput=paste("./sorties480/", spaceObject, shift*coeffMultiplying, ".png", sep="")
 png(nameOutput, width=480, height=270)
 topography.plot(x,y,z,-shift)
 dev.off()

 nameOutput=paste("./sorties1024/", spaceObject, shift*coeffMultiplying, ".png", sep="")
 png(nameOutput, width=1024, height=576)
 topography.plot(x,y,z,-shift)
 dev.off()

 nameOutput=paste("./sorties1280/", spaceObject, shift*coeffMultiplying, ".png", sep="")
 png(nameOutput, width=1280, height=720)
 topography.plot(x,y,z,-shift)
 dev.off()

 nameOutput=paste("./sorties1920/", spaceObject, shift*coeffMultiplying, ".png", sep="")
 png(nameOutput, width=1920, height=1080)
 topography.plot(x,y,z,-shift)
 dev.off()
}