https://www.rayshader.com/ https://rallydatajunkie.com/visualising-rally-stages/annotating-rayshader-maps.html

# carga los paquetes requeridos
x <- c("sf","ggplot2","rgdal","dplyr","raster","OpenStreetMap","tmap","tmaptools","broom","tidyverse","RStoolbox")
lapply(x, library, character.only = TRUE)
#generate a list of input rasters ("grids")
#pattern = "*.tif$" - filters for main raster files only and skips any associated files (e.g. world files)
grids <- list.files("D:/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales/CEM_V3_20170619_R15_E05_TIF/" , pattern = "*.tif$")

#create a raster stack from the input raster files 
COAH_RASTER <- 
  brick("D:/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales/CEM_V3_20170619_R15_E05_TIF/Coahuila_r15m.tif")

#######
# Reading raster files
rutaras<-"D:/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales/CEM_V3_20170619_R15_E05_TIF/"
rst <- lapply(list.files(rutaras , pattern = "*.tif$",full.names = T), FUN = brick)

rst_ext <- lapply(rst, FUN =  extent)

# Calculating full extent
xmin_rst <- c(); xmax_rst <- c(); ymin_rst <- c(); ymax_rst <- c();
for (i in 1:length(rst_ext)) {
    xmin_rst <- c(xmin_rst, rst_ext[[i]]@xmin)
    ymin_rst <- c(ymin_rst, rst_ext[[i]]@ymin)
    xmax_rst <- c(xmax_rst, rst_ext[[i]]@xmax)
    ymax_rst <- c(ymax_rst, rst_ext[[i]]@ymax)
}
full_extent <- extent(min(xmin_rst), max(xmax_rst),
                      min(ymin_rst), max(ymax_rst))

# Creating raster from full extent and first rasters' CRS and resolution
bounding_rst <- raster(full_extent, 
                   crs = crs(rst[[1]]),
                   res = res(rst[[1]]))

# Resampling rasters to match attributes of the bounding raster
rst_resampled <- lapply(X = rst, FUN = function(x) {
    target_rst <- crop(bounding_rst, x)
    resample(x, target_rst, method="bilinear")
})

# Creating mosaic
rst_mosaic <- do.call("mosaic", c(rst_resampled, fun = mean))

######
sp::proj4string(rst_mosaic) <- CRS("+proj=longlat +datum=WGS84 +no_defs")

plot(rst_mosaic)

rutacam<-"D:/Documents/Claudia/Midropbox/Investigacion y escritos/mapasnacionales/"
limitesbox<-read.csv(paste0(rutacam,"limites.txt"),header=TRUE, sep=",",encoding="latin")
names(limitesbox)<-c("long","lat")
p = Polygon(limitesbox)
ps = Polygons(list(p),1)
sps = SpatialPolygons(list(ps))
sp::proj4string(sps) <- CRS("+proj=longlat +datum=WGS84 +no_defs")
crs(sps)
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]
recorterasterinegi <- crop(rst_mosaic, extent(sps), snap="out") 

library(rayshader)
library(sp)
library(scales)
memory.limit(size=1000000)
## [1] 1e+06
recorterasterinegi <- crop(rst_mosaic, extent(sps), snap="out") 

crs(recorterasterinegi)
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]
crs(rst_mosaic)
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]
elmat = raster_to_matrix(recorterasterinegi) 

elmat %>%
  sphere_shade(texture = "desert", progbar=FALSE) %>%
  plot_map()

elev_img<-recorterasterinegi
heightmap = matrix(
  raster::extract(elev_img, raster::extent(elev_img), method = 'bilinear'),
  nrow = ncol(elev_img),
  ncol = nrow(elev_img)
)
plot(raster(heightmap))

heightmap_raster  = raster( t(heightmap) )

plot(heightmap_raster)

plot(heightmap_raster,
     col = grey.colors(10, start=0, end=1) )

extent(heightmap_raster) = extent(elev_img)
crs(heightmap_raster) = crs(elev_img)

plot(heightmap_raster)

ray_shade_array = ray_shade(heightmap, progbar=FALSE)

ray_shade_array %>% 
  plot_map()

rasterify_rayshade = function(rayshade, original_raster){
  # Transpose and flip the shade matrix
  reorient_rayshade = fliplr(t( rayshade ))
  # Generate the raster
  rayshade_raster  = raster( reorient_rayshade )
  # Set the extent
  extent(rayshade_raster) = extent(original_raster)
  # Set the CRS
  crs(rayshade_raster) = crs(original_raster)
  
  rayshade_raster
}
fliplr = function(x) {
  if(length(dim(x)) == 2) {
    x[,ncol(x):1]
  } else {
    x[,ncol(x):1,]
  }
}
ray_shade_raster = rasterify_rayshade(ray_shade_array, elev_img)
ray_shade_raster
## class      : RasterLayer 
## dimensions : 2490, 4101, 10211490  (nrow, ncol, ncell)
## resolution : 0.0001388889, 0.0001388889  (x, y)
## extent     : -103.9247, -103.3551, 25.22643, 25.57227  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)
plot(ray_shade_raster)

