library(dplyr)
library(DT)
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:raster':
## 
##     select
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(readxl) 

library(leaflet)  # libreria para graficar mapas interactivos
library(sf)  # manejo de informacion geografica 
# library(viridis)  # paletas de colores
# library(RColorBrewer)  # mas paletas de colores 
library(rio)  
## 
## Attaching package: 'rio'
## The following object is masked from 'package:plotly':
## 
##     export
library(bs4Dash)
## Warning: package 'bs4Dash' was built under R version 4.1.2
## 
## Attaching package: 'bs4Dash'
## The following object is masked from 'package:graphics':
## 
##     box
library(fresh)
## Warning: package 'fresh' was built under R version 4.1.2
library(sf)
library(rgl)
## Warning: package 'rgl' was built under R version 4.1.2
library(rgdal)
library(viridis)
## Loading required package: viridisLite
## 
## Attaching package: 'viridis'
## The following object is masked from 'package:scales':
## 
##     viridis_pal
library(viridisLite)
library(rayshader)
library(ggplot2)
library(rglwidget)
# image <- ee$Image("LANDSAT/LC08/C01/T1/LC08_232090_20211105")
# image <- image$select("B[2-4]")
# roi <- ee$Geometry$Point(list(-72.50630061337516, -42.83997103647681))$buffer(25000) 
# clipped <- image$clip(roi)
# 
# print(clipped$getInfo())
# Map$setCenter(-72.50630061337516, -42.83997103647681)
# 
# vizParams <- list(
#       bands = c("B4", "B3", "B2"),
#       min = 5000,
#       max = 15000,
#       gamma = 1.3
#     )
#     
# # Map$addLayer(image, vizParams, "Full Image", FALSE) +
# Map$addLayer(clipped, vizParams, "Clipped Image")
# ee_raster <- ee_as_raster(
#   image = clipped,
#   region = roi,
#   dsn = "D:/TEMP/glaciar_1136.tif",
#   scale = 30,
  # via = "drive"
# )
library(raster)
# img_1 = raster("D:/TEMP/glaciar_1136.tif")
# plot(img_1)
Band4 <- raster("D:/TEMP/glaciar_1136.tif",band = 3)
Band3 <- raster("D:/TEMP/glaciar_1136.tif",band = 2)
Band2 <- raster("D:/TEMP/glaciar_1136.tif",band = 1)

zion_rbg = raster::stack(Band2, Band4, Band3)
# raster::plotRGB(zion_rbg, scale=255^2)
# GDALinfo("D:/TEMP/glaciar_1136.tif")
# dataset <- ee$Image('CGIAR/SRTM90_V4')
# elevation <- dataset$select('elevation')
# zion_elevation <- elevation$clip(roi)

# Download raster
# ee_raster <- ee_as_raster(
#   image = zion_elevation,
#   region = roi,
#   dsn = "D:/TEMP/glaciar_1136_elevation.tif",
#   scale = 30,
#   via = "drive"
# )
zion_elevation = raster::raster("D:/TEMP/glaciar_1136_elevation.tif")
# zion_elevation
zion_rbg_corrected = sqrt(raster::stack(Band4, Band3, Band2))
# raster::plotRGB(zion_rbg_corrected)
zion_elevation_utm = raster::projectRaster(zion_elevation, crs = crs(Band4), method = "bilinear")
crs(zion_elevation_utm)
## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs

-42.717705319603624, -72.71819889148851 -42.93983872425845, -72.25107596486319

coordenadas <- extent(zion_elevation_utm)
coordenadas
## class      : Extent 
## xmin       : 678086.4 
## xmax       : 729720.4 
## ymin       : -4771983 
## ymax       : -4720226
coordenadas[1]
## [1] 678086.4
coordenadas[2]
## [1] 729720.4
coordenadas[3]
## [1] -4771983
coordenadas[4]
## [1] -4720226
# bottom_left = c(y=-72.71819889148851 , x=-42.717705319603624)
# top_right   = c(y=-72.25107596486319,  x=-42.93983872425845)
bottom_left = c(y=-72.53754801990983 , x=-42.74105240679547)
top_right   = c(y=-72.32279177110796,  x=-42.88092599413714)
# bottom_left = c(y=coordenadas[2] , x=coordenadas[3])
# top_right   = c(y=coordenadas[1],  x=coordenadas[4])

extent_latlong = sp::SpatialPoints(rbind(bottom_left, top_right), proj4string=sp::CRS("+proj=longlat +zone=18  +ellps=WGS84 +datum=WGS84 +units=m +no_dfs +towgs84=0,0,0"))
extent_utm = sp::spTransform(extent_latlong, raster::crs(zion_elevation_utm))

e = raster::extent(extent_utm)
e
## class      : Extent 
## xmin       : 701556.2 
## xmax       : 718641.7 
## ymin       : -4751069 
## ymax       : -4735000
zion_rgb_cropped = raster::crop(zion_rbg_corrected, e)
elevation_cropped = raster::crop(zion_elevation_utm, e)
names(zion_rgb_cropped) = c("r","g","b")

zion_r_cropped = rayshader::raster_to_matrix(zion_rgb_cropped$r)
zion_g_cropped = rayshader::raster_to_matrix(zion_rgb_cropped$g)
zion_b_cropped = rayshader::raster_to_matrix(zion_rgb_cropped$b)

zionel_matrix = rayshader::raster_to_matrix(elevation_cropped)

zion_rgb_array = array(0,dim=c(nrow(zion_r_cropped),ncol(zion_r_cropped),3))

zion_rgb_array[,,1] = zion_r_cropped/255 #Red layer
zion_rgb_array[,,2] = zion_g_cropped/255 #Blue layer
zion_rgb_array[,,3] = zion_b_cropped/255 #Green layer
# 
zion_rgb_array = aperm(zion_rgb_array, c(2,1,3))

# plot_map(zion_rgb_array)
zion_rgb_contrast = scales::rescale(zion_rgb_array,to=c(0,1))

# plot_map(zion_rgb_contrast)
plot_3d(zion_rgb_contrast, zionel_matrix, windowsize = c(1080,720), zscale = 14, shadowdepth = -50,
        zoom=0.5, phi=27,theta=-45,fov=70, background = "#FFFFFF", shadowcolor = "#523E2B")
render_snapshot(title_text = "Prueva vista de Iquique | Imagery: Landsat 8 | DEM: 30m SRTM",
                title_bar_color = "#1f5214", title_color = "white", title_bar_alpha = 1)

rglwidget()