labsMap.R

require(maptools)
require(dplyr)
require(googledrive)
require(leaflet)
require(leaflet.extras)

# https://stackoverflow.com/questions/31336898/how-to-save-leaflet-in-r-map-as-png-or-jpg-file
library(devtools)
# install_github("wch/webshot")
library(htmlwidgets)
library(htmltools)
library(webshot)
# webshot::install_phantomjs()

# Obteniendo la fecha actual
fecha <- Sys.time() %>% format.Date("%Y%m%d_%H%M%p")

#' https://gist.github.com/briatte/18a4d543d1ccca194b2a03ac512be2b4
#' https://stackoverflow.com/questions/19497652/reading-kml-files-into-r
#' Read Points out of a KML file.
#' 
#' @param x A KML file exported from Google Maps.
#' @param layer The name of the layer to extract from: defaults to \code{"d1"}.
#' @param verbose Whether to report invalid coordinates and/or altitudes below
#' sea level; defaults to \code{TRUE}. See \link{kml_coordinate}.
#' @return A \link[tibble:tibble]{tibble} containing the \code{name}, 
#' \code{description}, \code{styleUrl} and Point coordinates (\code{longitude}, 
#' \code{latitude} and \code{altitude}) of each Placemark element contained in 
#' the KML file. Placemark elements with no Point coordinates, such as Polygon
#' elements, will be discarded.
#' @seealso \url{https://developers.google.com/kml/documentation/kmlreference}
kml_points <- function(x, layer = "d1", verbose = TRUE) {
  
  require(dplyr)
  require(stringr)
  require(xml2)
  
  #' Extract Placemark fields.
  #' 
  #' @param x A nodeset of Placemarks.
  #' @param field The name of the field to extract, e.g. \code{"name"}.
  #' @param layer The name of the layer to extract from; defaults to \code{"d1"}.
  #' @return A character vector. Missing values, i.e. empty fields, will be
  #' returned as \code{NA} values.
  get_field <- function(x, field, layer = "d1") {
    
    # vectorization required to get missing values when field is xml_missing
    lapply(x, xml_find_first, str_c(layer, ":", field)) %>%
      sapply(xml_text)
    
  }
  
  x <- read_xml(x) %>%
    xml_find_all(str_c("//", layer, ":Point/.."))
  
  x <- data_frame(
    name = get_field(x, "name", layer),
    description = get_field(x, "description", layer),
    styleUrl = get_field(x, "styleUrl", layer),
    coordinates = get_field(x, str_c("Point/", layer, ":coordinates"), layer)
  )
  
  x$longitude <- kml_coordinate(x$coordinates, 1, verbose)
  x$latitude  <- kml_coordinate(x$coordinates, 2, verbose)
  x$altitude  <- kml_coordinate(x$coordinates, 3, verbose)
  
  return(select(x, -coordinates))
  
}

#' Read Polygons out of a KML file.
#' 
#' @param x A KML file exported from Google Maps.
#' @param layer The name of the layer to extract from: defaults to \code{"d1"}.
#' @param verbose Whether to report invalid coordinates and/or altitudes below
#' sea level; defaults to \code{TRUE}. See \link{kml_coordinate}.
#' @return A \link[tibble:tibble]{tibble} containing the \code{name}, 
#' \code{description}, \code{styleUrl} and Point coordinates (\code{longitude}, 
#' \code{latitude} and \code{altitude}) of each Placemark element contained in 
#' the KML file. Placemark elements with no Point coordinates, such as Polygon
#' elements, will be discarded.
#' @seealso \url{https://developers.google.com/kml/documentation/kmlreference}
kml_polygons <- function(x, layer = "d1", verbose = TRUE) {
  
  require(dplyr)
  require(stringr)
  require(xml2)
  
  #' Extract Placemark fields.
  #' 
  #' @param x A nodeset of Placemarks.
  #' @param field The name of the field to extract, e.g. \code{"name"}.
  #' @param layer The name of the layer to extract from; defaults to \code{"d1"}.
  #' @return A character vector. Missing values, i.e. empty fields, will be
  #' returned as \code{NA} values.
  get_field <- function(x, field, layer = "d1") {
    
    xml_find_first(x, str_c(layer, ":", field)) %>%
      xml_text
    
  }
  
  x <- read_xml(x) %>%
    xml_find_all(str_c("//", layer, ":Polygon/.."))
  
  x <- lapply(x, function(x) {
    data_frame(
      name = get_field(x, "name"),
      description = get_field(x, "description"),
      styleUrl = get_field(x, "styleUrl"),
      coordinates = get_field(x, str_c("Polygon//", layer, ":coordinates")) %>%
        str_split(" ") %>%
        unlist
    )
  }) %>%
    bind_rows
  
  x$longitude <- kml_coordinate(x$coordinates, 1)
  x$latitude  <- kml_coordinate(x$coordinates, 2)
  x$altitude  <- kml_coordinate(x$coordinates, 3)
  
  return(select(x, -coordinates))
  
}

