#distance in kilometers between two long/lat positions (from "fossil" package)
earth.dist <- function (long1, lat1, long2, lat2)
{
    rad <- pi/180
    a1 <- lat1 * rad
    a2 <- long1 * rad
    b1 <- lat2 * rad
    b2 <- long2 * rad
    dlon <- b2 - a2
    dlat <- b1 - a1
    a <- (sin(dlat/2))^2 + cos(a1) * cos(b1) * (sin(dlon/2))^2
    c <- 2 * atan2(sqrt(a), sqrt(1 - a))
    R <- 6378.145
    d <- R * c
    return(d)
}

#################################################################################
#conv lon lat
new.lon.lat <-
function (lon, lat, bearing, distance)
{
    rad <- pi/180
    a1 <- lat * rad
    a2 <- lon * rad
    tc <- bearing * rad
    d <- distance/6378.145
    nlat <- asin(sin(a1) * cos(d) + cos(a1) * sin(d) * cos(tc))
    dlon <- atan2(sin(tc) * sin(d) * cos(a1), cos(d) - sin(a1) *
        sin(nlat))
    nlon <- ((a2 + dlon + pi)%%(2 * pi)) - pi
    npts <- cbind(nlon/rad, nlat/rad)
    return(npts)
}


#################################################################################
#A function that is sometimes useful in determining the 
#coordinate(i.e. row and column number) of a matrix position
#(and vice-versa). 
#Either a vector of positions ("pos") 
#OR a 2 column matrix of matrix coordinates, ("coord", i.e. cbind(row,col)), 
#AND the matrix dimentions must be supplied (dim.mat, i.e. c(nrow,ncol)).
pos2coord<-function(pos=NULL, coord=NULL, dim.mat=NULL){
 if(is.null(pos) & is.null(coord) | is.null(dim.mat)){
  stop("must supply either 'pos' or 'coord', and 'dim.mat'")
 }
 if(is.null(pos) & !is.null(coord) & !is.null(dim.mat)){
  pos <- ((coord[,2]-1)*dim.mat[1])+coord[,1]
  return(pos)
 }
 if(!is.null(pos) & is.null(coord) & !is.null(dim.mat)){
  coord <- matrix(NA, nrow=length(pos), ncol=2)
  coord[,1] <- ((pos-1) %% dim.mat[1]) +1
  coord[,2] <- ((pos-1) %/% dim.mat[1]) +1
  return(coord)
 }
}