plot(ray_shade_raster,
     col = grey.colors(10, start=0, end=1))

auto_zscale<-8
raymat <- ray_shade(elmat, zscale = auto_zscale, lambert = TRUE)
# And "ambient occlusion shadow"
ambmat <- ambient_shade(elmat, zscale = auto_zscale)

shadow_layer = height_shade(elmat) %>% 
  add_overlay(sphere_shade(elmat, texture = "desert", 
                           zscale=auto_zscale, colorintensity = 5), alphalayer=0.5) 



watermap <- detect_water(elmat, zscale = 8)
elmat %>%
  #sphere_shade(texture = "imhof1") %>%
  # Optionally set sun angle
  sphere_shade(sunangle = -45, texture = "unicorn") %>%
  add_water(watermap, color = "bw") %>%
  plot_map()

library(geoviz)

# stamen: toner, watercolor, terrain
overlay_image_watercolor <-
  slippy_overlay(elev_img, 
                 image_source = "stamen",
                 image_type = "watercolor",
                 png_opacity = 0.5)

overlay_image_tonerlite <-
  slippy_overlay(elev_img, 
                 image_source = "stamen",
                 image_type = "toner",
                 png_opacity = 0.5)

overlay_image_terrain <-
  slippy_overlay(elev_img, 
                 image_source = "stamen",
                 image_type = "terrain",
                 png_opacity = 0.9)

scene <- elmat %>%
  sphere_shade(sunangle = -45, texture = "unicorn") %>% 
  add_overlay(overlay_image_watercolor)

scene  %>%
  plot_map()

scene <- elmat %>%
  sphere_shade(sunangle = -45, texture = "desert") %>% 
  add_overlay(overlay_image_terrain)

scene  %>%
  plot_map()

scene <- elmat %>%
  sphere_shade(sunangle = 270, texture = "unicorn") %>% 
  add_overlay(overlay_image_tonerlite)

scene  %>%
  plot_map()

overlay_image_tonerlite_low_alt <-
  elevation_transparency(overlay_image_terrain,
                        elev_img,
                        pct_alt_high = 0.5,
                        alpha_max = 0.9)

scene <- elmat %>%
  sphere_shade(sunangle = 270, texture = "bw") %>% 
  add_overlay(overlay_image_tonerlite_low_alt) %>%
  plot_map()

contour_color = "#ffffff" #"#7d4911"

contours_layer = generate_contour_overlay(elmat, color = contour_color, linewidth = 1,levels=seq(1000,                    1300,by=1))


contoured_map <- elmat %>%
  sphere_shade(sunangle = 270, texture = "bw") %>% 
  add_overlay(contours_layer, alphalayer = 0.9)
 

contoured_map %>%  
   plot_map()

library(osmdata)

stage_osm  = opq(sf::st_bbox(sps[1,])) %>% 
  add_osm_feature("highway") %>% 
  osmdata_sf()

stage_lines = stage_osm$osm_lines

stage_roads = stage_lines[stage_lines$highway %in%
                                 c("unclassified", "secondary", "tertiary", "residential", "service","primary"), ] 


roads_layer = generate_line_overlay(stage_roads,
                                    extent =extent(elev_img),
                                    color="green",
                                    heightmap = elmat)



stage_osm  = opq(sf::st_bbox(sps[1,])) %>% 
  add_osm_feature("waterway") %>% 
  osmdata_sf()

stage_lines = stage_osm$osm_lines

stage_roads = stage_lines[stage_lines$waterway %in%
          c("stream", "river", "canal"), ] 
                            #c("unclassified", "secondary", "tertiary", "residential", "service","primary"), ] 


roads_water = generate_line_overlay(stage_roads,
                                    extent =extent(elev_img),
                                    linewidth=10,
                                    color="navyblue",
                                    heightmap = elmat)

    #00004d
contoured_map <- elmat %>%
  sphere_shade(sunangle = 270, texture = "desert") %>% 
   add_overlay(overlay_image_tonerlite)%>% 
   add_overlay(roads_layer)  %>%  
   add_overlay(roads_water)  %>%
  plot_map()

contoured_map <- elmat %>%
  sphere_shade(sunangle = 270, texture = "desert") %>% 
   add_overlay(overlay_image_watercolor)%>% 
   add_overlay(roads_layer)  %>%
   add_overlay(roads_water)  %>%
  plot_map()

contoured_map <- elmat %>%
sphere_shade(sunangle = 270, texture = "bw") %>% 
   add_overlay(overlay_image_tonerlite_low_alt)%>% 
   add_overlay(roads_layer)  %>%
   add_overlay(roads_water)  %>%

  plot_map()

#<https://jollydata.blog/posts/2022-06-05-ny-times-bestsellers/>
#add_overlay(overlay_image_tonerlite)
#overlay_image_tonerlite_low_alt

??generate_line_overlay