#' Extract KML coordinates
#' 
#' @param x A character vector of KML coordinates, of the form
#' \code{"longitude,latitude,altitude"}.
#' @param coord Which coordinate to extract: either \code{1} (longitude), 
#' \code{2} (latitude) or \code{3} (altitude).
#' @param verbose Whether to report invalid coordinates and/or altitudes below
#' sea level; defaults to \code{TRUE}.
#' @return A numeric vector.
#' @seealso \url{https://developers.google.com/kml/documentation/kmlreference}
# kml_coordinate <- function(x, coord, verbose = TRUE) {
#   
#   require(stringr) # includes `%>%`
#   
#   x <- str_replace(x, "(.*),(.*),(.*)", str_c("\\", coord)) %>%
#     as.numeric
#   
#   if (verbose && coord == 1 && any(abs(x) > 180))
#     message("Some longitudes are not contained within [-180, 180].")
#   
#   if (verbose && coord == 2 && any(abs(x) > 90))
#     message("Some latitudes are not contained within [-90, 90].")
#   
#   if (verbose && coord == 3 && any(x < 0))
#     message("Some altitudes are below sea level.")
#   
#   return(x)
#   
# }
# Usando función para extraer puntos de un archivo kml.
# https://gist.github.com/briatte/18a4d543d1ccca194b2a03ac512be2b4
# puntos_kml <- kml_points("Map_01.kml")
# puntos_kml <- puntos_kml %>% 
#   select(c("name", "description", "longitude", "latitude")) %>% 
#   mutate(name = name %>% str_replace("\n", ""))

puntos_atonal <- kml_points("puntos_atonal.kml") %>% 
    select(c("name", "description", "longitude", "latitude")) %>%
    mutate(name = name %>% str_replace("\n", ""), comunidad = "atonal")

puntos_tutunichapa <- kml_points("puntos_tutunichapa.kml") %>% 
  select(c("name", "description", "longitude", "latitude")) %>%
  mutate(name = name %>% str_replace("\n", ""), comunidad = "tutunichapa")

puntos_asuncion <- kml_points("puntos_asuncion.kml") %>% 
  select(c("name", "description", "longitude", "latitude")) %>%
  mutate(name = name %>% str_replace("\n", ""), comunidad = "asuncion")

puntos_santa_fe <- kml_points("puntos_santa_fe.kml") %>% 
  select(c("name", "description", "longitude", "latitude")) %>%
  mutate(name = name %>% str_replace("\n", ""), comunidad = "santa_fe")

puntos_kml <- rbind(puntos_atonal, puntos_asuncion, puntos_santa_fe, puntos_tutunichapa)

# Mapeo preliminar de puntos
# mapa_puntos <- leaflet(puntos_kml) %>% 
# addTiles() %>% 
# addMarkers(~longitude, ~latitude)
mapa_puntos <- leaflet(puntos_kml) %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addMarkers(~longitude, ~latitude, popup = ~htmlEscape(name))