#################################################################################
#interpolation of xyz data and projection onto a map
#modif par guillaume pour ajouter une scale
map.xyz.scale <- function(xyz_data, file_name="map.xyz.png",
 
projection="stereographic", orientation=NULL,
 
par=NULL, fill=FALSE, nside=40, breaks=100, res=200,
 
xlim = NULL, ylim = NULL, legend= "", titlelab="" 
 
){
 
 
 if(length(which(rowSums(is.na(xyz_data)) != 0)) > 0){
 
  xyz_data <- xyz_data[-which(rowSums(is.na(xyz_data)) != 0),]
 
 }
temp<-interp(
 
  xyz_data[,1], xyz_data[,2], xyz_data[,3],  
 
  xo=seq(min(xyz_data[,1]), max(xyz_data[,1]), length = nside), 
 
  yo=seq(min(xyz_data[,2]), max(xyz_data[,2]), length = nside, extrap=TRUE)
 
 )
 
polys<-vector("list", length(as.vector(temp$z)))
 
 for(i in 1:length(polys)){
 
  lonx <- pos2coord(pos=i, dim.mat=dim(temp$z))[1]
 
  laty <- pos2coord(pos=i, dim.mat=dim(temp$z))[2]
 
  ifelse(laty < length(temp$y), neigh_y<-c(laty+1,lonx), neigh_y<-c(laty-1,lonx))
 
  ifelse(lonx < length(temp$x), neigh_x<-c(laty,lonx+1), neigh_x<-c(laty,lonx-1)) 
 
  dist_y <- earth.dist(temp$x[lonx], temp$y[laty], temp$x[neigh_y[2]], temp$y[neigh_y[1]])
 
  dist_x <- earth.dist(temp$x[lonx], temp$y[laty], temp$x[neigh_x[2]], temp$y[neigh_x[1]])
 
  s1 = new.lon.lat(temp$x[lonx], temp$y[laty], 180, dist_y/2)
 
  s3 = new.lon.lat(temp$x[lonx], temp$y[laty], 0, dist_y/2)
 
  s2 = new.lon.lat(temp$x[lonx], temp$y[laty], 270, dist_x/2)
 
  s4 = new.lon.lat(temp$x[lonx], temp$y[laty], 90, dist_x/2)
 
  polys[[i]] = cbind(c(s2[1], s2[1], s4[1], s4[1]), c(s1[2], s3[2], s3[2], s1[2]))
 
 }

 M.range=range(temp$z, na.rm=TRUE)
 
 M.breaks <- pretty(M.range, n=breaks)
 
 M.cols=colorRampPalette(c("blue","cyan", "yellow", "red"),space="rgb")
 
 colorlut <- M.cols(length(M.breaks)) # color lookup table
 
 colorvalues <- colorlut[((as.vector(temp$z)-M.breaks[1])/(range(M.breaks)[2]-range(M.breaks)[1])*length(M.breaks))+1] # assign colors to heights for each point
 
 if(is.null(xlim)){xlim <- range(xyz_data[,1], na.rm=TRUE)}
 
 if(is.null(ylim)){ylim <- range(xyz_data[,2], na.rm=TRUE)}
 
#figure settings
heights=c(4)
widths=c(4,1)

 png(filename=file_name, res=res,width = sum(widths), height = sum(heights), units="in")
 par(omi=c(0.1, 0.1, 0.1, 0.1), ps=12)
 layout(matrix(c(1,2),nrow=1,ncol=2,byrow=TRUE), widths = widths, heights = heights, respect=TRUE)
 layout.show(2)

 par(mai=c(0.2, 0.2, 0.2, 0.2)) 

 map("world",projection=projection, orientation=orientation, par=par, fill=fill, xlim=xlim, ylim=ylim)
 
 for(i in 1:length(polys)){
 
  polygon(mapproject(x=polys[[i]][,1], y=polys[[i]][,2]), col=colorvalues[i], border=NA)
 
#  print(i)
 
 }

 map("world", add=TRUE, projection="", par=par, fill=fill, xlim=xlim, ylim=ylim)
 
 map.grid(c(-180, 180, -80, 80), nx=10, ny=18, labels=FALSE, col="grey")

 title(main=titlelab)

 pal=M.cols
 #add scale
 par(mai=c(0.2, 0, 0.2, 0.6))
 image.scale(temp$z, col=pal(100), horiz=FALSE, yaxt="n")
 axis(4,las=1)
 mtext(legendlab, side=4, line=2.5)
 box()
 dev.off()
 
}

#################################################################################
#interpolation of xyz data and projection onto a map
map.xyz <- function(xyz_data, file_name="map.xyz.png",
 
projection="stereographic", orientation=NULL,
 
par=NULL, fill=FALSE, nside=40, breaks=100, res=100,
 
xlim = NULL, ylim = NULL
 
){
 
 
 if(length(which(rowSums(is.na(xyz_data)) != 0)) > 0){
 
  xyz_data <- xyz_data[-which(rowSums(is.na(xyz_data)) != 0),]
 
 }

 temp<-interp(
 
  xyz_data[,1], xyz_data[,2], xyz_data[,3],  
 
  xo=seq(min(xyz_data[,1]), max(xyz_data[,1]), length = nside), 
 
  yo=seq(min(xyz_data[,2]), max(xyz_data[,2]), length = nside, extrap=TRUE)
 
 )
 
 
 
 polys<-vector("list", length(as.vector(temp$z)))
 
 for(i in 1:length(polys)){
 
  lonx <- pos2coord(pos=i, dim.mat=dim(temp$z))[1]
 
  laty <- pos2coord(pos=i, dim.mat=dim(temp$z))[2]
 
  ifelse(laty < length(temp$y), neigh_y<-c(laty+1,lonx), neigh_y<-c(laty-1,lonx))
 
  ifelse(lonx < length(temp$x), neigh_x<-c(laty,lonx+1), neigh_x<-c(laty,lonx-1)) 
 
  dist_y <- earth.dist(temp$x[lonx], temp$y[laty], temp$x[neigh_y[2]], temp$y[neigh_y[1]])
 
  dist_x <- earth.dist(temp$x[lonx], temp$y[laty], temp$x[neigh_x[2]], temp$y[neigh_x[1]])
 
  s1 = new.lon.lat(temp$x[lonx], temp$y[laty], 180, dist_y/2)
 
  s3 = new.lon.lat(temp$x[lonx], temp$y[laty], 0, dist_y/2)
 
  s2 = new.lon.lat(temp$x[lonx], temp$y[laty], 270, dist_x/2)
 
  s4 = new.lon.lat(temp$x[lonx], temp$y[laty], 90, dist_x/2)
 
  polys[[i]] = cbind(c(s2[1], s2[1], s4[1], s4[1]), c(s1[2], s3[2], s3[2], s1[2]))
 
 }
M.range=range(temp$z, na.rm=TRUE)
 
 M.breaks <- pretty(M.range, n=breaks)
 
 M.cols=colorRampPalette(c("blue","cyan", "yellow", "red"),space="rgb")
 
 colorlut <- M.cols(length(M.breaks)) # color lookup table
 
 colorvalues <- colorlut[((as.vector(temp$z)-M.breaks[1])/(range(M.breaks)[2]-range(M.breaks)[1])*length(M.breaks))+1] # assign colors to heights for each point
 
 if(is.null(xlim)){xlim <- range(xyz_data[,1], na.rm=TRUE)}
 
 if(is.null(ylim)){ylim <- range(xyz_data[,2], na.rm=TRUE)}
 
 png(filename=file_name, res=res)
 
 map("world",projection=projection, orientation=orientation, par=par, fill=fill, xlim=xlim, ylim=ylim)
 
 for(i in 1:length(polys)){
 
  polygon(mapproject(x=polys[[i]][,1], y=polys[[i]][,2]), col=colorvalues[i], border=NA)
 
  print(i)
 
 }
 
 map("world", add=TRUE, projection="", par=par, fill=fill, xlim=xlim, ylim=ylim)
 
 map.grid(c(-180, 180, -80, 80), nx=10, ny=18, labels=FALSE, col="grey")
 
 dev.off()
 
}




