library(OpenStreetMap)
library(png)
library(raster)
library(rgl)
library(classInt)

clouds  <- function(ra, offset, mypal=c("#ffffff00","#0b13ff7f","#050773ff"), n = 10, style = "equal", fixed = NULL){
  
  z <- rotate(matrix(ra[], nrow = ra@nrows, ncol = ra@ncols, byrow = T))
  x <- seq(ra@extent@xmin, ra@extent@xmax, length.out = ra@ncols)
  y <- seq(ra@extent@ymin, ra@extent@ymax, length.out = ra@nrows)
  
  h5 <- classIntervals(as.numeric(z[]), n=n, style = style, fixedBreaks = fixed)
  h5Colours <- findColours(h5, mypal)
  opac <- as.numeric(z[])
  opac[is.na(opac)] <- 0
  opac <- opac/max(opac)
  
  surface3d(x, y, offset*z/z, color = h5Colours, alpha = opac, back = "lines")
}

get.CellWidth <- function(dem){ 
  ((xres(dem)*pi/180) * 6371000) #* cos(mean(c(dem@extent@ymin,dem@extent@ymax)) * pi/180)
}

rotate <- function(x) t(apply(x, 2, rev))
rotate_cc <- function(x) apply(t(x), 2, rev)

pan3d <- function(button) {
  start <- list()
  
  begin <- function(x, y) {
    start$userMatrix <<- par3d("userMatrix")
    start$viewport <<- par3d("viewport")
    start$scale <<- par3d("scale")
    start$projection <<- rgl.projection()
    start$pos <<- rgl.window2user( x/start$viewport[3], 1 - y/start$viewport[4], 0.5, 
                                   projection=start$projection)
  }
  
  update <- function(x, y) {
    xlat <- (rgl.window2user( x/start$viewport[3], 1 - y/start$viewport[4], 0.5,
                              projection = start$projection) - start$pos)*start$scale
    mouseMatrix <- translationMatrix(xlat[1], xlat[2], xlat[3])
    par3d(userMatrix = start$userMatrix %*% t(mouseMatrix) )
  }
  rgl.setMouseCallbacks(button, begin, update)
 # cat("Callbacks set on button", button, "of rgl device",rgl.cur(),"\n")
}

download.srtm <- function(region, v=1, dem_temp=NULL){
  if(is.null(dem_temp)) dem_temp <- tempfile(fileext=".tif")
  
  url <- paste0("http://ot-data1.sdsc.edu:9090/otr/getdem?demtype=SRTMGL",v,
                "&west=",region@xmin,
                "&south=",region@ymin,
                "&east=",region@xmax,
                "&north=",region@ymax,"&outputFormat=GTiff")
  download.file(url = url, destfile = dem_temp)
  dem <- raster(dem_temp)
  return(dem)
}

dem3d <- function(dem, sf=1, texture=F, zoom=NULL, type="bing", tmp_texture=NULL, col=NULL){

  # z <- rotate(as.matrix(dem))
  z <- rotate(matrix(dem[], nrow = dem@nrows, ncol = dem@ncols, byrow = T))
  x <- seq(dem@extent@xmin, dem@extent@xmax, length.out = dem@ncols)
  y <- seq(dem@extent@ymin, dem@extent@ymax, length.out = dem@nrows)

  # open3d()
  if(!texture){
    if(is.null(col)){
      z[z[]<=0 | is.na(z[])] <- 1
      zlim <- range(z[])
      zlen <- zlim[2] - zlim[1] + 1
      colorlut <- terrain.colors(zlen) # height color lookup table
      col <- colorlut[z] # assign colors to heights for each point
    }
    # surface3d(x, y, z*xres(dem)/sf, color = col, back = "lines")
    persp3d(x, y, z, color = col, back = "lines",
            aspect=c(1,length(y)/length(x),diff(range(z))/(get.CellWidth(dem)/sf)/length(x)),
            axes=F, xlab="", ylab = "", zlab="")
  }else{
    if(!file.exists(paste0(tmp_texture,""))){
      if(is.null(zoom))stop("zoom = ??? Specify zoomlevel!")
      tmp_texture <- load.texture(dem, zoom, type, tmp_texture)
    }
    # surface3d(x, y, z*xres(dem)/sf, color = "white", texture=tmp_texture, texmipmap=T, back = "lines")
    persp3d(x, y, z, color = "white", back = "lines",
            # aspect=c(length(x),length(y),1/xres(dem)/sf),
            aspect=c(1,length(y)/length(x),diff(range(z, na.rm = T))/(get.CellWidth(dem)/sf)/length(x)),
            axes=F, xlab="", ylab = "", zlab="",
            texture=tmp_texture, texmipmap=T)
  }
  pan3d(3)
}