# Guarda archivo .csv de los puntos
write.csv(puntos_kml, file="puntos_kml_final.csv")

# Generación de mapa web
mapafilename01 <- paste0("mapa_puntos",fecha)
saveWidget(mapa_puntos, paste0(mapafilename01,".html"), selfcontained = FALSE)
webshot(paste0(mapafilename01,".html"), file = paste0(mapafilename01,".png"),
        cliprect = "viewport")

# Mapa de calor Leaflet

# Observaciones
url1 <-  ""
url2 <-  ""
url3 <-  ""

drive_download(file = as_id(url1),
               path = "observaciones.csv",
               overwrite = TRUE)

drive_download(file = as_id(url2),
               path = "leyenda.csv",
               overwrite = TRUE)

drive_download(file = as_id(url3),
               path = "escenarios.csv",
               overwrite = TRUE)

# Procesando observaciones
situacion <- read.csv("situacion.csv")
situacion <- tbl_df(situacion)

observaciones <- read.csv("observaciones.csv")
observaciones <- tbl_df(observaciones)
observaciones <- observaciones %>% 
  select(c("comunidad", "nombre", "ubicacion", "puntaje", "participantes", "situacion"))

observaciones <- full_join(observaciones,situacion)
# Procesando leyenda
puntos_final <- read.csv("puntos_kml_final.csv")

# Recolectando los puntos relacionados a riesgos naturales
riesgos_naturales <- puntos_final %>% 
  filter(lugar %in% c("Sismo", "Fuego", "Desechos", "Inundacion"))
riesgos_naturales <- tbl_df(riesgos_naturales)

leyenda <- read.csv("leyenda.csv",header=T,encoding="UTF-8")
leyenda <- tbl_df(leyenda)
leyenda <- leyenda %>% 
    select(comunidad,ubicacion,nombre,longitude,latitude)

leyenda_full <- left_join(leyenda,puntos_final) %>% 
  select(comunidad,ubicacion,nombre,longitude,latitude)

# Uniendo datos
datos_joined <- full_join(observaciones, leyenda_full) %>% 
  mutate(punto = paste0(comunidad,ubicacion)) %>% 
  mutate(calificacion = puntaje/participantes)

datos_joined <- full_join(datos_joined, situacion)

datos_avg <- datos_joined %>% 
 group_by(punto, tipo, nombre) %>% 
 summarize(
   promedio = mean(calificacion),
   latitude=mean(latitude),
   longitude=mean(longitude)
 )

# range01 <- function(x){(x-min(x, na.rm = T))/(max(x, na.rm = T)-min(x, na.rm = T))}

# datos_avg <- datos_avg %>% mutate(rango = range01(datos_avg$promedio))


# https://rpubs.com/bhaskarvk/leaflet-heatmap
# 
datos_avg <- datos_avg %>% na.omit

# https://stackoverflow.com/questions/38274803/heatmap-colors-on-leaflet-in-r
m <- leaflet(datos_avg) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addWebGLHeatmap(
    lng=~longitude,
    lat=~latitude,
    intensity=~promedio,
    size = 182,
    opacity = 0.4,
  )

m_asuncion <- leaflet(datos_avg %>% filter(nombre=="atonal")) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addWebGLHeatmap(
    lng=~longitude,
    lat=~latitude,
    intensity=~promedio,
    size = 120,
    opacity = 0.4,
  )

m_tutunichapa <- leaflet(datos_avg %>% filter(nombre=="tutunichapa")) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addWebGLHeatmap(
    lng=~longitude,
    lat=~latitude,
    intensity=~promedio,
    size = 120,
    opacity = 0.4,
  )

m_santa_fe <- leaflet(datos_avg %>% filter(nombre=="santa_fe")) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addWebGLHeatmap(
    lng=~longitude,
    lat=~latitude,
    intensity=~promedio,
    size = 120,
    opacity = 0.4,
  )