#################################################################################
#this function converts a vector of values("z") to a vector of color
#levels. One must define the number of colors. The limits of the color
#scale("zlim") or the break points for the color changes("breaks") can 
#also be defined. when breaks and zlim are defined, breaks overrides zlim.

val2col<-function(z, zlim, col = heat.colors(12), breaks){
 if(!missing(breaks)){
  if(length(breaks) != (length(col)+1)){stop("must have one more break than colour")}
 }
 if(missing(breaks) & !missing(zlim)){
  zlim[2] <- zlim[2]+c(zlim[2]-zlim[1])*(1E-3)#adds a bit to the range in both directions
  zlim[1] <- zlim[1]-c(zlim[2]-zlim[1])*(1E-3)
  breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1)) 
 }
 if(missing(breaks) & missing(zlim)){
  zlim <- range(z, na.rm=TRUE)
  zlim[2] <- zlim[2]+c(zlim[2]-zlim[1])*(1E-3)#adds a bit to the range in both directions
  zlim[1] <- zlim[1]-c(zlim[2]-zlim[1])*(1E-3)
  breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1))
 }
 colorlevels <- col[((as.vector(z)-breaks[1])/(range(breaks)[2]-range(breaks)[1]))*(length(breaks)-1)+1] # assign colors to heights for each point
 colorlevels
}


#################################################################################
#This function creates a color scale for use with the image()
#function. Input parameters should be consistent with those
#used in the corresponding image plot. The "horiz" argument
#defines whether the scale is horizonal(=TRUE) or vertical(=FALSE).
image.scale <- function(z, zlim, col = heat.colors(12),
breaks, horiz=TRUE, ...){
 if(!missing(breaks)){
  if(length(breaks) != (length(col)+1)){stop("must have one more break than colour")}
 }
 if(missing(breaks) & !missing(zlim)){
  breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1)) 
 }
 if(missing(breaks) & missing(zlim)){
  zlim <- range(z, na.rm=TRUE)
  zlim[2] <- zlim[2]+c(zlim[2]-zlim[1])*(1E-3)#adds a bit to the range in both directions
  zlim[1] <- zlim[1]-c(zlim[2]-zlim[1])*(1E-3)
  breaks <- seq(zlim[1], zlim[2], length.out=(length(col)+1))
 }
 poly <- vector(mode="list", length(col))
 for(i in seq(poly)){
  poly[[i]] <- c(breaks[i], breaks[i+1], breaks[i+1], breaks[i])
 }
 xaxt <- ifelse(horiz, "s", "n")
 yaxt <- ifelse(horiz, "n", "s")
 if(horiz){ylim<-c(0,1); xlim<-range(breaks)}
 if(!horiz){ylim<-range(breaks); xlim<-c(0,1)}
 plot(1,1,t="n",ylim=ylim, xlim=xlim, xaxt=xaxt, yaxt=yaxt, xaxs="i", yaxs="i", ...)  
 for(i in seq(poly)){
  if(horiz){
   polygon(poly[[i]], c(0,0,1,1), col=col[i], border=NA)
  }
  if(!horiz){
   polygon(c(0,0,1,1), poly[[i]], col=col[i], border=NA)
  }
 }
}