Mapterhorn + R: LOS analysis in seconds! 🚀

Posted byAlbertoPosted in5G, 5G telco, AI, analisis, Application, Cloud based computation, cnig, code/código, CV, geovisualizacion, GIS / SIG, gis analysis, IA, madrid, Modelado 3D, monitoring, news, noticias, open data, R, remote sensingTags:, , , , , , , , , , ,

The analysis of LOS (Line of Sight) in telecommunications is the study that determines whether a clear, unobstructed path exists between a transmitting antenna and a receiver. While this calculation was traditionally reserved for large microwave links over long distances or in rural environments, the arrival of 5G networks and the horizon of 6G have turned it into an absolute priority for urban deployment, completely transforming how networks are planned in major cities.

Diagram 1 – Mapterhorn, R and Geovisualization.net. A good team!!

This radical shift is due to the nature of the high frequencies used by next-generation technologies, particularly millimeter waves. Unlike older 4G networks, which operated on bands capable of bending around buildings and penetrating walls with relative ease, these new frequencies offer massive speeds but behave in an almost optical manner, similar to the beam of a laser pointer. Any physical obstacle, no matter how small, completely absorbs or blocks the signal. This forces the industry to shift away from a model of large antennas on rooftops toward a hyper-dense network of micro-antennas distributed at street level on elements like lampposts, traffic lights, or bus shelters.

This radical shift is due to the nature of the high frequencies used by next-generation technologies, particularly millimeter waves, which behave in an almost optical manner and require a hyper-dense network of micro-antennas distributed at street level. To map this dense urban mesh in central Madrid, telecommunications engineers rely heavily on geographic information systems and open data, where the Centro Nacional de Información Geográfica (CNIG) plays a crucial role. CNIG offers high-resolution Digital Terrain Models (5m in this case) and LiDAR flights that replicate every building and tree in the capital down to the centimeter. A key factor for the telecommunications industry is that CNIG operates under an open data license that explicitly permits commercial use. This legal framework gives operators and tech companies the absolute freedom to integrate these government-funded geographical datasets into their commercial simulation software and network planning pipelines without legal friction, significantly reducing the cost and accelerating the deployment of new infrastructure.

Diagram 2 – Visible vs Blocked

To bridge the gap between these massive raw geographic datasets and practical engineering solutions, the industry relies on specialized open-source tools, where Mapterhorn stands out as an essential link. Mapterhorn acts as a powerful bridge between R—the statistical and data computing language—and spatial data, allowing engineers to ingest, clean, and manipulate heavy LiDAR files and raster data with high performance. Through this integration, data analysts can programmatically automate the extraction of urban profiles and line-of-sight paths directly within their data pipelines.

Diagram 3 – LOS in Madrid downtown

This synergy between CNIG’s open commercial-ready data and the processing capabilities of Matterhorn within the R ecosystem enables the calculation of thousands of simultaneous ray-tracing simulations. It allows engineers to predict exactly how the signal will bounce off a glass facade on Gran Vía or diffract around a narrow street corner, making the connectivity of the future viable and stable at street level, despite the architectural difficulties of central Madrid.

Diagram 4 – Fotostock mate

Please let me know if this likes you!