load.texture <- function(dem, zoom, type="bing", tmp_texture=NULL){
  map <- openmap(c(dem@extent@ymax, dem@extent@xmin),c(dem@extent@ymin, dem@extent@xmax), zoom, type)
  map <- openproj(map)
  
  map <- matrix(map$tiles[[1]]$colorData, nrow = map$tiles[[1]]$yres, ncol = map$tiles[[1]]$xres)
  map <- map[rowSums(is.na(map))!=ncol(map), ]
  map <- map[ ,colSums(is.na(map))!=nrow(map)]
  map <- rotate_cc(map)
  
  if(is.null(tmp_texture)) tmp_texture <- tempfile(fileext=".png")
  
  rgba <- array(dim = c(dim(map), 4))
  rgb.data <- col2rgb(map)
  rgba[, , 4] <- 1
  rgba[, , 4][is.na(map)] <- 0
  rgba[, , 1] <- rgb.data[1, ]/255
  rgba[, , 2] <- rgb.data[2, ]/255
  rgba[, , 3] <- rgb.data[3, ]/255
  
  
  writePNG(image=rgba,target = tmp_texture)
  return(tmp_texture)
}

shape3d <- function(shape,name,pal=c("darkgreen","orange","red")){
  
  n <- length(shape@data[,1])
  
  class <- classIntervals(shape@data[,name], n=n, style="equal")
  colf <- findColours(class, pal=pal)
  
  open3d()
  
  
  for(i in 1:n){
    
    lat <- shape@polygons[[i]]@Polygons[[1]]@coords[,2]
    long <- shape@polygons[[i]]@Polygons[[1]]@coords[,1]
    
    lat.m <- matrix(lat, ncol=length(lat), nrow=2, byrow=T)
    long.m <- matrix(long, ncol=length(lat), nrow=2, byrow=T)
    z.m <- matrix(c(rep(0,length(lat)),rep(shape@data[i,name],length(lat))), nrow=2, byrow=T)
    if(i==1){
      persp3d(long.m, lat.m, z.m, col=colf[i], aspect=c(1,diff(range(lat))/diff(range(long)),1),
              xlab="", ylab="", zlab="", axes=F)
    }else{
      persp3d(long.m, lat.m, z.m, col=colf[i], add=TRUE)
    }
    
    # use gpc.poly for fast triangulation!
    xy <- as(cbind(long,lat), "gpc.poly")
    triangles <- triangulate(xy)
    triangles3d(triangles[,1],triangles[,2], z.m[2,1], col=colf[i])
    
    # shorter code, but to slow
    # shade3d( extrude3d(lat, long, shape@data[i,name]/100000), col = "red" )
    
  }
  pan3d(3)
}

raster.to.png <- function (raster, name = "test", mypal = c("wheat1", "red3", "green"),
                           n = 10, style = "equal", fixed = NULL, trans = 1, legend = FALSE, flip=F){
  png.path <- getOption("png.path", default = getwd())
  matrix <- matrix(raster[], nrow = raster@nrows, ncol = raster@ncols, byrow = TRUE)
  
  if(flip) matrix <- matrix[nrow(matrix):1,]
  
  h5 <- classIntervals(as.numeric(matrix[]), n = n, style = style, fixedBreaks = fixed)
  h5Colours <- findColours(h5, mypal)
  if (legend == TRUE) {
    pdf(paste(png.path, "/", name, "_legend.pdf", sep = ""))
    plot(1, axes = FALSE, ann = F, type="n")
    l.text <- gsub(","," - ",gsub("\\]","",gsub("\\)","",gsub("\\[","",names(attr(h5Colours,"table"))))))
    legend("top", col = attr(h5Colours, "palette"), pch=15, pt.cex = 2.5,
           legend = l.text, text.col ="grey10",
           bg = "grey90", box.col = "grey80")
    dev.off()
  }
  
  rgba <- array(dim = c(dim(matrix), 4))
  rgba[, , 4] <- trans
  rgba[, , 4][is.na(matrix)] <- 0
  rgb.data <- col2rgb(h5Colours)
  rgba[, , 1] <- rgb.data[1, ]/255
  rgba[, , 2] <- rgb.data[2, ]/255
  rgba[, , 3] <- rgb.data[3, ]/255
  writePNG(rgba, target = paste0(png.path, "/", name, ".png"))
  
  write.table(c(xres(raster),0,0,-yres(raster),extent(raster)[1]+xres(raster)/2,extent(raster)[1]-yres(raster)/2),
              paste0(png.path, "/", name, ".pngw"), col.names = F, row.names = F)
}