#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
#	Druksisme, portrait topographique
#	d'après une idée originale de Michael Druks (1974-1975)
#	langage : R 2.9
#	auteur : procrastin@fxdm.org
#	date de dernière modification : 06 avril 2010
#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

###
### chargement des packages
###
#install.packages('rgdal')
library(rgdal)  
#install.packages('akima')
library(akima)

###
### importation de l'image carrée 288x288 pixels, conversion et rotation
###
imdata <- readGDAL("procrastin.ppm")
green.matrix <- as.matrix(imdata[2])
# rotation de la matrice
green.matrix <- green.matrix[,ncol(green.matrix):1]

###
### conversion de la matrice en coordonnées X, Y, Z
###
X <- numeric(288*288)
Y <- numeric(288*288)
Z <- numeric(288*288)
ligne <- 1
for ( x in 1:288) {
	for (y in 1:288) {
		X[ligne] <- x
		Y[ligne] <- y
		Z[ligne] <- green.matrix[x,y]
		ligne <- ligne + 1
	}
}
carte.li <- interp(X, Y, Z)

###
### tracé de la carte
###
windows(width = 20, height = 21.38, pointsize = 12)
colcode <- c(heat.colors(25),"#cbdddf")
filled.contour(carte.li$x,carte.li$y,carte.li$z, 
	asp = 1, xlim = c(5,282), ylim = c(12,350),	
	col = colcode, xaxs="i", yaxs="i",	
	plot.axes={ 
		## dessin des axes x et y # axis(1); axis(2); 
		axis(1, at = 150, lab = "174°O", cex = 2);
		axis(2, at = 200, lab = "47°S", srt = 45, cex = 5)
		for (i in seq(0, 350, 50)) { 
			axis(1, at = i, lab = "", cex = 2);
			axis(2, at = i, lab = "", cex = 2);
		}
		## coloriage des marges en dehors de l'image proprement dite
		polygon(c(0, 0, 288,288) , c(288, 350, 350, 288),  col ="#cbdddf", lty = 0, border = NULL);
		polygon(c(-20, -20, 1,1) , c(-10, 350, 350, -10),  col ="#cbdddf", lty = 0, border = NULL);
		polygon(c(288,288, 350, 350) , c(-10, 350, 350, -10),  col ="#cbdddf", lty = 0, border = NULL);
		## dessin des lignes de niveaux
		contour(carte.li$x,carte.li$y,carte.li$z, lwd = 0.5, asp = 1, add=TRUE);		
		## quelques sommets
		points(180, 195, pch = 24, col = "black", bg = "black", cex = 0.6);
		points(172, 126, pch = 24, col = "black", bg = "black", cex = 0.6);	
		## repères marins aléatoires
		pas <- 25; x <- 0; y <- 0; i <- 1;
		while (i <=100) {
			x <- sample(seq(0, 282, 5), 1)
			y <- sample(seq(0, 350, 5), 1)
			if ((x>288)||(y>288)||(length(matr[x,y]) == 0)||( matr[x,y] == 255)) {
				text(x, y, sample(10:80, 1), col ="grey", cex = 0.5); i <- i +1;
			}
		}		
		## tracé des liaisons maritimes 
		x <- c(221, 225, 250, 300); y <- c(75, 75, 25, 20); xspline(x, y, open = T, shape = 1, border = "#666b94", lty = 3)
		x <- c(221, 250, 300); y <- c(75, 90, 150); xspline(x, y, open = T, shape = 1, border = "#666b94", lty = 3)
		x <- c(-10, 50, 100); y <- c(250, 270, 255); xspline(x, y, open = T, shape = 1, border = "#666b94", lty = 3)
		x <- c(100,10, 150,200); y <- c(255, 300, 375, 400); xspline(x, y, open = T, shape = 1, border = "#666b94", lty = 3)
		## échelle noire et blanche le long des axes
		largeur <- 1; long <- 50; i <- 0; coulech = "black";
		while (i <= 350){
			polygon(c(0, 0, largeur, largeur), c(i, i+long, i+long, i),  col =coulech, lwd = 0.25); 
			polygon(c(i, i+long, i+long, i), c(12, 12, 12 +largeur, 12 + largeur),  col =coulech, lwd = 0.25);
			polygon(c(i, i+long, i+long, i), c(350, 350, 350-largeur, 350- largeur),  col =coulech, lwd = 0.25);
			x <- 287; polygon(c(x, x, x-largeur, x-largeur), c(i, i+long, i+long, i),  col =coulech, lwd = 0.25); 
			if (coulech =="black") { coulech = "white"} else {coulech = "black"}; i <- i + long;
		}
		## ticks le long des axes
		longueur <- 1; espacement <- 2; i <- 0; coulech = "grey";
		polygon(c(largeur, largeur, largeur + longueur,  largeur + longueur), c(0, 350, 350, 0), col = "white", lwd = 0.25)
		polygon( c(0, 350, 350, 0), c(12+largeur, 12+largeur, 12+largeur + longueur,  12+largeur + longueur), col = "white", lwd = 0.25)
		polygon( c(0, 350, 350, 0), c(350-largeur, 350-largeur, 350-largeur -longueur,  350-largeur -longueur), col = "white", lwd = 0.25)
		polygon(c(284+largeur, 284+largeur, 284+largeur + longueur,  284+largeur + longueur), c(0, 350, 350, 0), col = "white", lwd = 0.25)
		while (i <= 350){
			segments(largeur, i, largeur + longueur, i, lwd = 0.25) 
			segments(i, 12+largeur, i, 12+largeur + longueur, lwd = 0.25)
			segments(i, 350-largeur, i, 350-largeur - longueur, lwd = 0.25)
			segments(284+largeur, i, 284+largeur + longueur, i, lwd = 0.25) 
			i <- i + espacement;
		}			
		## parallele et longitude
		abline(v = 150, col = "#666b94");
		abline(h = 200, col = "#666b94");
	}
)