# Desarrollado por Alberto Concejal / Geovisualization.net
# LOS Analyzer — Madrid | MDT05 IGN/CNIG via WCS
# Dependencias: shiny, terra, leaflet, sf
# install.packages(c("shiny","terra","leaflet","sf"))
library(shiny)
library(terra)
library(leaflet)
library(sf)
library(jsonlite)
# ─────────────────────────────────────────────────────────────────────────────
# FUNCIONES AUXILIARES
# ─────────────────────────────────────────────────────────────────────────────
#' Descarga MDT05 del IGN via WCS 1.0.0 para un bounding box dado (EPSG:4326)
#' Devuelve un SpatRaster en EPSG:25830 (UTM30) — sin reproyectar para evitar
#' conflicto con proj.db de PostgreSQL/PostGIS
descargar_mdt05 <- function(lon_min, lat_min, lon_max, lat_max, buffer_deg = 0.005) {
  # Buffer en grados
  xmin <- lon_min - buffer_deg
  ymin <- lat_min - buffer_deg
  xmax <- lon_max + buffer_deg
  ymax <- lat_max + buffer_deg
  # Convertir bbox a UTM30 via sf (sf usa su propia PROJ, funciona bien)
  bbox_geo <- st_bbox(c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax),
                      crs = 4326)
  bbox_utm <- st_bbox(st_transform(st_as_sfc(bbox_geo), 25830))
  # Resolución MDT05: 5 m/pixel — limitar a 800px máx
  res_m     <- 5
  width_px  <- min(ceiling((bbox_utm["xmax"] - bbox_utm["xmin"]) / res_m), 800)
  height_px <- min(ceiling((bbox_utm["ymax"] - bbox_utm["ymin"]) / res_m), 800)
  url <- paste0(
    "https://servicios.idee.es/wcs-inspire/mdt",
    "?SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCoverage",
    "&COVERAGE=Elevacion25830_5",
    "&CRS=EPSG:25830",
    "&BBOX=", bbox_utm["xmin"], ",", bbox_utm["ymin"], ",",
               bbox_utm["xmax"], ",", bbox_utm["ymax"],
    "&WIDTH=",  width_px,
    "&HEIGHT=", height_px,
    "&FORMAT=GEOTIFFINT16"
  )
  tmp <- tempfile(fileext = ".tif")
  tryCatch(
    download.file(url, tmp, mode = "wb", quiet = TRUE),
    error = function(e) stop("Error descargando MDT05: ", conditionMessage(e))
  )
  # Verificar que no es un XML de error
  first_line <- readLines(tmp, n = 1, warn = FALSE)
  if (grepl("xml|ServiceException", first_line, ignore.case = TRUE)) {
    stop("El WCS devolvió un error: ", paste(readLines(tmp), collapse = " "))
  }
  r <- tryCatch(
    rast(tmp),
    error = function(e) stop("Error leyendo GeoTIFF del WCS: ", conditionMessage(e))
  )
  # Asignar extensión UTM30 manualmente (el WCS devuelve ENGCRS)
  # No reproyectamos — trabajamos en UTM30 durante todo el análisis
  ext(r)  <- c(bbox_utm["xmin"], bbox_utm["xmax"],
               bbox_utm["ymin"], bbox_utm["ymax"])
  return(r)  # SpatRaster en coordenadas UTM30 (metros)
}
#' Extrae perfil de elevación entre dos puntos WGS84
#' El DEM está en UTM30 — se convierten los puntos con sf, no el raster
#' Devuelve data.frame: distance, elevation, lon, lat
perfil_elevacion <- function(dem, lon_a, lat_a, lon_b, lat_b, n = 300) {
  # Interpolar puntos en WGS84 y convertir a UTM30
  lons <- seq(lon_a, lon_b, length.out = n)
  lats <- seq(lat_a, lat_b, length.out = n)
  pts_sf  <- st_as_sf(data.frame(lon = lons, lat = lats),
                       coords = c("lon", "lat"), crs = 4326)
  pts_utm <- st_transform(pts_sf, 25830)
  coords  <- st_coordinates(pts_utm)
  # Extraer elevación del DEM en coordenadas UTM30
  # extract devuelve 1 columna con nombre = nombre del archivo temporal
  elev_raw <- terra::extract(dem, coords)
  elev <- elev_raw[, ncol(elev_raw)]
  # Distancia acumulada (m)
  dx   <- diff(coords[, 1])
  dy   <- diff(coords[, 2])
  dist <- c(0, cumsum(sqrt(dx^2 + dy^2)))
  data.frame(
    idx       = seq_len(n),
    distance  = dist,
    elevation = elev,
    lon       = lons,
    lat       = lats
  )
}
#' Calcula LOS entre A y B sobre un perfil dado
#' Devuelve lista con: visible, bloqueo_dist, bloqueo_elev, angle_deg
calcular_los <- function(perfil, h_a, h_b) {
  n     <- nrow(perfil)
  elev  <- perfil$elevation
  dist  <- perfil$distance
  elev_a <- elev[1] + h_a
  elev_b <- elev[n] + h_b
  d_total <- dist[n]
  # NA check
  if (any(is.na(elev))) {
    elev[is.na(elev)] <- min(elev, na.rm = TRUE)
  }
  # Pendiente de la línea de visión directa (en elevación/distancia)
  slope <- (elev_b - elev_a) / d_total
  # Elevación de la LOS en cada punto intermedio
  elev_los <- elev_a + slope * dist
  # Diferencia: terreno - LOS en puntos intermedios (excluir A y B)
  diff_elev <- elev[2:(n-1)] - elev_los[2:(n-1)]
  visible       <- isTRUE(all(diff_elev <= 0))
  bloqueo_idx   <- if (!visible) which(diff_elev > 0)[1] + 1 else NA_integer_
  bloqueo_dist  <- if (!visible && !is.na(bloqueo_idx)) dist[bloqueo_idx] else NA_real_
  bloqueo_elev  <- if (!visible && !is.na(bloqueo_idx)) elev[bloqueo_idx] else NA_real_
  exceso_max    <- if (!visible) max(diff_elev, na.rm = TRUE) else 0
  # Ángulo de elevación A→B (grados)
  angle_deg <- atan2(elev_b - elev_a, d_total) * 180 / pi
  list(
    visible      = visible,
    bloqueo_dist = bloqueo_dist,
    bloqueo_elev = bloqueo_elev,
    exceso_max   = exceso_max,
    elev_a       = elev[1],
    elev_b       = elev[n],
    elev_los     = elev_los,
    angle_deg    = angle_deg,
    dist_total   = d_total
  )
}
#' Descarga edificios y vías OSM via Overpass en un buffer de la línea A-B
#' Devuelve lista: $edificios (sf polygons), $roads (sf lines)
descargar_osm <- function(lon_a, lat_a, lon_b, lat_b, buffer_m = 200) {
  linea     <- st_sfc(st_linestring(rbind(c(lon_a, lat_a), c(lon_b, lat_b))), crs = 4326)
  buf_wgs   <- st_transform(st_buffer(st_transform(linea, 25830), buffer_m), 4326)
  bb        <- st_bbox(buf_wgs)
  resultado <- list(edificios = NULL, roads = NULL)
  # Query 1: ways (buildings + highways)
  q_ways <- paste0(
    "[out:json][timeout:30];(",
    "way[building](", bb["ymin"],",",bb["xmin"],",",bb["ymax"],",",bb["xmax"],");",
    "way[highway](", bb["ymin"],",",bb["xmin"],",",bb["ymax"],",",bb["xmax"],");",
    ");out body;>;out skel qt;"
  )
  # Query 2: POI nodes
  fetch_overpass <- function(query) {
    tmp <- tempfile(fileext = ".json")
    download.file(
      paste0("https://overpass-api.de/api/interpreter?data=",
             utils::URLencode(query, reserved = TRUE)),
      tmp, mode = "wb", quiet = TRUE
    )
    jsonlite::fromJSON(readLines(tmp, warn = FALSE), simplifyVector = FALSE)$elements
  }
  tmp <- tempfile(fileext = ".json")
  tryCatch({
    elems_ways <- tryCatch(fetch_overpass(q_ways), error = function(e) list())
    # Ways y nodes de geometría — solo de elems_ways
    nodes_geom <- Filter(function(e) e$type == "node", elems_ways)
    ways_all   <- Filter(function(e) e$type == "way" && !is.null(e$nodes), elems_ways)
    # Node lookup table para reconstruir geometrías de ways
    node_df <- NULL
    if (length(nodes_geom) > 0) {
      m <- do.call(rbind, lapply(nodes_geom, function(n)
        c(id = as.numeric(n$id), lon = as.numeric(n$lon), lat = as.numeric(n$lat))))
      node_df <- as.data.frame(m)
    }
    # Helper: way nodes -> matrix of coords
    way_coords <- function(w) {
      if (is.null(node_df)) return(NULL)
      ids <- as.numeric(unlist(w$nodes))
      idx <- match(ids, node_df$id)
      idx <- idx[!is.na(idx)]
      if (length(idx) < 2) return(NULL)
      cbind(node_df$lon[idx], node_df$lat[idx])
    }
    # ── Edificios ──
    ways_bld <- Filter(function(w) !is.null(w$tags$building), ways_all)
    if (length(ways_bld) > 0) {
      polys <- lapply(ways_bld, function(w) {
        co <- way_coords(w)
        if (is.null(co) || !is.matrix(co) || nrow(co) < 4) return(NULL)
        tryCatch(st_polygon(list(co)), error = function(e) NULL)
      })
      polys <- Filter(Negate(is.null), polys)
      if (length(polys) > 0) {
        sf_bld <- st_sf(
          osm_id   = sapply(ways_bld[seq_along(polys)], function(w) w$id),
          name     = sapply(ways_bld[seq_along(polys)], function(w)
                      if (isTRUE(length(w$tags$name) > 0)) as.character(w$tags$name)[1] else ""),
          geometry = st_sfc(polys, crs = 4326)
        )
        resultado$edificios <- tryCatch(st_intersection(sf_bld, buf_wgs),
                                        error = function(e) sf_bld)
      }
    }
    # ── Vías ──
    ways_hwy <- Filter(function(w) !is.null(w$tags$highway), ways_all)
    if (length(ways_hwy) > 0) {
      lns <- lapply(ways_hwy, function(w) {
        co <- way_coords(w)
        if (is.null(co)) return(NULL)
        tryCatch(st_linestring(co), error = function(e) NULL)
      })
      lns <- Filter(Negate(is.null), lns)
      if (length(lns) > 0) {
        sf_hwy <- st_sf(
          osm_id   = sapply(ways_hwy[seq_along(lns)], function(w) w$id),
          name     = sapply(ways_hwy[seq_along(lns)], function(w)
                      if (isTRUE(length(w$tags$name) > 0)) as.character(w$tags$name)[1] else ""),
          highway  = sapply(ways_hwy[seq_along(lns)], function(w) as.character(w$tags$highway)[1]),
          geometry = st_sfc(lns, crs = 4326)
        )
        resultado$roads <- tryCatch(st_intersection(sf_hwy, buf_wgs),
                                    error = function(e) sf_hwy)
      }
    }
  }, error = function(e) message("Overpass error: ", conditionMessage(e)))
  return(resultado)
}
# ─────────────────────────────────────────────────────────────────────────────
# UI
# ─────────────────────────────────────────────────────────────────────────────
ui <- fluidPage(
  tags$head(
    tags$style(HTML("
      @import url('https://fonts.googleapis.com/css2?family=IBM+Plex+Mono:wght@400;600&family=IBM+Plex+Sans:wght@400;500&display=swap');
      /* ── Paleta coherente: gris medio uniforme ── */
      :root {
        --bg:        #1e2023;
        --bg2:       #252830;
        --bg3:       #2d3138;
        --border:    #363b44;
        --txt:       #c8cdd6;
        --txt-dim:   #6b7280;
        --accent-a:  #b8e830;
        --accent-b:  #30b8e8;
        --danger:    #e85030;
      }
      * { box-sizing: border-box; margin: 0; padding: 0; }
      body { font-family: 'IBM Plex Sans', sans-serif; background: var(--bg); color: var(--txt); }
      .app-header {
        background: var(--bg2); border-bottom: 1px solid var(--border);
        padding: 12px 20px; display: flex; align-items: baseline; gap: 14px;
      }
      .app-title {
        font-family: 'IBM Plex Mono', monospace; font-size: 14px;
        font-weight: 600; color: var(--accent-a); letter-spacing: .08em;
      }
      .app-sub { font-family: 'IBM Plex Mono', monospace; font-size: 11px; color: var(--txt-dim); }
      .main-wrap { display: flex; height: calc(100vh - 49px); overflow: hidden; }
      /* ── Paneles laterales ── */
      .lpanel {
        width: 272px; min-width: 272px; background: var(--bg2);
        border-right: 1px solid var(--border); overflow-y: auto;
        display: flex; flex-direction: column; position: relative;
      }
      .rpanel {
        width: 320px; min-width: 220px; background: var(--bg2);
        border-left: 1px solid var(--border); overflow-y: auto;
        display: flex; flex-direction: column; position: relative;
      }
      .sec { padding: 13px 14px; border-bottom: 1px solid var(--border); }
      .sec-lbl {
        font-family: 'IBM Plex Mono', monospace; font-size: 9px; font-weight: 600;
        color: var(--txt-dim); letter-spacing: .12em; text-transform: uppercase; margin-bottom: 10px;
      }
      /* Bloques A / B */
      .pt-block {
        background: var(--bg3); border: 1px solid var(--border); border-radius: 5px;
        padding: 10px; margin-bottom: 8px;
      }
      .pt-block.pta { border-left: 3px solid var(--accent-a); }
      .pt-block.ptb { border-left: 3px solid var(--accent-b); }
      .pt-lbl { font-family: 'IBM Plex Mono', monospace; font-size: 10px;
                font-weight: 600; margin-bottom: 7px; }
      .pta .pt-lbl { color: var(--accent-a); }
      .ptb .pt-lbl { color: var(--accent-b); }
      /* inputs */
      .form-group { margin-bottom: 6px; }
      label { font-size: 10px; color: var(--txt-dim); font-family: 'IBM Plex Mono', monospace;
              margin-bottom: 2px; display: block; }
      .form-control {
        background: var(--bg) !important; border: 1px solid var(--border) !important;
        color: var(--txt) !important; font-family: 'IBM Plex Mono', monospace !important;
        font-size: 11px !important; border-radius: 3px !important;
        padding: 4px 7px !important; height: 27px !important; width: 100% !important;
      }
      .form-control:focus { border-color: var(--accent-a) !important; box-shadow: none !important; }
      select.form-control { height: 27px !important; }
      .row-2 { display: flex; gap: 6px; }
      .row-2 .form-group { flex: 1; }
      /* botones */
      .btn-go {
        width: 100%; background: var(--accent-a); color: #0e0e0e; border: none;
        border-radius: 4px; font-family: 'IBM Plex Mono', monospace;
        font-size: 12px; font-weight: 600; padding: 9px; cursor: pointer;
        letter-spacing: .05em; transition: background .12s;
      }
      .btn-go:hover:not(:disabled) { background: #caf050; }
      .btn-go:disabled { background: var(--border); color: var(--txt-dim); cursor: not-allowed; }
      .btn-sec {
        width: 100%; background: transparent; color: var(--accent-b);
        border: 1px solid var(--accent-b); border-radius: 4px;
        font-family: 'IBM Plex Mono', monospace; font-size: 11px;
        padding: 6px; cursor: pointer; margin-top: 5px; transition: background .12s;
      }
      .btn-sec:hover { background: rgba(48,184,232,.08); }
      /* resultado */
      .res-block { background: var(--bg3); border: 1px solid var(--border); border-radius: 5px; padding: 12px; }
      .res-status { display: flex; align-items: center; gap: 9px; margin-bottom: 8px; }
      .dot { width: 9px; height: 9px; border-radius: 50%; flex-shrink: 0; }
      .dot.ok  { background: var(--accent-a); }
      .dot.ko  { background: var(--danger); }
      .dot.na  { background: var(--border); }
      .res-txt { font-family: 'IBM Plex Mono', monospace; font-size: 13px; font-weight: 600; }
      .res-txt.ok { color: var(--accent-a); }
      .res-txt.ko { color: var(--danger); }
      .res-txt.na { color: var(--txt-dim); }
      .kv { display: flex; justify-content: space-between; padding: 4px 0;
            border-top: 1px solid var(--border); font-size: 11px; }
      .k { color: var(--txt-dim); font-family: 'IBM Plex Mono', monospace; }
      .v { color: var(--txt); font-family: 'IBM Plex Mono', monospace; font-weight: 600; }
      /* hint */
      .hint {
        background: #1e2a10; border: 1px solid #3a5018; border-radius: 4px;
        padding: 7px 9px; font-family: 'IBM Plex Mono', monospace;
        font-size: 10px; color: #8ab040; margin-top: 7px; line-height: 1.5;
      }
      #click_st {
        font-family: 'IBM Plex Mono', monospace; font-size: 10px; color: var(--accent-a);
        padding: 3px 0; min-height: 16px; display: block;
      }
      /* mapa */
      .map-wrap { flex: 1; position: relative; }
      #mapa { width: 100%; height: 100%; }
      .profile-wrap { padding: 14px; overflow-y: auto; }
      /* handles redimensionar */
      .resize-handle-r {
        position: absolute; top: 0; right: -4px; width: 8px; height: 100%;
        cursor: ew-resize; z-index: 100; background: transparent;
      }
      .resize-handle-l {
        position: absolute; top: 0; left: -4px; width: 8px; height: 100%;
        cursor: ew-resize; z-index: 100; background: transparent;
      }
      .resize-handle-r:hover, .resize-handle-r.dragging,
      .resize-handle-l:hover, .resize-handle-l.dragging {
        background: rgba(184,232,48,0.2);
      }
      /* scrollbar */
      ::-webkit-scrollbar { width: 4px; }
      ::-webkit-scrollbar-track { background: var(--bg); }
      ::-webkit-scrollbar-thumb { background: var(--border); border-radius: 2px; }
      .shiny-output-error { color: var(--danger); font-size: 11px; font-family: monospace; padding: 8px; }
      .shiny-output-error:before { content: '✕ '; }
      .footer { padding: 12px 14px; margin-top: auto; border-top: 1px solid var(--border); }
      .footer p { font-family: 'IBM Plex Mono', monospace; font-size: 10px; margin-bottom: 2px; }
      .btn-export { width: 100%; }
    "))
  ),
  # Header
  div(class = "app-header",
    tags$span(class = "app-title", "LOS ANALYZER"),
    tags$span(class = "app-sub",   "Madrid · MDT05 CNIG/IGN · 5 m/px")
  ),
  div(class = "main-wrap",
    # ── Panel izquierdo ──────────────────────────────────────────────────────
    div(class = "lpanel", id = "lpanel",
      tags$div(class = "resize-handle", id = "resize-handle"),
      # Puntos
      div(class = "sec",
        tags$p(class = "sec-lbl", "Puntos"),
        div(class = "pt-block pta",
          tags$p(class = "pt-lbl", "PUNTO A — Observador"),
          div(class = "row-2",
            div(class = "form-group",
              tags$label("Latitud"),
              numericInput("lat_a", NULL, 40.4168, step = 0.00001)
            ),
            div(class = "form-group",
              tags$label("Longitud"),
              numericInput("lon_a", NULL, -3.7038, step = 0.00001)
            )
          ),
          div(class = "form-group",
            tags$label("Altura obs. (m s/suelo)"),
            numericInput("h_a", NULL, 1.7, min = 0, max = 300, step = 0.1)
          )
        ),
        div(class = "pt-block ptb",
          tags$p(class = "pt-lbl", "PUNTO B — Objetivo"),
          div(class = "row-2",
            div(class = "form-group",
              tags$label("Latitud"),
              numericInput("lat_b", NULL, 40.4200, step = 0.00001)
            ),
            div(class = "form-group",
              tags$label("Longitud"),
              numericInput("lon_b", NULL, -3.6900, step = 0.00001)
            )
          ),
          div(class = "form-group",
            tags$label("Altura obj. (m s/suelo)"),
            numericInput("h_b", NULL, 1.7, min = 0, max = 300, step = 0.1)
          )
        ),
        actionButton("toggle_clic", "⊕  Definir puntos en el mapa", class = "btn-sec"),
        tags$span(id = "click_st", textOutput("click_estado", inline = TRUE)),
        div(class = "hint",
          "Activa el modo clic → haz clic en el mapa para definir A y luego B."
        )
      ),
      # Parámetros
      div(class = "sec",
        tags$p(class = "sec-lbl", "Parámetros"),
        div(class = "row-2",
          div(class = "form-group",
            tags$label("Nº muestras perfil"),
            numericInput("n_pts", NULL, 300, min = 100, max = 1000, step = 50)
          ),
          div(class = "form-group",
            tags$label("Buffer descarga (°)"),
            numericInput("buf", NULL, 0.005, min = 0.002, max = 0.02, step = 0.001)
          )
        )
      ),
      # Calcular + Reset
      div(class = "sec",
        actionButton("calcular", "▶  CALCULAR LOS", class = "btn-go"),
        tags$div(style = "margin-top: 6px;",
          actionButton("resetear", "↺  Resetear vista", class = "btn-sec",
                       style = "color: var(--txt-dim); border-color: var(--border);")
        )
      ),
      # Resultado
      div(class = "sec",
        tags$p(class = "sec-lbl", "Resultado LOS"),
        div(class = "res-block",
          uiOutput("resultado_ui")
        )
      ),
      # Exportar
      div(class = "sec",
        tags$p(class = "sec-lbl", "Exportar"),
        downloadButton("exportar_gpkg", "⬇  Exportar GPKG",
          style = "width:100%; background: var(--bg3); color: var(--txt);
                   border: 1px solid var(--border); border-radius: 4px;
                   font-family: IBM Plex Mono, monospace; font-size: 11px;
                   padding: 7px; cursor: pointer; text-align: center;")
      ),
      div(class = "footer",
        tags$p(style = "color: var(--accent-a); margin-bottom: 3px;",
               "Developed by Alberto Concejal"),
        tags$p(style = "color: var(--txt-dim);",
               "Geovisualization.net")
      )
    ),
    # ── Mapa ─────────────────────────────────────────────────────────────────
    div(class = "map-wrap",
      leafletOutput("mapa", width = "100%", height = "100%")
    ),
    # ── Panel perfil ─────────────────────────────────────────────────────────
    div(class = "rpanel", id = "rpanel",
      tags$div(class = "resize-handle-l", id = "resize-handle-r"),
      div(class = "sec",
        tags$p(class = "sec-lbl", "Perfil de elevación")
      ),
      div(class = "profile-wrap",
        plotOutput("perfil_plot", height = "340px", width = "100%"),
        uiOutput("perfil_stats_ui")
      )
    )
  ),
  tags$script(HTML("
    (function() {
      var activeHandle = null, activePanel = null, dragging = false;
      var startX, startW, direction;
      function initResize(handleId, panelId, dir) {
        var h = document.getElementById(handleId);
        var p = document.getElementById(panelId);
        if (!h || !p) return;
        h.addEventListener('mousedown', function(e) {
          dragging     = true;
          activeHandle = h;
          activePanel  = p;
          direction    = dir;
          startX       = e.clientX;
          startW       = p.getBoundingClientRect().width;
          h.classList.add('dragging');
          e.preventDefault();
        });
      }
      document.addEventListener('DOMContentLoaded', function() {
        initResize('resize-handle',   'lpanel', 1);
        initResize('resize-handle-r', 'rpanel', -1);
      });
      document.addEventListener('mousemove', function(e) {
        if (!dragging || !activePanel) return;
        var delta = (e.clientX - startX) * direction;
        var newW  = Math.min(Math.max(startW + delta, 180), 700);
        activePanel.style.width    = newW + 'px';
        activePanel.style.minWidth = newW + 'px';
      });
      document.addEventListener('mouseup', function() {
        dragging = false;
        if (activeHandle) activeHandle.classList.remove('dragging');
        activeHandle = null;
        activePanel  = null;
      });
    })();
    Shiny.addCustomMessageHandler('leaflet_resize', function(msg) {
      setTimeout(function() {
        var map = HTMLWidgets.find('#' + msg.id);
        if (map) map.getMap().invalidateSize();
      }, 100);
    });
  "))
)
# ─────────────────────────────────────────────────────────────────────────────
# SERVER
# ─────────────────────────────────────────────────────────────────────────────
server <- function(input, output, session) {
  rv <- reactiveValues(
    dem          = NULL,
    perfil       = NULL,
    los          = NULL,
    click_modo   = FALSE,
    next_pt      = "A",
    dist_m       = NA,
    edificios    = NULL,
    roads        = NULL,
    lon_a_calc   = NULL,
    lat_a_calc   = NULL,
    lon_b_calc   = NULL,
    lat_b_calc   = NULL
  )
  # ── Mapa base ──────────────────────────────────────────────────────────────
  output$mapa <- renderLeaflet({
    leaflet() %>%
      addProviderTiles("Esri.WorldGrayCanvas",  group = "Gris") %>%
      addProviderTiles("CartoDB.DarkMatter",    group = "Oscuro") %>%
      addProviderTiles("Esri.WorldImagery",     group = "Satélite") %>%
      addProviderTiles("CartoDB.Positron",      group = "Claro") %>%
      addProviderTiles("OpenStreetMap.Mapnik",  group = "OSM") %>%
      addLayersControl(
        baseGroups   = c("Gris", "Oscuro", "Satélite", "Claro", "OSM"),
        overlayGroups = c("Edificios OSM", "Vías OSM"),
        options      = layersControlOptions(collapsed = TRUE)
      ) %>%
      setView(lng = -3.703, lat = 40.418, zoom = 13)
  })
  # ── Observe: dibujar mapa + OSM cuando rv$los cambia ─────────────────────
  observeEvent(rv$los, {
    req(rv$los, rv$lon_a_calc, rv$lat_a_calc, rv$lon_b_calc, rv$lat_b_calc)
    los       <- rv$los
    lon_a     <- rv$lon_a_calc
    lat_a     <- rv$lat_a_calc
    lon_b     <- rv$lon_b_calc
    lat_b     <- rv$lat_b_calc
    col_linea <- if (isTRUE(los$visible)) "#b8e830" else "#e83030"
    # ── Línea LOS
    p <- leafletProxy("mapa") %>%
      clearGroup("LOS") %>%
      clearGroup("Obstruccion") %>%
      clearGroup("Edificios OSM") %>%
      clearGroup("Vías OSM") %>%
      addPolylines(
        lng = c(lon_a, lon_b), lat = c(lat_a, lat_b),
        group = "LOS", color = col_linea, weight = 2.5, opacity = .85,
        dashArray = if (isTRUE(los$visible)) NULL else "7,5"
      )
    # ── Marcador obstrucción
    if (!isTRUE(los$visible) && !is.na(los$bloqueo_dist)) {
      frac    <- los$bloqueo_dist / los$dist_total
      lon_blq <- lon_a + frac * (lon_b - lon_a)
      lat_blq <- lat_a + frac * (lat_b - lat_a)
      p <- p %>% addCircleMarkers(
        lng = lon_blq, lat = lat_blq, group = "Obstruccion",
        radius = 8, color = "#e83030", fillColor = "#e83030",
        fillOpacity = 0.9, weight = 2, opacity = 1,
        label = paste0("⚠ Obstrucción · ", round(los$bloqueo_dist, 0),
                       " m · cota ", round(los$bloqueo_elev, 0), " m"),
        labelOptions = labelOptions(
          permanent = TRUE, direction = "top", offset = c(0, -12),
          style = list(
            "background-color" = "#1e2023", "color" = "#e83030",
            "font-family" = "IBM Plex Mono, monospace", "font-size" = "10px",
            "border-color" = "#e83030", "font-weight" = "600", "padding" = "3px 6px"
          )
        )
      )
    }
    # ── Zoom
    pad <- 0.004
    p %>% fitBounds(
      lng1 = min(lon_a, lon_b) - pad, lat1 = min(lat_a, lat_b) - pad,
      lng2 = max(lon_a, lon_b) + pad, lat2 = max(lat_a, lat_b) + pad
    )
    # ── Descargar OSM y dibujar en observe separado via rv
    osm <- tryCatch(
      descargar_osm(lon_a, lat_a, lon_b, lat_b, buffer_m = 200),
      error = function(e) list(edificios = NULL, roads = NULL)
    )
    rv$edificios <- osm$edificios
    rv$roads     <- osm$roads
  })
  # ── Observe: dibujar capas OSM cuando rv$edificios cambia ─────────────────
  observeEvent(rv$edificios, {
    req(rv$lon_a_calc)
    p <- leafletProxy("mapa") %>%
      clearGroup("Edificios OSM") %>%
      clearGroup("Vías OSM")
    if (!is.null(rv$edificios) && nrow(rv$edificios) > 0) {
      p <- p %>% addPolygons(
        data = rv$edificios, group = "Edificios OSM",
        color = "#f0a030", weight = 1, opacity = .9,
        fillColor = "#f0a030", fillOpacity = .25,
        label = ~paste0("id: ", osm_id),
        labelOptions = labelOptions(style = list(
          "background-color" = "#1e2023", "color" = "#f0a030",
          "font-family" = "IBM Plex Mono", "font-size" = "10px"
        ))
      )
    }
    if (!is.null(rv$roads) && nrow(rv$roads) > 0) {
      p <- p %>% addPolylines(
        data = rv$roads, group = "Vías OSM",
        color = "#a0b8d0", weight = 1.5, opacity = .7,
        label = ~name
      )
    }
  })
  # ── Toggle modo clic ───────────────────────────────────────────────────────
  observeEvent(input$toggle_clic, {
    rv$click_modo <- !rv$click_modo
    rv$next_pt    <- "A"
  })
  output$click_estado <- renderText({
    if (rv$click_modo)
      paste0("Clic activo — define Punto ", rv$next_pt, " en el mapa")
    else ""
  })
  # ── Capturar clic en mapa ──────────────────────────────────────────────────
  observeEvent(input$mapa_click, {
    if (!rv$click_modo) return()
    cl <- input$mapa_click
    if (rv$next_pt == "A") {
      updateNumericInput(session, "lat_a", value = round(cl$lat, 5))
      updateNumericInput(session, "lon_a", value = round(cl$lng, 5))
      rv$next_pt <- "B"
    } else {
      updateNumericInput(session, "lat_b", value = round(cl$lat, 5))
      updateNumericInput(session, "lon_b", value = round(cl$lng, 5))
      rv$click_modo <- FALSE
      rv$next_pt    <- "A"
    }
  })
  # ── Actualizar marcadores reactivamente ───────────────────────────────────
  observe({
    req(input$lat_a, input$lon_a, input$lat_b, input$lon_b)
    leafletProxy("mapa") %>%
      clearGroup("Puntos") %>%
      addCircleMarkers(
        lng = input$lon_a, lat = input$lat_a, group = "Puntos",
        radius = 7, color = "#b8e830", fillColor = "#b8e830", fillOpacity = .9,
        weight = 2, opacity = 1, label = "A — Observador"
      ) %>%
      addCircleMarkers(
        lng = input$lon_b, lat = input$lat_b, group = "Puntos",
        radius = 7, color = "#30b8e8", fillColor = "#30b8e8", fillOpacity = .9,
        weight = 2, opacity = 1, label = "B — Objetivo"
      )
  })
  # ── Calcular LOS ───────────────────────────────────────────────────────────
  observeEvent(input$calcular, {
    req(input$lat_a, input$lon_a, input$lat_b, input$lon_b)
    # Reset
    rv$dem    <- NULL
    rv$perfil <- NULL
    rv$los    <- NULL
    withProgress(message = "Descargando MDT05 del IGN…", value = .1, {
      tryCatch({
        lon_min <- min(input$lon_a, input$lon_b)
        lon_max <- max(input$lon_a, input$lon_b)
        lat_min <- min(input$lat_a, input$lat_b)
        lat_max <- max(input$lat_a, input$lat_b)
        dem <- descargar_mdt05(lon_min, lat_min, lon_max, lat_max,
                               buffer_deg = input$buf)
        rv$dem <- dem
        incProgress(.4, message = "Extrayendo perfil de elevación…")
        perfil <- perfil_elevacion(dem,
                                   input$lon_a, input$lat_a,
                                   input$lon_b, input$lat_b,
                                   n = as.integer(input$n_pts))
        rv$perfil  <- perfil
        rv$dist_m  <- perfil$distance[nrow(perfil)]
        incProgress(.3, message = "Calculando LOS…")
        los <- calcular_los(perfil, input$h_a, input$h_b)
        rv$los <- los
        incProgress(.05, message = "Listo.")
        # Guardar coords para el observe de mapa
        rv$lon_a_calc <- input$lon_a
        rv$lat_a_calc <- input$lat_a
        rv$lon_b_calc <- input$lon_b
        rv$lat_b_calc <- input$lat_b
      }, error = function(e) {
        showNotification(paste("Error:", conditionMessage(e)),
                         type = "error", duration = 10)
      })
    })
  })
  # ── Reset ─────────────────────────────────────────────────────────────────
  observeEvent(input$resetear, {
    rv$dem       <- NULL
    rv$perfil    <- NULL
    rv$los       <- NULL
    rv$edificios <- NULL
    rv$dist_m    <- NA
    leafletProxy("mapa") %>%
      clearShapes() %>%
      clearMarkers() %>%
      clearGroup("Edificios OSM") %>%
      clearGroup("Vías OSM") %>%
      setView(lng = -3.703, lat = 40.418, zoom = 13)
  })
  # ── Exportar GPKG ────────────────────────────────────────────────────────
  output$exportar_gpkg <- downloadHandler(
    filename = function() {
      paste0("LOS_Madrid_", format(Sys.time(), "%Y%m%d_%H%M%S"), ".gpkg")
    },
    content = function(file) {
      req(rv$perfil, rv$los)
      # ── Punto A
      pt_a <- st_sf(
        id    = "A",
        tipo  = "Observador",
        cota  = rv$los$elev_a,
        geometry = st_sfc(st_point(c(isolate(input$lon_a), isolate(input$lat_a))), crs = 4326)
      )
      # ── Punto B
      pt_b <- st_sf(
        id    = "B",
        tipo  = "Objetivo",
        cota  = rv$los$elev_b,
        geometry = st_sfc(st_point(c(isolate(input$lon_b), isolate(input$lat_b))), crs = 4326)
      )
      # ── Línea LOS
      linea_los <- st_sf(
        visible      = rv$los$visible,
        dist_m       = round(rv$dist_m, 1),
        angulo_deg   = round(rv$los$angle_deg, 3),
        exceso_m     = round(rv$los$exceso_max, 2),
        bloqueo_dist = ifelse(is.na(rv$los$bloqueo_dist), NA_real_,
                              round(rv$los$bloqueo_dist, 1)),
        geometry = st_sfc(
          st_linestring(rbind(
            c(isolate(input$lon_a), isolate(input$lat_a)),
            c(isolate(input$lon_b), isolate(input$lat_b))
          )),
          crs = 4326
        )
      )
      # ── Perfil de elevación como puntos con estado LOS por punto
      perfil_df <- rv$perfil[, c("idx","distance","elevation","lon","lat")]
      # Calcular elevación de la LOS en cada punto del perfil
      elev_a_obs <- rv$los$elev_a + isolate(input$h_a)
      elev_b_obs <- rv$los$elev_b + isolate(input$h_b)
      slope_los  <- (elev_b_obs - elev_a_obs) / rv$los$dist_total
      elev_los_pts <- elev_a_obs + slope_los * perfil_df$distance
      # Punto bloqueado = terreno supera la LOS
      perfil_df$elev_LOS   <- round(elev_los_pts, 2)
      perfil_df$exceso_m   <- round(perfil_df$elevation - elev_los_pts, 2)
      perfil_df$bloqueado  <- perfil_df$elevation > elev_los_pts
      perfil_df$estado     <- ifelse(perfil_df$bloqueado, "BLOQUEADO", "VISIBLE")
      perfil_sf <- st_as_sf(perfil_df, coords = c("lon","lat"), crs = 4326)
      # ── Escribir GPKG — capas vectoriales
      st_write(pt_a,      file, layer = "punto_A",          driver = "GPKG", quiet = TRUE)
      st_write(pt_b,      file, layer = "punto_B",          driver = "GPKG", append = TRUE, quiet = TRUE)
      st_write(linea_los, file, layer = "linea_LOS",        driver = "GPKG", append = TRUE, quiet = TRUE)
      st_write(perfil_sf, file, layer = "perfil_elevacion", driver = "GPKG", append = TRUE, quiet = TRUE)
      # ── Edificios OSM
      if (!is.null(rv$edificios) && nrow(rv$edificios) > 0) {
        # Quedarse solo con columnas simples (sin geometrías anidadas de st_intersection)
        edif_export <- tryCatch({
          edif <- rv$edificios
          edif <- edif[, c("osm_id", "name", "geometry")]
          edif$osm_id <- as.character(edif$osm_id)
          # Forzar solo polígonos/multipolígonos
          edif <- edif[st_geometry_type(edif) %in% c("POLYGON","MULTIPOLYGON"), ]
          if (nrow(edif) > 0) edif else NULL
        }, error = function(e) NULL)
        if (!is.null(edif_export) && nrow(edif_export) > 0) {
          st_write(edif_export, file, layer = "edificios_OSM",
                   driver = "GPKG", append = TRUE, quiet = TRUE)
        }
      }
      # ── Vías OSM
      if (!is.null(rv$edificios)) {
        # roads están en rv implícitamente — los regeneramos del osm guardado
        # Para no re-descargar, guardamos roads en rv también
        if (!is.null(rv$roads) && nrow(rv$roads) > 0) {
          roads_export <- tryCatch({
            r <- rv$roads[, c("osm_id","name","highway","geometry")]
            r$osm_id <- as.character(r$osm_id)
            r[st_geometry_type(r) %in% c("LINESTRING","MULTILINESTRING"), ]
          }, error = function(e) NULL)
          if (!is.null(roads_export) && nrow(roads_export) > 0) {
            st_write(roads_export, file, layer = "vias_OSM",
                     driver = "GPKG", append = TRUE, quiet = TRUE)
          }
        }
      }
      # ── Punto de obstrucción (si existe)
      if (!isTRUE(rv$los$visible) && !is.na(rv$los$bloqueo_dist)) {
        frac    <- rv$los$bloqueo_dist / rv$los$dist_total
        lon_blq <- isolate(input$lon_a) + frac * (isolate(input$lon_b) - isolate(input$lon_a))
        lat_blq <- isolate(input$lat_a) + frac * (isolate(input$lat_b) - isolate(input$lat_a))
        pt_blq  <- st_sf(
          tipo         = "Obstrucción",
          dist_desde_A = round(rv$los$bloqueo_dist, 1),
          cota         = round(rv$los$bloqueo_elev, 1),
          exceso_m     = round(rv$los$exceso_max, 2),
          geometry     = st_sfc(st_point(c(lon_blq, lat_blq)), crs = 4326)
        )
        st_write(pt_blq, file, layer = "obstruccion",
                 driver = "GPKG", append = TRUE, quiet = TRUE)
      }
    }
  )
  # ── UI Resultado ──────────────────────────────────────────────────────────
  output$resultado_ui <- renderUI({
    if (is.null(rv$los)) {
      return(tagList(
        div(class = "res-status",
          div(class = "dot na"),
          span(class = "res-txt na", "Sin calcular")
        )
      ))
    }
    los     <- rv$los
    visible <- isTRUE(los$visible)
    cls     <- if (isTRUE(visible)) "ok" else "ko"
    lbl     <- if (isTRUE(visible)) "VISIBLE" else "BLOQUEADO"
    tagList(
      div(class = "res-status",
        div(class = paste("dot", cls)),
        span(class = paste("res-txt", cls), lbl)
      ),
      div(class = "kv",
        span(class = "k", "Distancia 2D"),
        span(class = "v", paste0(round(rv$dist_m, 0), " m"))
      ),
      div(class = "kv",
        span(class = "k", "Cota Punto A"),
        span(class = "v", paste0(round(los$elev_a, 1), " m"))
      ),
      div(class = "kv",
        span(class = "k", "Cota Punto B"),
        span(class = "v", paste0(round(los$elev_b, 1), " m"))
      ),
      div(class = "kv",
        span(class = "k", "Ángulo elevación"),
        span(class = "v", paste0(round(los$angle_deg, 2), "°"))
      ),
      if (!isTRUE(visible)) div(class = "kv",
        span(class = "k", "Obstrucción a"),
        span(class = "v", paste0(round(los$bloqueo_dist, 0), " m"))
      ),
      if (!isTRUE(visible)) div(class = "kv",
        span(class = "k", "Exceso sobre LOS"),
        span(class = "v", paste0(round(los$exceso_max, 1), " m"))
      ),
      div(class = "kv",
        span(class = "k", "Edificios ±200 m"),
        span(class = "v",
          if (is.null(rv$edificios)) "—"
          else paste0(nrow(rv$edificios), " edificios")
        )
      )
    )
  })
  # ── Plot perfil ───────────────────────────────────────────────────────────
  output$perfil_plot <- renderPlot({
    req(rv$perfil, rv$los)
    df  <- rv$perfil
    los <- rv$los
    n   <- nrow(df)
    col_los <- if (isTRUE(los$visible)) "#b8e830" else "#e83030"
    par(
      bg = "#252830", col.axis = "#444", col.lab = "#555",
      fg = "#1e1e1e", mar = c(3.5, 4, 1, 1), family = "mono"
    )
    y_min <- min(df$elevation, na.rm = TRUE) - 5
    y_max <- max(max(df$elevation, na.rm = TRUE),
                 max(los$elev_los, na.rm = TRUE)) + 10
    plot(df$distance, df$elevation,
         type = "n", ylim = c(y_min, y_max),
         xlab = "Distancia (m)", ylab = "Elevación (m s.n.m.)",
         las = 1, cex.axis = .72, cex.lab = .8,
         col.lab = "#555", col.axis = "#444")
    # 1. Relleno terreno (base)
    polygon(c(df$distance, rev(df$distance)),
            c(df$elevation, rep(y_min, n)),
            col = "#1c2a12", border = NA)
    # 2. Zona bloqueada: solo donde terreno > LOS, encima del terreno
    if (!isTRUE(los$visible)) {
      elev_clipped <- pmin(df$elevation, los$elev_los)
      # Polígono entre LOS y terreno donde terreno > LOS
      over <- df$elevation > los$elev_los
      if (any(over, na.rm = TRUE)) {
        for (i in seq_len(n - 1)) {
          if (isTRUE(over[i]) || isTRUE(over[i+1])) {
            xs <- c(df$distance[i], df$distance[i+1],
                    df$distance[i+1], df$distance[i])
            ys <- c(los$elev_los[i], los$elev_los[i+1],
                    df$elevation[i+1], df$elevation[i])
            polygon(xs, ys, col = "#5a1010cc", border = NA)
          }
        }
      }
      # Línea vertical en punto de primer bloqueo
      if (!is.na(los$bloqueo_dist)) {
        abline(v = los$bloqueo_dist, col = "#e83030aa", lty = 3, lwd = 1.2)
      }
    }
    # 3. Línea terreno (encima del relleno)
    lines(df$distance, df$elevation, col = "#4a7020", lwd = 1.8)
    # 4. Línea de visión directa
    lines(df$distance, los$elev_los,
          col = col_los, lwd = 2,
          lty = if (isTRUE(los$visible)) 1 else 2)
    # Puntos A y B
    points(c(0, df$distance[n]),
           c(df$elevation[1] + input$h_a, df$elevation[n] + input$h_b),
           col = c("#b8e830", "#30b8e8"), pch = 19, cex = 1.3)
    legend("topright",
           legend = c("Terreno", "LOS"),
           col    = c("#4a7020", col_los),
           lwd    = c(1.8, 2),
           lty    = c(1, if (isTRUE(los$visible)) 1 else 2),
           bty = "n", text.col = "#555", cex = .73)
  }, bg = "#252830")
  # ── Stats perfil ──────────────────────────────────────────────────────────
  output$perfil_stats_ui <- renderUI({
    req(rv$perfil)
    df <- rv$perfil
    tagList(
      tags$br(),
      div(class = "kv",
        span(class = "k", "Δ Elevación A→B"),
        span(class = "v", paste0(round(df$elevation[nrow(df)] - df$elevation[1], 1), " m"))
      ),
      div(class = "kv",
        span(class = "k", "Máx. perfil"),
        span(class = "v", paste0(round(max(df$elevation, na.rm=TRUE), 0), " m"))
      ),
      div(class = "kv",
        span(class = "k", "Mín. perfil"),
        span(class = "v", paste0(round(min(df$elevation, na.rm=TRUE), 0), " m"))
      )
    )
  })
}
# ─────────────────────────────────────────────────────────────────────────────
shinyApp(ui = ui, server = server)

Alberto C.
GeoAI analyst

https://posit.co/download/rstudio-desktop
https://mapterhorn.com/
https://www.bluemarblegeo.com/global-mapper/
https://www.tekonelectronics.com/en/news/tekon-blog/line-of-sight-real-range-wireless-communications/

Leave a comment