m_atonal <- leaflet(datos_avg %>% filter(nombre=="atonal")) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addWebGLHeatmap(
    lng=~longitude,
    lat=~latitude,
    intensity=~promedio,
    size = 120,
    opacity = 0.4,
  )
# 
mapafilename02 <- paste0("mapa_calor",fecha)

saveWidget(m, paste0(mapafilename02, ".html"), selfcontained = FALSE)
saveWidget(m_atonal, paste0(mapafilename02, "_atonal.html"), selfcontained = FALSE)
saveWidget(m_tutunichapa, paste0(mapafilename02, "_tutunichapa.html"), selfcontained = FALSE)
saveWidget(m_santa_fe, paste0(mapafilename02, "_santa_fe.html"), selfcontained = FALSE)
saveWidget(m_asuncion, paste0(mapafilename02, "_asuncion.html"), selfcontained = FALSE)

webshot(paste0(mapafilename02, ".html"), file = paste0(mapafilename02, ".png"),
        cliprect = "viewport")
webshot(paste0(mapafilename02, "_atonal.html"), file = paste0(mapafilename02, "_atonal.png"),
        cliprect = "viewport")
webshot(paste0(mapafilename02, "_tutunichapa.html"), file = paste0(mapafilename02, "_tutunichapa.png"),
        cliprect = "viewport")
webshot(paste0(mapafilename02, "_santa_fe.html"), file = paste0(mapafilename02, "_santa_fe.png"),
        cliprect = "viewport")
webshot(paste0(mapafilename02, "_asuncion.html"), file = paste0(mapafilename02, "_asuncion.png"),
        cliprect = "viewport")

# 
riesgos_asuncion <- riesgos_naturales %>% filter(nombre=="asuncion")
color_asuncion <- getColor(riesgos_asuncion)
riesgos_asuncion <- cbind(riesgos_asuncion,color_asuncion)
m_asuncion_n <- leaflet(riesgos_asuncion) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircles(~longitude, ~latitude, color=color_asuncion, radius=3, fillOpacity=0.6, stroke=F, label=~as.character(lugar))

riesgos_atonal <- riesgos_naturales %>% filter(nombre=="atonal")
color_atonal <- getColor(riesgos_atonal)
riesgos_atonal <- cbind(riesgos_atonal,color_atonal)
m_atonal_n <- leaflet(riesgos_atonal) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircles(~longitude, ~latitude, color=color_atonal, radius=3, fillOpacity=0.6, stroke=F, label=~as.character(lugar))

riesgos_tutunichapa <- riesgos_naturales %>% filter(nombre=="tutunichapa")
color_tutunichapa <- getColor(riesgos_tutunichapa)
riesgos_tutunichapa <- cbind(riesgos_tutunichapa,color_tutunichapa)
m_tutunichapa_n <- leaflet(riesgos_tutunichapa) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircles(~longitude, ~latitude, color=color_tutunichapa, fillOpacity=0.6, radius=3, stroke=F, label=~as.character(lugar))

riesgos_santa_fe <- riesgos_naturales %>% filter(nombre=="santa_fe")
color_santa_fe <- getColor(riesgos_santa_fe)
riesgos_santa_fe <- cbind(riesgos_santa_fe,color_santa_fe)
m_santa_fe_n <- leaflet(riesgos_santa_fe) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircles(~longitude, ~latitude, color=color_santa_fe, radius=4, fillOpacity=0.6, stroke=F, label=~as.character(lugar))

mapafilename03 <- paste0("mapa_natural",fecha)

saveWidget(m_atonal_n, paste0(mapafilename03, "_atonal.html"), selfcontained = FALSE)
saveWidget(m_tutunichapa_n, paste0(mapafilename03, "_tutunichapa.html"), selfcontained = FALSE)
saveWidget(m_santa_fe_n, paste0(mapafilename03, "_santa_fe.html"), selfcontained = FALSE)
saveWidget(m_asuncion_n, paste0(mapafilename03, "_asuncion.html"), selfcontained = FALSE)