POI Intelligence for Urban Asset Analysis: assetIQ

When analysing urban assets, there is genuine value in moving beyond generic neighborhood scores. The density of a coffee shop cluster, the proximity to a financial hub, or the concentration of accommodation around a transport node are signals that traditional datasets flatten into averages — or ignore entirely. assetIQ was built to change that.

Diagram 1 – Turning the mess into value by using assetIQ

What is assetIQ?

assetIQ is an R application powered by DuckDB and Overture Maps that extracts, classifies, and scores Points of Interest (POIs) for any location on Earth. You define a city and a search radius — from 100 meters to 25 kilometers — and the tool queries the Overture Maps Places dataset in real time, classifying each POI into thematic groups: Food & Drink, Retail, Health, Education, Transport, Accommodation, Financial Services, Leisure & Culture, Sport, and more.

Diagram 2 – Retail POIs in Madrid downtown

The core output is an attribute value called POIQ: a normalized 0–1 score assigned to every building footprint within the area of interest, derived from a Kernel Density Estimation of the selected thematic group. A building in a dense retail corridor scores close to 1. An isolated residential block far from any commerce scores close to 0. This transforms thousands of individual points — which in raw form tell you very little — into a single, interpretable attribute per building, ready for downstream modelling, valuation, or site selection.

Diagram 3 – Accommodation POIs in Madrid downtown

A companion Median Center marker identifies the geographic centroid of maximum concentration for the selected group, giving a precise, reproducible anchor for the activity zone rather than a subjective description.

Diagram 4 – Median center of the POI distribution

Who is it for?

The most immediate use cases sit in industries where location is a core value driver:

Real Estate & Property Valuation — retail proximity, hospitality density, and financial services concentration are established drivers of commercial and residential yield. POIQ provides a quantified, reproducible variable to include directly in hedonic pricing models or investment scoring frameworks.

Diagram 5 – The output: POIQ value in the attribute table

Telecommunications & Network Planning — coverage prioritization, site acquisition for small cells or retail outlets, and churn modelling all benefit from understanding which buildings sit inside high-activity commercial ecosystems versus which are isolated. POIQ adds a demand-side spatial signal to network infrastructure decisions.

  • Carrousel A

Retail & Franchise Expansion — identifying whether a candidate site is surrounded by complementary commerce (a food cluster that attracts footfall) or by competing category saturation is exactly the kind of micro-level distinction assetIQ makes visible.

Urban Planning & Consultancy — tracking how thematic POI density shifts across a corridor over time, or comparing two neighborhoods before and after an infrastructure intervention, becomes a structured, repeatable workflow.

Insurance & Risk — the presence of hospitality, nightlife, or financial services clusters correlates with foot traffic, vandalism exposure, and commercial risk profiles. POIQ can serve as a spatial covariate in risk scoring models.

Data sources and technology

Overture Maps Foundation is a collaborative open-source mapping initiative founded by AWS, Meta, Microsoft, and TomTom under the Linux Foundation, combining brand-verified location data with open contributions to produce a globally consistent POI dataset. Recent releases have incorporated data from Foursquare Open Source Places, adding millions of new POIs to expand coverage, with each record licensed under either Apache 2.0 or CDLA 2.0 depending on source (building footprints come from the Overture Buildings theme).

  • Carrousel B

The application is written entirely in R, using Shiny for the interactive interface, DuckDB with the httpfs extension for serverless remote querying of GeoParquet files directly from S3, sf for spatial operations, and MASS for kernel density estimation. No data is downloaded to disk during analysis — DuckDB streams only the rows that fall within the bounding box of interest, making even large-radius queries tractable in under two minutes.

Outputs are exported as GeoPackage (.gpkg) — an OGC-standard format readable by QGIS, ArcGIS, and any GIS-capable pipeline — with separate layers for POIs (full attribute table including category taxonomy, confidence score, brand, contact, and address fields), building footprints with the POIQ field in the attribute table, and the Median Center point. GeoJSON export is straightforward from any of these layers for web or API integration.

Limitations

The quality of any POIQ score depends directly on the completeness and classification accuracy of the underlying POI data. Independent analysis of the Overture Places dataset has found that while coverage of major branded locations is strong in high-income markets, smaller independent businesses and coverage in less-mapped regions can be uneven — with some brand-country combinations showing coverage ratios well below the ideal threshold. The confidence field provided by Overture (0–1, representing the estimated probability that a place actually exists at the reported location) is preserved in the output and can be used to filter or weight results downstream.

  • Carrousel C

The thematic classification layer in assetIQ — mapping Overture’s basic_category taxonomy to eleven groups — is a deliberate simplification. Edge cases exist: a hospital pharmacy classified under Health rather than Retail, a hotel gym that could sit in either Accommodation or Sport. These classification boundaries should be treated as configurable starting points rather than ground truth.

Diagram 6 – Disaggregation of accomodation Madrid downtown

Finally, POIQ is a relative index within a given area of interest, not an absolute score. A POIQ of 0.8 in a 500m radius around Puerta del Sol is not directly comparable to a POIQ of 0.8 in a 5km radius around a suburban retail park — the denominator changes with the search parameters.

Interested?

assetIQ is built for analysts, developers, and decision-makers who want to move from “there are lots of restaurants nearby” to “this building sits at the 94th percentile of Food & Drink density within its competitive set.” Ping me if you want to know more. Download this sample, simbolize it on your own and let me know if it suits you!

Built with R · DuckDB · Overture Maps

Diagram 7 – Raw output, what you have just downloaded. Now imagine it over your own AOI
sample_assetIQ_20260528Descarga

A Geospatial solution for everybody!

Alberto C.
GIS analyst

https://overturemaps.org
https://www.openstreetmap.org/
https://posit.co/download/rstudio-desktop
https://www.bluemarblegeo.com/global-mapper/

Mapterhorn + R: LOS analysis in seconds! 🚀

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.

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/

SWEEPLY! HOME TASKS, SHARED

🏠 “I feel like I do everything around here…”

We’ve all been there. Domestic load is real, and tracking who does what at home can be a challenge—whether you want to balance tasks with your partner, get your kids to collaborate, or simply prove that chores don’t just “do themselves.”

That is why I created SWEEPLY! 🧹 Our motto says it all: HOME TASKS, SHARED.
https://sweeply-home-tasks-shared-v3.netlify.app/

With SWEEPLY, you can:

  • 🗺️ Visualize your home: Generate and use your own floor plan to see tasks room by room.
  • 📊 Quantify effort: Assign point values to chores based on difficulty and frequency so everyone’s contribution is clear.
  • 🚨 Never forget a chore: Set up smart alerts for pending tasks or upcoming deadlines.
  • 📈 Track trends: Monitor charts and weekly insights to see who is leading the scoreboard (shoutout to Paige in my video demo!).

It’s time to turn chores into a fair, visual, and engaging family effort.

I’d love to hear your thoughts! How do you handle the chore split at home?👇

#AppDevelopment #AI #Productivity #GIS #FamilyLife #TechForGood #Sweeply

As a GIS analyst, I spend my days working with XYZ coordinates and mapping spaces. Thanks to AI technology, I’ve finally been able to bring a long-time idea to life and apply that spatial logic to home organization! So hopefully you like this off-topic stuff 🙂

Alberto C.
GIS analyst and new developer of funny stuff 🙂

Detecting Potential Mobile Coverage Gaps Using OpenCellID, GHSL and Overture Maps: Case study over TUNIS

Mobile connectivity has become a fundamental component of modern infrastructure, yet significant spatial inequalities in network access still persist across both urban peripheries and rural environments. Using openly available geospatial datasets, this analysis explores potential mobile coverage gaps by combining OpenCellID cellular infrastructure observations, GHSL population layers and vector data extracted from Overture Maps. The objective is not to reproduce real telecom propagation models, but to generate a simplified spatial estimation of coverage capable of identifying populated areas potentially located outside the influence of nearby cellular infrastructure.

Diagram 1 – In Global Mapper I buffered by RANGE all possible technologies available over Tunis

OpenCellID provides crowdsourced observations of mobile network cells, including location, radio technology and an estimated operational range. The dataset includes GSM (2G), UMTS (3G), LTE (4G) and, in some regions, NR (5G) deployments. In practical terms, the reported range represents an approximate distance over which a cell may provide service under average conditions, although real coverage depends on many additional factors such as topography, antenna orientation, carrier frequency, urban density and atmospheric conditions. In this workflow, the “range” attribute was used to generate simplified coverage buffers around each observed cell, allowing the construction of continuous spatial representations of network availability by technology.

Diagram 2 – After applying my methology in Tunis, in yellow areas >100 inhabs. not covered

Population distribution was derived from the Global Human Settlement Layer (GHSL), developed by the Joint Research Centre of the European Commission. GHSL provides spatially explicit population estimates derived from satellite imagery and census modeling, making it possible to compare settlement patterns directly against estimated cellular coverage.

Overture Maps vector layers, including buildings, roads and hydrographic features, were integrated to provide geographic context and improve interpretation of the spatial structure of the analyzed regions. The workflow was implemented using javascript in GEE and R, where cellular infrastructure datasets were imported as custom assets and processed interactively.

Diagram 3 – After applying my methology in Tunis, in yellow areas >100 inhabs. not covered (zoom in)

Two study areas were selected to illustrate different territorial and infrastructural conditions: Tunis and Luxembourg. Tunis represents a dense and rapidly expanding metropolitan environment with heterogeneous infrastructure distribution, while Luxembourg provides an example of a smaller and highly connected European territory with dense transport and telecommunications networks. Preliminary results show strong LTE and UMTS concentration around urban cores and major transport corridors, whereas fragmented or sparse coverage appears more frequently in peripheral settlements and lower-density areas. Several populated clusters identified through GHSL remain outside the simplified coverage envelopes, highlighting potential candidates for further connectivity assessment.

Diagram 4 – After applying my methodology in Tunis, in yellow areas >100 inhabs. not covered (more zoom in!!)

This analysis contains important limitations. OpenCellID is a crowdsourced dataset with variable completeness across regions, and the reported cell ranges are only approximate estimates rather than engineering measurements. The simplified buffering methodology does not account for terrain effects, signal attenuation, antenna directionality or radio interference, while GHSL population values are themselves modeled estimates. Consequently, the results should be interpreted as exploratory indicators of possible connectivity gaps rather than authoritative coverage maps. Nevertheless, the integration of open telecom observations, global population datasets and large-scale vector mapping frameworks demonstrates the growing potential of open geospatial ecosystems for large-scale infrastructure and digital accessibility analysis.

The final output of the workflow is a vector geometry representing areas where population and built-up structures exist but where no nearby cellular coverage has been identified from the technologies available in OpenCellID. These geometries can be exported and directly visualized in platforms such as Google Earth, where they highlight zones that potentially should have mobile coverage but appear to lack it according to the observed GSM, UMTS, LTE and NR datasets used in this analysis.

Diagram 5 – Exporting to google Earth to share results internally (and externally!)

As I had introduced before, additional contextual information from Overture Maps can also be integrated into the same output, including addresses, points of interest, buildings, roads, railways, hydrography, land use and other infrastructure layers, providing a richer spatial framework to interpret the detected gaps. This entire data extraction and integration process can be reproduced in minutes using the R-based application I developed last week, which automates the retrieval and preprocessing of Overture Maps vector datasets for geospatial analysis workflows.

From Overture Maps to GPKG in minutes: Building a Geospatial Data Extractor with R and DuckDB
Diagram 6 – Enhancing thr model by extracing Overture Maps using the application I built last week using R

Now we can automatize the whole procedure switching to a different ASSET and that’d take just a few minutes to get to know if we have potential coverage issues we might address right away!. This below has been adapted to do the same over Luxembourg. Only to know these guys are completely covered!!

var gsm  = ee.FeatureCollection('projects/analysis202601/assets/dissolved_GSM_lux');
var lte  = ee.FeatureCollection('projects/analysis202601/assets/dissolved_LTE_lux');
var umts = ee.FeatureCollection('projects/analysis202601/assets/dissolved_UMTS_lux');
Map.centerObject(gsm, 9);
Map.setOptions('HYBRID');
// ── SRTM hillshade ────────────────────────────────────────────
var srtm = ee.Image('USGS/SRTMGL1_003');
var hillshade = ee.Terrain.hillshade(srtm, 315, 45);
Map.addLayer(hillshade, {min: 0, max: 255, opacity: 0.35}, 'Hillshade SRTM', true);
// ── GHSL población ────────────────────────────────────────────
var ghsl = ee.ImageCollection('JRC/GHSL/P2023A/GHS_POP')
  .filter(ee.Filter.date('2020-01-01', '2020-12-31'))
  .select('population_count')
  .mosaic();
var POP_MIN = 0;
var POP_MAX = 1000;
var POP_THRESHOLD = 100;
var popMask = ghsl.gt(POP_THRESHOLD);
var ghslVis = {min: POP_MIN, max: POP_MAX, palette: ['ffffff', 'ffcccc', 'ff8888', 'ff3333', 'cc0000']};
Map.addLayer(ghsl.updateMask(ghsl.gt(0)), ghslVis, 'GHSL población', true);
// ── Coberturas outlined ───────────────────────────────────────
var empty = ee.Image().byte();
Map.addLayer(empty.paint({featureCollection: gsm,  color: 1, width: 2}), {palette: ['e040fb'], opacity: 0.9}, 'GSM outline',  true);
Map.addLayer(empty.paint({featureCollection: lte,  color: 1, width: 2}), {palette: ['00e676'], opacity: 0.9}, 'LTE outline',  true);
Map.addLayer(empty.paint({featureCollection: umts, color: 1, width: 2}), {palette: ['ff3d3d'], opacity: 0.9}, 'UMTS outline', true);
// ── Huecos de cobertura ───────────────────────────────────────
var allCoverage = gsm.merge(lte).merge(umts);
var gapsImg = ee.Image(1).paint(allCoverage, 0).selfMask().updateMask(popMask);
var roi = gsm.geometry().bounds().buffer(50000);
var gapsVec = gapsImg.reduceToVectors({geometry: roi, scale: 100, maxPixels: 1e9});
Map.addLayer(empty.paint({featureCollection: gapsVec, color: 1, width: 3}), {palette: ['ffff00'], opacity: 1}, 'Huecos de cobertura', true);
Export.table.toDrive({
  collection: gapsVec.map(function(f) { return f.setGeometry(f.geometry().simplify(200)); }),
  description: 'huecos_cobertura_luxemburgo',
  fileFormat: 'GeoJSON',
  folder: 'GEE_exports'
});
print('Total huecos:', gapsVec.size());
print('GSM polígonos:',  gsm.size());
print('LTE polígonos:',  lte.size());
print('UMTS polígonos:', umts.size());
// ── Leyenda ───────────────────────────────────────────────────
var palette = ['ffffff', 'ffcccc', 'ff8888', 'ff3333', 'cc0000'];
var colorbar = ui.Thumbnail({
  image: ee.Image.pixelLonLat().select('longitude')
    .multiply((POP_MAX - POP_MIN) / 360.0)
    .add(POP_MAX / 2.0)
    .visualize({min: POP_MIN, max: POP_MAX, palette: palette}),
  params: {bbox: [-180, -5, 180, 5], dimensions: '200x16'},
  style: {stretch: 'horizontal', margin: '4px 0px'}
});
var legendPanel = ui.Panel({
  style: {
    position: 'bottom-right',
    padding: '10px 14px',
    backgroundColor: 'rgba(255,255,255,0.92)',
    width: '230px'
  }
});
legendPanel.add(ui.Label('Población GHSL 2020', {
  fontWeight: 'bold', fontSize: '13px', margin: '0 0 4px 0'
}));
legendPanel.add(ui.Label('habitantes / 100m pixel', {
  fontSize: '10px', color: '#888888', margin: '0 0 6px 0'
}));
legendPanel.add(colorbar);
var labelsRow = ui.Panel({
  layout: ui.Panel.Layout.flow('horizontal'),
  style: {stretch: 'horizontal', margin: '2px 0 0 0'}
});
labelsRow.add(ui.Label(String(POP_MIN), {fontSize: '10px', margin: '0'}));
labelsRow.add(ui.Label('▲ ' + String(POP_THRESHOLD) + ' umbral huecos', {
  fontSize: '10px', color: '#ffcc00', textAlign: 'center', stretch: 'horizontal'
}));
labelsRow.add(ui.Label(POP_MAX + '+', {fontSize: '10px', margin: '0', textAlign: 'right'}));
legendPanel.add(labelsRow);
legendPanel.add(ui.Panel({
  style: {backgroundColor: '#dddddd', height: '1px', margin: '8px 0 6px 0', stretch: 'horizontal'}
}));
legendPanel.add(ui.Label('Cobertura de red', {
  fontWeight: 'bold', fontSize: '13px', margin: '0 0 6px 0'
}));
var techs = [
  {label: 'GSM',             color: '#e040fb'},
  {label: 'LTE',             color: '#00e676'},
  {label: 'UMTS',            color: '#ff3d3d'},
  {label: 'Hueco cobertura', color: '#ffff00'}
];
techs.forEach(function(t) {
  var row = ui.Panel({layout: ui.Panel.Layout.flow('horizontal'), style: {margin: '2px 0'}});
  row.add(ui.Label('━━', {color: t.color, fontSize: '14px', margin: '0 8px 0 0'}));
  row.add(ui.Label(t.label, {fontSize: '11px', margin: '2px 0 0 0'}));
  legendPanel.add(row);
});
Map.add(legendPanel);
Diagram 7 – A happy person making sure no gap of GSM /UMTS /GSM coverture in Luxembourg!

Alberto C.
GeoAI analyst

https://www.opencellid.org
https://human-settlement.emergency.copernicus.eu/
https://earthengine.google.com/
https://www.google.es/intl/es/earth/index.html

CHANGE DETECTION ARCGIS PRO AND LIVING ATLAS 2017-2025

The quantification of land-use dynamics necessitates a spatiotemporal framework that ensures categorical stability over long-term observation windows. The ESRI 10-Meter Global Land Cover time series, accessible through the ArcGIS Living Atlas, provides a harmonized baseline for this purpose, derived from the dense temporal stack of the ESA Sentinel-2 mission.

Diagram 1 – you can download any region in the world 2017-2025 from Living Atlas

Since 2017, this dataset has utilized a consistent Deep Learning (U-Net) classification architecture, which is critical for mitigating the thematic noise often introduced by changing sensors or shifting methodologies. By maintaining a uniform 10-meter resolution and a standardized 9-class nomenclature, the series allows for a direct pixel-to-pixel comparison across years, enabling the identification of subtle urban expansion and environmental degradation without the requirement for secondary normalization or complex post-processing.

Diagram 2 – First preliminary checks to get to fully understand data

This repository is exposed via ArcGIS Pro as a series of multidimensional Image Services, a format that prioritizes computational efficiency through server-side processing. Rather than downloading massive raster tiles, the service allows for the direct execution of Pixel-Level Change Detection algorithms against the cloud-hosted imagery.

Diagram 3 – The actual ArcGIS Pro interface

This “Open Data” approach facilitates a seamless transition from raw spectral data to actionable socio-spatial indicators. For researchers focusing on long-term urban trajectories, the continuity of this data from 2017 to the present ensures that the detected transitions reflect genuine terrestrial shifts rather than artifacts of the classification process. In an ArcGIS Pro environment, this translates to the ability to run Change Detection Raster Functions on-the-fly, producing transition matrices that quantify the conversion of permeable surfaces into built environments with high geographical precision and statistical rigor.

  • 2017
  • 2025
  • CHANGE DETECTION
Video – MP4 from the website livingatlas.arcgis.com/landcoverexplorer/

Analytical Implementation: Categorical Transition Matrices

The core of the change detection workflow in ArcGIS Pro involves the application of the Compute Change Raster Function directly against the multidimensional Image Service slices. Unlike simple pixel subtraction used in spectral indices, categorical change detection requires the generation of a transition matrix to identify specific class-to-class trajectories. By defining the 2017 slice as the initial state and the current year as the final state, the software computes a thematic output where each pixel value represents a unique transition (e.g., Value 5 to Value 7, indicating the conversion of Crops to Built Area).

This process is optimized through the use of Raster Function Templates (RFTs), which allow for server-side processing of the Sentinel-2 derived layers. This avoids the latency of local data managed in geodatabases while maintaining the 10-meter spatial fidelity. For urban vulnerability research, this allows for the isolation of specific “vulnerability triggers”—such as the rapid expansion of impervious surfaces (Built Area) in areas previously classified as Flooded Vegetation (Class 4) or Water (Class 1). The resulting change raster is not merely a visual aid but a spatially explicit database capable of producing area-weighted statistics and frequency distributions that characterize the rate of urban transformation over the last decade.

The integration of the ESRI 10-Meter Annual Land Cover series within ArcGIS Pro establishes a high-fidelity analytical pipeline for quantifying terrestrial transitions. By utilizing that consistent Deep Learning (U-Net) classification architecture across the Sentinel-2 archive from 2017 to the present, the dataset eliminates the methodological variances that typically compromise multi-temporal change detection. This longitudinal consistency is vital for isolating genuine socio-spatial shifts—such as the conversion of Rangeland (Class 11) or Crops (Class 5) into Built Area (Class 7)—from sensor-induced artifacts or classification drift.

Diagram 4 – Change analysis step by step

Accessible as Open Data through the ArcGIS Living Atlas, these layers are delivered as multidimensional Image Services. This cloud-native architecture facilitates direct pixel-level interrogation without the prerequisite for local data ingestion or atmospheric correction. In a technical workflow, these services support the execution of Change Detection Raster Functions, enabling the generation of categorical transition matrices that intersect the 10-meter spatial resolution with the temporal depth of the Sentinel-2 mission. The result is a robust, replicable framework for monitoring urban sprawl, environmental vulnerability, and land-use governance with unprecedented statistical rigor.

Here is the full step-by-step procedure:

Step 1 — Load the rasters

  1. Open ArcGIS Pro and create a new project or open an existing one
  2. In the Catalog panel, navigate to your files
  3. Drag and drop both files into the map:
    • 30T_20170101-20180101_noclouds
    • 30T_20250101-20251231_noclouds

Step 2 — Project to a metric CRS

  1. Go to Analysis → Tools → search for Project Raster
  2. Run it twice (once per raster):
    • Input Raster: 30T_20170101-20180101_noclouds
    • Output Raster: 30T_2017_projected.tif
    • Output Coordinate System: select your UTM zone (e.g. WGS 1984 UTM Zone 30N)
    • Resampling: Nearest (mandatory for classified rasters — never use Bilinear)
  3. Repeat for 30T_20250101-20251231_noclouds → output 30T_2025_projected.tif

Step 3 — Run Compute Change Raster

  1. Go to Analysis → Tools → search for Compute Change Raster
  2. Fill in the parameters:
ParameterValue
Input From Raster30T_2017_projected.tif
Input To Raster30T_2025_projected.tif
Compute Change MethodCategorical Difference
Filter MethodChanged pixels only
From Classes(leave empty)
To Classes7
Output Rasterchange_urban_2017_2025.tif
  1. Click Run

Step 4 — Symbolize the change layer

  1. In the Contents panel, right-click change_urban_2017_2025.tif → Symbology
  2. Set Primary symbology to Unique Values
  3. Select all values in the color table → Ctrl+A
  4. Right-click → Properties for all selected
  5. Set color to #E82020 (or your preferred red)
  6. Set outline to No Color
  7. Close Symbology

Step 5 — Add hillshade background

  1. In the Contents panel, right-click the map → Add Data
  2. Add a basemap: Map → Basemap → Light Gray Canvas (or any hillshade basemap)
  3. In the Contents panel, drag change_urban_2017_2025.tif above the basemap
  4. Right-click the change layer → Properties → Appearance
  5. Set Transparency to 20%
  • 2017
  • 2025
  • any class to Built Area
  • Built Area highlighted
  • Built Area highlighted zoom in
  • Built Area highlighted zoom in and labels

Geographic Information Systems have transformed the way we understand territorial change, turning complex spatial phenomena into clear, measurable evidence. Tools like ArcGIS Pro consolidate in a single environment what once required multiple specialized software packages, reducing analysis time from days to minutes. The change detection workflow presented here — from raw classified rasters to a quantified urban growth map — was completed in under an hour, with no programming required.

This accessibility is key: it allows technicians, planners, and decision-makers to work with the same data and the same conclusions. Urban land monitoring is not a technical exercise in isolation; it feeds directly into land use planning, infrastructure investment, environmental impact assessment, and climate adaptation strategies. Knowing exactly where and how fast built-up area is expanding between two dates is information that governments, developers, and conservation agencies all need. The combination of freely available global products like the Esri Land Cover dataset with the analytical power of ArcGIS Pro and the Sentinel-2 satellite constellation has made this type of analysis accessible at any scale, from a municipality to an entire country. Change detection is not just a methodology — it is a way of making time visible on a map.

Diagram 5 – A potential interactive interface with your own data ser

Interested in these change detection analysis? Let me know!

Alberto
GeoAI Analyst

https://www.arcgis.com/home/item.html?id=cfcb7609de5f478eb7666240902d4d3d

https://livingatlas.arcgis.com/landcoverexplorer/#mapCenter=147.74397%2C-37.18990%2C11.00&mode=step&timeExtent=2017%2C2025&year=2025

Who Gets to See the War? Satellite Imagery Censorship and the Copernicus Alternative

On April 5, 2026, Planet Labs announced it was complying with a US government request to indefinitely withhold high-resolution satellite imagery of a broad swathe of the Middle East — retroactively from March 9 — covering Iran, Iraq, Israel, Lebanon and the Persian Gulf States. The company moved to a “managed access” model, releasing images only on a case-by-case basis.

Diagram 1 – Sentinel 2 before and after Super-Resolution 1m is applied. Mehrabad Airport. 21/4/2026

A private company, operating under a US federal license, was effectively told by the Trump administration to go dark over an active war zone. Planet complied. Vantor and BlackSky — both heavily dependent on US defense revenue — said they hadn’t even received such a request, because they were already tightly controlling access. The selective pressure lands precisely on the most open, most commercially accessible provider. That is not a coincidence.

This is not transparency. This is state-managed information blackout dressed up as voluntary cooperation.

The Asymmetry Is the Message

The US government did not ask all providers equally. It asked Planet — the one with the open-access model, the one serving journalists, researchers, NGOs, and ordinary citizens who want to understand what is happening in their name. Companies deeply embedded in the defense contracting ecosystem didn’t need asking. They were already compliant by design.

Diagram 2 – PLANET vs SENTINEL 2 (STD AND SR)

This asymmetry reveals the real intent: not security, but narrative control. When Bloomberg’s own analysis — using radar imagery from non-US sources — documents over 7,600 damaged or destroyed buildings in Iran including schools and health facilities, and the White House flatly denies it, you understand why high-resolution optical imagery had to disappear.

Enter Copernicus — and Why It Changes Everything

Here is what the US government apparently forgot, or hopes you don’t know: the European Union’s Copernicus programme is not under Washington’s jurisdiction (isn’t it?).

Sentinel-2 provides multispectral imagery at 10-metre resolution. More critically, the Sentinel-1 radar constellation captures the ground through clouds and at night — and it is precisely this data that researchers at Oregon State University used to map the destruction the US government tried to hide.

https://gamayos.github.io/gamma-earth-api/s2dr3-demo-20250305.html?ds=IR-T39SWV-8278809df-20240307#15.7/35.689707/51.298125

Diagram 3 – Mehrabad Airport. Sentinel 2 SR. 07/03/2026 and 21/04/2026

https://gamayos.github.io/gamma-earth-api/s2dr3-demo-20250305.html?ds=IR-T39SWV-8278809df-20260421#15.44/35.691551/51.297791

But the real game-changer is super-resolution. Applying deep learning models to Sentinel-2 data can produce outputs approaching 1-metre effective resolution — sharper than Planet’s native 3–5 m commercial imagery, and entirely outside the reach of any US export control or NOAA licensing mechanism. Copernicus data is free, open, and governed by the European Union, not by an administration that has a direct interest in what gets seen.

A Question of Democratic Accountability

The Copernicus programme was built on a principle: Earth observation is a public good. Citizens, scientists, journalists, and civil society deserve access to independent, verifiable information about what is happening on this planet — including when their governments are the ones doing it.

The United States is now openly demonstrating that it views commercial satellite operators as an extension of its information warfare capability. That should alarm everyone — not just geospatial professionals, but every citizen who believes that a functioning democracy requires an informed public.

The EU has a responsibility here. Copernicus must be defended, funded, and expanded — not as a scientific programme, but as a democratic infrastructure. And the open-source intelligence community, GIS professionals, researchers, and journalists have both the tools and the obligation to use it.

Restricted areas from Planet Labs

Alberto C.
Geospatial analyst

https://newsletter.terrawatchspace.com/the-defense-paradox-in-earth-observation-the-enabler-and-the-hindrance/?ref=terrawatch-space-newsletter

https://www.bloomberg.com/graphics/2026-satellite-imagery-restrictions-iran/

https://www.timesofisrael.com/satellite-firm-says-its-indefinitely-withholding-iran-war-images-at-us-request/

https://www.incibe.es/incibe-cert/blog/osint-la-informacion-es-poder

#Copernicus #RemoteSensing #OpenData #Sentinel #SuperResolution #EarthObservation #Transparency #OSIN

From Overture Maps to GPKG in minutes: Building a Geospatial Data Extractor with R and DuckDB

Why I Developed This Tool

Modern geospatial workflows increasingly depend on fast, reliable access to city-scale vector data — building footprints, road networks, land use polygons, points of interest, address databases. Whether you are designing a 5G radio network, modelling urban heat islands, planning last-mile logistics, or simulating emergency response coverage, you almost always start from the same question: “How do I get clean, structured geodata for this city, right now, without spending two days on it?”

The Overture Maps Extractor is my answer to that question. It is a Shiny application written in R that lets any GIS professional extract multiple thematic layers from the Overture Maps Foundation dataset — for any city in the world — in a matter of minutes, with zero command-line interaction and zero manual data wrangling.

Diagram 1 – Application interface. Final output /Layers available

What Is Overture Maps?

The Overture Maps Foundation is an open data initiative backed by Amazon, Meta, Microsoft and TomTom. It publishes a continuously updated, globally consistent vector dataset covering:

  • Buildings — footprints with height and floor count attributes where available
  • Transportation — roads, railways, and hydrography segments
  • Places (POIs) — points of interest with categories and confidence scores
  • Addresses — street-level address points with number, postcode and city
  • Land Use / LULC — polygonal land cover classification
  • Administrative divisions — municipal and regional boundaries
  • Base land features — parks, forests, water bodies as polygons
  • Infrastructure — energy towers, antennas, utilities

The data is published as GeoParquet files on Amazon S3, updated with every release cycle. The latest release as of this writing is 2026-04-15.0.

The Technical Stack: R + DuckDB + Overture Maps

The core of the extractor is a single R Shiny application. No Python, no Docker, no cloud infrastructure required — just a running R session.

Diagram 2 – Application interface

DuckDB as the Query Engine

The key enabler is DuckDB, an in-process analytical database that can query remote Parquet files over HTTP range requests without downloading the entire dataset. A city-scale extraction that would otherwise require gigabytes of data transfer is reduced to fetching only the relevant row groups.

r

dbExecute(con, "INSTALL httpfs; LOAD httpfs;")
dbExecute(con, "INSTALL spatial; LOAD spatial;")
dbExecute(con, "SET s3_region='us-west-2';")

The spatial filter is pushed directly into the Parquet scan via bounding box predicates on the bbox struct columns — a pattern that Overture’s file layout is specifically optimised for. Overture’s latest schema exposes native GEOMETRY('OGC:CRS84') types, which DuckDB Spatial handles natively via ST_AsText() for clean WKT export to R sf objects.

Diagram 3 – Naming convention and final result of layers’ counts

In this particular case of a spanish scenario over Madrid, we have DTM/DSM Open Data from CNIG, then we can easily get AGL heights from a raster difference (DSM-DTM).

Diagram 4 – Loading CNIG DTM 1m and a DSM extracted from LIDAR 5 samples per sqm (Global Mapper)

Yes, you are right, I did it using global Mapper.

Diagram 5 – adding AGL to the buildings’ footprints I got from Overture Maps (Global Mapper). Gran Via, Madrid

Workflow in Three Steps

1. Locate — type a city name (or drag the interactive marker anywhere on the map) and set a radius in kilometres. The application geocodes via OpenStreetMap Nominatim and computes a circular extraction area using a proper metric buffer in EPSG:3857, avoiding the distortion of a simple bounding box.

2. Select layers — toggle any combination of the ten available thematic layers with checkboxes.

3. Extract — a single button press launches the DuckDB queries, clips the results to the circular area, and writes everything to a single OGC GeoPackage (.gpkg) with standardised layer names: BLDG, ROADS, POIS, ADDRESSES, LULC, LAND, DIVISIONS, INFRA.

The output is immediately ready for QGIS, ArcGIS Pro, PostGIS, or any OGC-compliant GIS tool.

Applications Across Domains

Smart Cities & Urban Analytics Building footprints combined with LULC and POI density form the backbone of urban morphology analysis — walkability indices, urban heat island modelling, green space accessibility, population exposure estimation.

Telecommunications & Radio Planning For 5G and LTE network planning, building height data directly feeds propagation models (ITU-R P.1411, COST-Hata). Address layers enable subscriber geocoding and coverage gap analysis. Road network topology supports site accessibility routing. The circular extraction pattern maps directly to cell sector coverage radii.

Emergency Management & Resilience First-responder route optimisation, shelter location analysis, flood exposure by building footprint — all start from the same data layers this tool produces in minutes.

Real Estate & Infrastructure Investment Rapid assessment of urban fabric density, street connectivity, and amenity proximity for any candidate city worldwide, with a consistent schema regardless of country.

Diagram 6 – Once we import the GPKG into Global Mapper (or any other GIS software)

OGC GeoPackage: The Right Output Format

The output format is deliberately OGC GeoPackage — a single SQLite-based file that is fully OGC-compliant, contains all layers, preserves CRS metadata (EPSG:4326), and requires no proprietary software to open. It is the format of choice for field deployment, data exchange, and reproducible workflows.

Diagram 5 – RStudio interface of the project. As simple as that.

Do you want to code it yourself?. I herewith give you my code:

r

# =============================================================================
#  OVERTURE MAPS EXTRACTOR
#  Dependencies: shiny, bslib, duckdb, DBI, sf, tidygeocoder, leaflet,
#                leaflet.extras, dplyr, glue, shinyjs, waiter, jsonlite
#  Developed by Alberto Concejal / Geovisualization.net 
# =============================================================================
library(shiny); library(bslib); library(duckdb); library(DBI)
library(sf); library(tidygeocoder); library(leaflet); library(leaflet.extras)
library(dplyr); library(glue); library(shinyjs); library(waiter); library(jsonlite)
overture_theme <- bs_theme(
  version=5, bg="#0d1117", fg="#e6edf3", primary="#58a6ff",
  secondary="#3d4451", success="#3fb950", danger="#f85149", warning="#d29922",
  base_font=font_google("IBM Plex Mono"), heading_font=font_google("Space Mono"),
  font_scale=0.9
)
# Capas Overture -> nombre GPKG (railways + waterways se fusionan en ROADS)
LAYER_NAMES <- c(
  buildings="BLDG", roads="ROADS", railways="ROADS", waterways="ROADS",
  pois="POIS", addresses="ADDRESSES", lulc="LULC",
  land="LAND", divisions="DIVISIONS", infrastructure="INFRA"
)
LAYERS <- list(
  buildings = list(
    label="🏛  Edificios", theme="buildings", type="building",
    fields='id, names."primary" AS name, height, num_floors, class, ST_AsText(geometry) AS geom_wkt'
  ),
  roads = list(
    label="🛣  Carreteras", theme="transportation", type="segment",
    fields='id, names."primary" AS name, class, subtype, ST_AsText(geometry) AS geom_wkt',
    filter="AND subtype = 'road'"
  ),
  railways = list(
    label="🚆  Ferrocarril", theme="transportation", type="segment",
    fields='id, names."primary" AS name, class, subtype, ST_AsText(geometry) AS geom_wkt',
    filter="AND subtype = 'rail'"
  ),
  waterways = list(
    label="🌊  Ríos / agua", theme="base", type="water",
    fields='id, names."primary" AS name, subtype, class, ST_AsText(geometry) AS geom_wkt'
  ),
  pois = list(
    label="📍  POIs", theme="places", type="place",
    fields='id, names."primary" AS name, basic_category AS category, confidence, ST_AsText(geometry) AS geom_wkt'
  ),
  addresses = list(
    label="📬  Direcciones", theme="addresses", type="address",
    fields='id, street, number, postcode, postal_city, country, ST_AsText(geometry) AS geom_wkt'
  ),
  lulc = list(
    label="🌿  LULC", theme="base", type="land_use",
    fields='id, names."primary" AS name, subtype, class, ST_AsText(geometry) AS geom_wkt'
  ),
  land = list(
    label="🏞  Parques / tierra", theme="base", type="land",
    fields='id, names."primary" AS name, subtype, class, ST_AsText(geometry) AS geom_wkt'
  ),
  divisions = list(
    label="🗺  Límites admin.", theme="divisions", type="division_area",
    fields='id, names."primary" AS name, subtype, admin_level, ST_AsText(geometry) AS geom_wkt'
  ),
  infrastructure = list(
    label="⚡  Infraestructura", theme="base", type="infrastructure",
    fields='id, names."primary" AS name, subtype, class, ST_AsText(geometry) AS geom_wkt'
  )
)
# Colores por capa (también usados en el control de capas del mapa)
LAYER_COLORS <- c(
  buildings="#f97316", roads="#facc15", railways="#a78bfa",
  waterways="#38bdf8", pois="#4ade80", addresses="#fb7185",
  lulc="#86efac", land="#6ee7b7", divisions="#c084fc", infrastructure="#fbbf24"
)
# ---- Helpers -----------------------------------------------------------------
get_latest_release <- function() {
  tryCatch(jsonlite::fromJSON("https://stac.overturemaps.org/catalog.json")$latest,
           error=function(e) "2026-04-15.0")
}
geocode_place <- function(place_name) {
  res <- tidygeocoder::geo(address=place_name, method="osm", full_results=FALSE, quiet=TRUE)
  if (nrow(res)==0 || is.na(res$long)) return(NULL)
  c(lon=res$long, lat=res$lat)
}
bbox_from_center <- function(lon, lat, buffer_km) {
  dlat <- buffer_km / 111.32
  dlon <- buffer_km / (111.32 * cos(lat * pi / 180))
  list(xmin=lon-dlon, ymin=lat-dlat, xmax=lon+dlon, ymax=lat+dlat)
}
circle_from_center <- function(lon, lat, buffer_km) {
  st_sfc(st_point(c(lon, lat)), crs=4326) |>
    st_transform(3857) |>
    st_buffer(buffer_km * 1000, nQuadSegs=32) |>
    st_transform(4326)
}
extract_layer <- function(con, layer_id, bbox, release) {
  cfg   <- LAYERS[[layer_id]]
  s3    <- glue("s3://overturemaps-us-west-2/release/{release}/theme={cfg$theme}/type={cfg$type}/*")
  extra <- if (!is.null(cfg$filter)) cfg$filter else ""
  sql <- glue("
    SELECT {cfg$fields}
    FROM read_parquet('{s3}', filename=true, hive_partitioning=1)
    WHERE bbox.xmin <= {bbox$xmax} AND bbox.xmax >= {bbox$xmin}
      AND bbox.ymin <= {bbox$ymax} AND bbox.ymax >= {bbox$ymin}
      {extra}
  ")
  tryCatch({
    df <- dbGetQuery(con, sql)
    if (nrow(df)==0) return(NULL)
    st_as_sf(df, wkt="geom_wkt", crs=4326)
  }, error=function(e) { message("[",layer_id,"] ",conditionMessage(e)); NULL })
}
clip_to_circle <- function(sf_obj, circle_sf) {
  tryCatch({
    cl <- st_intersection(st_make_valid(sf_obj), circle_sf)
    if (nrow(cl)==0) NULL else cl
  }, error=function(e) { message("clip: ",conditionMessage(e)); sf_obj })
}
write_gpkg_layer <- function(sf_obj, out_path, layer_name) {
  sf_obj <- sf_obj |> select(where(~!is.list(.x)))
  sf::st_write(sf_obj, out_path, layer=layer_name, driver="GPKG",
               delete_layer=TRUE, quiet=TRUE)
}
# Nombre de fichero por defecto basado en lugar y radio
default_filename <- function(place, km) {
  slug <- place |>
    tolower() |>
    gsub("[^a-z0-9áéíóúñ ]", "", x=_) |>
    trimws() |>
    gsub("\\s+", "_", x=_)
  paste0("~/overture_", slug, "_", km, "km.gpkg")
}
# ---- UI ----------------------------------------------------------------------
ui <- fluidPage(
  theme=overture_theme, useShinyjs(), useWaiter(),
  tags$head(tags$style(HTML("
    html, body { background:#0d1117; height:100%; overflow:hidden; }
    .container-fluid { height:100%; }
    .app-header { border-bottom:1px solid #21262d; padding:0.5rem 1.2rem 0.4rem; margin-bottom:0.5rem }
    .app-title { font-family:'Space Mono',monospace; font-size:1.2rem; color:#58a6ff; letter-spacing:-0.5px }
    .app-subtitle { color:#8b949e; font-size:0.7rem; margin-top:1px }
    /* columna izquierda: scroll interno si no cabe */
    .left-col {
      height: calc(100vh - 70px);
      overflow-y: auto;
      overflow-x: hidden;
      padding-right: 6px;
      scrollbar-width: thin;
      scrollbar-color: #30363d #0d1117;
    }
    .left-col::-webkit-scrollbar { width:4px }
    .left-col::-webkit-scrollbar-track { background:#0d1117 }
    .left-col::-webkit-scrollbar-thumb { background:#30363d; border-radius:2px }
    .card { background:#161b22!important; border:1px solid #21262d!important; margin-bottom:0.4rem!important }
    .card-header { background:#1c2128!important; border-bottom:1px solid #21262d!important;
      font-family:'Space Mono',monospace; font-size:0.72rem; color:#8b949e;
      text-transform:uppercase; letter-spacing:1px; padding:0.28rem 0.65rem!important }
    .card-body { padding:0.4rem 0.65rem!important }
    .layer-check .form-check { margin-bottom:0.12rem }
    .layer-check .form-check-label { font-size:0.78rem }
    .form-label { font-size:0.78rem; margin-bottom:0.12rem }
    .form-control { font-size:0.78rem; padding:0.18rem 0.4rem }
    .form-range { padding:0 }
    .btn-sm { font-size:0.75rem; padding:0.15rem 0.4rem }
    #geocode_btn { font-size:0.75rem; padding:0.18rem 0.4rem }
    #extract_btn { font-family:'Space Mono',monospace; background:#1f6feb; border:none;
                   font-size:0.78rem; padding:0.25rem 0.5rem; white-space:nowrap }
    #extract_btn:hover { background:#388bfd }
    #reset_btn { font-family:'Space Mono',monospace; background:#21262d;
                 border:1px solid #30363d; font-size:0.75rem; padding:0.22rem 0.5rem;
                 color:#8b949e; white-space:nowrap }
    #reset_btn:hover { background:#30363d; color:#e6edf3 }
    .btn-row { display:flex; gap:5px; margin-top:5px }
    .btn-row .btn { flex:1; min-width:0 }
    .log-box { background:#010409; border:1px solid #21262d; border-radius:6px;
      padding:0.4rem 0.65rem; font-size:0.67rem; color:#8b949e;
      font-family:'IBM Plex Mono',monospace;
      height:78px; overflow-y:auto; white-space:pre-wrap }
    .stat-badge { display:inline-block; background:#1c2128; border:1px solid #30363d;
      border-radius:4px; padding:1px 5px; font-size:0.63rem; color:#3fb950;
      margin:1px; font-family:'IBM Plex Mono',monospace }
    .release-badge { display:inline-block; background:#1f2937; border:1px solid #374151;
      border-radius:12px; padding:1px 7px; font-size:0.63rem; color:#d29922;
      font-family:'IBM Plex Mono',monospace }
    .leaflet-container { border-radius:6px }
    .irs--shiny .irs-bar { height:3px }
    .irs--shiny .irs-line { height:3px }
    .irs--shiny .irs-handle { top:22px }
    .irs-min, .irs-max, .irs-single { font-size:0.67rem }
  "))),
  
  div(class="app-header",
      div(class="app-title", "◈ OVERTURE MAPS EXTRACTOR"),
      div(class="app-subtitle",
          "by ",
          tags$a("Geovisualization.net", href="https://geovisualization.net",
                 target="_blank",
                 style="color:#58a6ff;text-decoration:none;"),
          " · DuckDB · GeoPackage · latest release: ",
          uiOutput("release_badge", inline=TRUE)
      )
  ),
  
  fluidRow(style="margin:0",
           # ---- Columna izquierda ---------------------------------------------------
           column(3, class="left-col", style="padding-left:6px",
                  
                  card(card_header("Área de extracción"), card_body(
                    textInput("place_name", "Lugar / ciudad", value="Madrid, España"),
                    sliderInput("buffer_km", "Radio (km)", min=0.5, max=25, value=3, step=0.5),
                    actionButton("geocode_btn", "📍 Geocodificar",
                                 class="btn btn-outline-primary btn-sm w-100 mb-1"),
                    div(style="font-size:0.7rem;color:#8b949e;line-height:1.4;margin-top:2px",
                        textOutput("bbox_info"))
                  )),
                  
                  card(card_header("Capas"), card_body(padding="0.35rem",
                                                       div(class="layer-check",
                                                           checkboxGroupInput("layers_sel", label=NULL,
                                                                              choices=setNames(names(LAYERS), sapply(LAYERS,`[[`,"label")),
                                                                              selected=c("buildings","roads","pois")))
                  )),
                  
                  card(card_header("Salida"), card_body(
                    textInput("filename_base", "Nombre del fichero", value="",
                              placeholder="auto (lugar_radio.gpkg)"),
                    textInput("out_dir", "Directorio", value="~"),
                    div(style="font-size:0.68rem;color:#8b949e;margin-bottom:4px",
                        textOutput("out_path_preview")),
                    div(class="btn-row",
                        actionButton("extract_btn", "⬇  EXTRAER", class="btn btn-primary"),
                        actionButton("reset_btn",   "↺  RESET",   class="btn")
                    )
                  )),
                  
                  card(card_header("Resumen"), card_body(uiOutput("stats_ui")))
           ),
           
           # ---- Columna derecha -----------------------------------------------------
           column(9, style="padding-left:6px",
                  card(card_header("Vista previa"),
                       card_body(padding=0,
                                 leafletOutput("preview_map", height="470px"))),
                  card(class="mt-2", card_header("Log"),
                       card_body(padding="0.35rem",
                                 div(class="log-box", id="log_box", textOutput("log_text"))))
           )
  )
)
# ---- Server ------------------------------------------------------------------
server <- function(input, output, session) {
  
  rv <- reactiveValues(
    lon=(-3.7038), lat=40.4168, bbox=NULL, circle=NULL,
    logs=character(0), results=list(), release="…"
  )
  
  output$release_badge <- renderUI(tags$span(class="release-badge", rv$release))
  observe({ rv$release <- isolate(get_latest_release()) })
  
  # Logger
  add_log <- function(msg) {
    rv$logs <- c(rv$logs, paste0("[",format(Sys.time(),"%H:%M:%S"),"] ",msg))
  }
  output$log_text <- renderText(paste(tail(rv$logs, 60), collapse="\n"))
  observe({
    rv$logs
    shinyjs::runjs("var l=document.getElementById('log_box');if(l)l.scrollTop=l.scrollHeight;")
  })
  
  # Ruta de salida calculada
  get_out_path <- reactive({
    base <- trimws(input$filename_base)
    dir  <- trimws(input$out_dir)
    if (dir=="") dir <- "~"
    if (base=="") {
      base <- default_filename(input$place_name, input$buffer_km)
      base <- basename(base)
    }
    if (!grepl("\\.gpkg$", base, ignore.case=TRUE)) base <- paste0(base, ".gpkg")
    path.expand(file.path(dir, base))
  })
  
  output$out_path_preview <- renderText({ get_out_path() })
  
  # Mapa base con control de capas
  output$preview_map <- renderLeaflet({
    leaflet(options=leafletOptions(preferCanvas=TRUE)) %>%
      # Basemap oscuro con etiquetas en inglés
      addProviderTiles("Stadia.AlidadeSmoothDark",
                       options=tileOptions(opacity=0.95)) %>%
      setView(lng=rv$lon, lat=rv$lat, zoom=13) %>%
      # Marker arrastrable en posición inicial
      addMarkers(
        lng=rv$lon, lat=rv$lat,
        layerId="center_marker",
        options=markerOptions(draggable=TRUE),
        label="Drag to reposition"
      ) %>%
      addLayersControl(
        overlayGroups = names(LAYERS),
        options = layersControlOptions(collapsed=TRUE)
      )
  })
  
  # Cuando el marker se suelta en nueva posición
  observeEvent(input$preview_map_marker_dragend, {
    info <- input$preview_map_marker_dragend
    req(info$id == "center_marker")
    rv$lon <- info$lng
    rv$lat <- info$lat
    rv$bbox   <- bbox_from_center(rv$lon, rv$lat, input$buffer_km)
    rv$circle <- circle_from_center(rv$lon, rv$lat, input$buffer_km)
    add_log(sprintf("📌 Marcador movido → %.5f, %.5f", rv$lon, rv$lat))
    
    # Reverse geocode para actualizar el campo de texto
    tryCatch({
      res <- tidygeocoder::reverse_geo(
        lat=rv$lat, long=rv$lon, method="osm", full_results=FALSE, quiet=TRUE
      )
      if (nrow(res)>0 && !is.na(res$address)) {
        updateTextInput(session, "place_name", value=res$address)
        add_log(paste("  ↳ lugar:", res$address))
      }
    }, error=function(e) NULL)
    
    updateTextInput(session, "filename_base", value="")
    update_map_circle()
  })
  
  # Geocodificar
  observeEvent(input$geocode_btn, {
    req(nchar(trimws(input$place_name))>0)
    add_log(paste("Geocodificando:", input$place_name))
    coords <- withProgress(message="Geocodificando…", value=0.5,
                           geocode_place(input$place_name))
    if (is.null(coords)) {
      add_log("⚠ Sin coordenadas. Prueba con nombre más específico")
      add_log("  ej: 'León, Castilla y León, España'  o arrastra el marcador")
      showNotification("Geocodificación fallida — prueba a arrastrar el marcador",
                       type="warning", duration=5)
      return()
    }
    rv$lon    <- coords["lon"]; rv$lat <- coords["lat"]
    rv$bbox   <- bbox_from_center(rv$lon, rv$lat, input$buffer_km)
    rv$circle <- circle_from_center(rv$lon, rv$lat, input$buffer_km)
    add_log(sprintf("✓ %.5f, %.5f  |  radio: %g km", rv$lon, rv$lat, input$buffer_km))
    updateTextInput(session, "filename_base", value="")
    update_map_circle()
  })
  
  observeEvent(input$buffer_km, {
    if (!is.null(rv$bbox)) {
      rv$bbox   <- bbox_from_center(rv$lon, rv$lat, input$buffer_km)
      rv$circle <- circle_from_center(rv$lon, rv$lat, input$buffer_km)
      update_map_circle()
    }
  })
  
  update_map_circle <- function() {
    b <- rv$bbox
    leafletProxy("preview_map") %>%
      clearShapes() %>%
      # Redibujar círculo
      addCircles(lng=rv$lon, lat=rv$lat, radius=input$buffer_km*1000,
                 fillColor="transparent", color="#58a6ff",
                 weight=1.5, dashArray="4") %>%
      # Mover marker arrastrable a nueva posición
      removeMarker("center_marker") %>%
      addMarkers(lng=rv$lon, lat=rv$lat,
                 layerId="center_marker",
                 options=markerOptions(draggable=TRUE),
                 label="Drag to reposition") %>%
      flyToBounds(b$xmin, b$ymin, b$xmax, b$ymax)
  }
  
  output$bbox_info <- renderText({
    if (is.null(rv$bbox)) return("No area set")
    sprintf("center: %.4f, %.4f  |  radius: %g km", rv$lon, rv$lat, input$buffer_km)
  })
  
  # Reset
  observeEvent(input$reset_btn, {
    rv$lon <- -3.7038; rv$lat <- 40.4168
    rv$bbox <- bbox_from_center(rv$lon, rv$lat, 3)
    rv$circle <- circle_from_center(rv$lon, rv$lat, 3)
    rv$results <- list()
    rv$logs <- character(0)
    updateTextInput(session, "place_name",    value="Madrid, España")
    updateTextInput(session, "filename_base", value="")
    updateTextInput(session, "out_dir",       value="~")
    updateSliderInput(session, "buffer_km",   value=3)
    updateCheckboxGroupInput(session, "layers_sel",
                             selected=c("buildings","roads","pois"))
    leafletProxy("preview_map") %>%
      clearShapes() %>%
      removeMarker("center_marker") %>%
      clearGroup(names(LAYERS)) %>%
      addMarkers(lng=-3.7038, lat=40.4168,
                 layerId="center_marker",
                 options=markerOptions(draggable=TRUE),
                 label="Drag to reposition") %>%
      setView(lng=-3.7038, lat=40.4168, zoom=13)
    add_log("↺ Reset — ready for new extraction")
  })
  
  # Extracción
  observeEvent(input$extract_btn, {
    req(!is.null(rv$bbox), !is.null(rv$circle), length(input$layers_sel)>0)
    
    layers   <- input$layers_sel
    bbox     <- rv$bbox
    circle   <- rv$circle
    release  <- rv$release
    out_path <- get_out_path()
    
    add_log("── Iniciando extracción ──────────────────────")
    add_log(paste("Release:", release))
    add_log(paste("Capas:", paste(layers, collapse=", ")))
    add_log(paste("Salida:", out_path))
    
    rv$results <- list()
    
    # ── DuckDB: conexión explícita con cierre garantizado ──────────────────
    # Usamos un bloque local + tryCatch para garantizar dbDisconnect
    # aunque falle cualquier query individual
    con <- NULL
    tryCatch({
      con <- dbConnect(duckdb::duckdb(), dbdir=":memory:")
      dbExecute(con, "INSTALL httpfs;  LOAD httpfs;")
      dbExecute(con, "INSTALL spatial; LOAD spatial;")
      dbExecute(con, "SET s3_region='us-west-2';")
      dbExecute(con, "SET memory_limit='3GB';")
      dbExecute(con, "SET threads=2;")          # evitar saturar RAM con paralelismo
      dbExecute(con, "SET http_timeout=120000;") # 120 s timeout S3
      add_log("✓ DuckDB listo")
    }, error=function(e) {
      add_log(paste("✗ DuckDB init error:", conditionMessage(e)))
      if (!is.null(con)) try(dbDisconnect(con, shutdown=TRUE), silent=TRUE)
      return()
    })
    
    # Limpiar capas anteriores del mapa
    leafletProxy("preview_map") %>% clearGroup(names(LAYERS))
    
    withProgress(message="Extracting from Overture…", value=0, {
      for (i in seq_along(layers)) {
        lid <- layers[i]
        add_log(paste0("→ ", LAYERS[[lid]]$label, "…"))
        incProgress(1/length(layers), detail=LAYERS[[lid]]$label)
        
        # ── 1. Query DuckDB (protegida individualmente) ──────────────────
        sf_raw <- tryCatch(
          extract_layer(con, lid, bbox, release),
          error=function(e) {
            add_log(paste0("  ✗ query error: ", conditionMessage(e)))
            NULL
          }
        )
        if (is.null(sf_raw)) {
          if (lid == "addresses") {
            add_log("  ↳ 0 features — addresses coverage is limited (Europe/US mainly)")
          } else {
            add_log("  ↳ 0 features (query)")
          }
          next
        }
        add_log(sprintf("  ↳ bbox: %s features", format(nrow(sf_raw), big.mark=".")))
        
        # ── 2. Clip al círculo (protegido; si falla usa bbox raw) ─────────
        sf_obj <- tryCatch({
          cl <- st_intersection(st_make_valid(sf_raw), circle)
          if (nrow(cl)==0) NULL else cl
        }, error=function(e) {
          add_log(paste0("  ⚠ clip falló, usando bbox: ", conditionMessage(e)))
          sf_raw   # fallback: sin clip
        })
        rm(sf_raw); gc()   # liberar memoria inmediatamente
        
        if (is.null(sf_obj)) { add_log("  ↳ 0 features (clip)"); next }
        add_log(sprintf("  ↳ clip: %s features", format(nrow(sf_obj), big.mark=".")))
        
        # ── 3. Guardar esta capa al GPKG inmediatamente ───────────────────
        # (no acumulamos en RAM — si el proceso muere lo ya guardado queda)
        gpkg_lyr <- LAYER_NAMES[[lid]]
        saved <- tryCatch({
          # Si la capa GPKG ya existe (ej: roads + railways → ROADS) hay que append
          existing_layers <- tryCatch(st_layers(out_path)$name, error=function(e) character(0))
          if (gpkg_lyr %in% existing_layers) {
            # Leer la existente, combinar, reescribir
            prev <- st_read(out_path, layer=gpkg_lyr, quiet=TRUE)
            combined <- bind_rows(prev, sf_obj |> select(where(~!is.list(.x))))
            rm(prev)
            sf::st_write(combined, out_path, layer=gpkg_lyr, driver="GPKG",
                         delete_layer=TRUE, quiet=TRUE)
            rm(combined)
          } else {
            write_gpkg_layer(sf_obj, out_path, gpkg_lyr)
          }
          TRUE
        }, error=function(e) {
          add_log(paste0("  ✗ GPKG write error: ", conditionMessage(e)))
          FALSE
        })
        
        if (saved) add_log(sprintf("  ✓ '%s' guardado", gpkg_lyr))
        
        # ── 4. Pintar en mapa y guardar referencia en rv$results ─────────
        # Pintamos ANTES de rm() para que el objeto esté en memoria
        tryCatch(
          paint_layer(lid, sf_obj),
          error=function(e) add_log(paste0("  ⚠ paint error: ", conditionMessage(e)))
        )
        rv$results[[lid]] <- sf_obj
        rm(sf_obj); gc()
      }
    })
    
    # Cerrar DuckDB limpiamente
    tryCatch(dbDisconnect(con, shutdown=TRUE), silent=TRUE)
    
    add_log("── Completado ────────────────────────────────")
    add_log(paste("GPKG:", out_path))
    showNotification(paste0("Saved: ", out_path), type="message", duration=8)
  })
  
  # Pintar una sola capa en el mapa (grupo = lid permite toggle)
  paint_layer <- function(lid, obj) {
    proxy <- leafletProxy("preview_map")
    col   <- LAYER_COLORS[[lid]]
    
    # Geometrías mixtas tras clip — forzar tipo dominante
    gtype <- as.character(sf::st_geometry_type(obj))
    is_point <- any(gtype %in% c("POINT","MULTIPOINT"))
    is_line  <- any(gtype %in% c("LINESTRING","MULTILINESTRING"))
    
    # Popup seguro: construir como vector de caracteres antes de pasar a leaflet
    popup_txt <- if ("name" %in% names(obj)) {
      paste0("<b>", LAYERS[[lid]]$label, "</b><br>", ifelse(is.na(obj$name), "", obj$name))
    } else {
      rep(paste0("<b>", LAYERS[[lid]]$label, "</b>"), nrow(obj))
    }
    
    if (is_point) {
      proxy %>% addCircleMarkers(data=obj, group=lid,
                                 radius=3, color=col, fillColor=col, fillOpacity=0.85, stroke=FALSE,
                                 popup=popup_txt)
    } else if (is_line) {
      proxy %>% addPolylines(data=obj, group=lid,
                             color=col, weight=1.4, opacity=0.85,
                             popup=popup_txt)
    } else {
      proxy %>% addPolygons(data=obj, group=lid,
                            color=col, weight=0.5, fillColor=col, fillOpacity=0.22,
                            popup=popup_txt)
    }
  }
  
  output$stats_ui <- renderUI({
    if (length(rv$results)==0)
      return(tags$span(style="color:#8b949e;font-size:0.75rem;", "Sin datos aún"))
    tags$div(lapply(names(rv$results), function(lid) {
      col <- LAYER_COLORS[[lid]]
      tags$div(tags$span(class="stat-badge",
                         style=paste0("border-left:3px solid ",col,";padding-left:5px;"),
                         sprintf("%s → %s  (%s)",
                                 LAYERS[[lid]]$label, LAYER_NAMES[[lid]],
                                 format(nrow(rv$results[[lid]]), big.mark="."))))
    }))
  })
  
  isolate({
    rv$bbox   <- bbox_from_center(rv$lon, rv$lat, 3)
    rv$circle <- circle_from_center(rv$lon, rv$lat, 3)
  })
}
shinyApp(ui, server)

The code is open, the data is open, and the workflow is as simple as it should be. You’re welcome! Built with R · Shiny · DuckDB · Overture Maps Foundation Desarrollado por Alberto Concejal / Geovisualization.net

Alberto C.

https://www.geopackage.org/
https://es.wikipedia.org/wiki/RStudio
https://overturemaps.org/
https://duckdb.org/

Dedicated to my fellow colleagues in Tunis. Five years, five brilliant minds, countless maps, and a shared drive to make every single day better than the last. You motivated me more than you know.

Latest projects developed by Alberto C. by 20260417

1. Precision Agriculture: NDVI/NDWI/ECI/SAVI. GEE /Javascript
https://geovisualization.net/2026/03/20/agricultura-de-precision-ii-app-para-integracion-con-catastro-rural-en-espana/
2. Remote Sensing. From 10m to 1m (x100 times in few lines of code). Python /Colab
https://geovisualization.net/2026/03/18/super-resolution-1-m-a-cadalso-de-los-vidrios-madrid-avec-sentinel-2-10m-magique/
3. Handling High resolution Terrain models in Mapterhorn. R
https://geovisualization.net/2026/03/17/setting-up-mapterhorn-terrain-in-rstudio/
4. RADAR Flood assessment. GEE /Javascript
https://geovisualization.net/2026/02/23/flood-assessment-using-radar-sentinel1-srtm-v3-and-ghsl-pop-estimation-in-spain/
5. From Lidar to DTM/DSM in using code. R
https://geovisualization.net/2026/02/16/from-lidar-usgs-to-dsm-in-a-few-lines-of-code-the-magic-of-r/
6. Big Data analysis of any subject (running). R
https://geovisualization.net/2026/01/27/con-r-de-running/
7. Flood assessment and Open Data. GEE /Javascript
https://geovisualization.net/2025/11/18/mapping-something-unthinkable-flood-risk-in-madrid-using-open-data/
8. Spatial patterns’ analysis on any subject (urban gambling). ArcGIS
https://geovisualization.net/2025/12/17/analyzing-spatial-correlation-between-purchase-power-index-and-gambling-stores-2/
9. LIDAR QC DTM/DSM/Open Data. Global Mapper
https://geovisualization.net/2025/09/05/precision-elevation-data-for-forest-giants-lidar-vs-eth-global-canopy-height-in-mata-do-bucaco-portugal/
10. WildFire assessment in NRT. GEE /Javascript
https://geovisualization.net/2025/08/20/se-nos-quema-la-peninsula-este-2025/
11. RADAR usage for Flood assessment. GEE /Javascript
https://geovisualization.net/2025/07/07/sentinel-1-sar-un-aliado-indispensable-para-el-analisis-y-seguimiento-de-inundaciones-derna-libia-2023/

1656 followers in LinkedIn is something (update 20260521)

Una red con nombre y apellidos. A veces, las métricas de las redes sociales nos hacen olvidar que detrás de cada clic hay una persona. Ver este camino de 1531 profesionales a vista de helicóptero no es un ejercicio de ego, sino de gratitud.

Seguidores de mi perfil de LinkedIN a mediados de Abril 2026

Si estás en esta imagen, es porque en algún momento hemos cruzado ideas, compartido un proyecto o simplemente encontrado un punto de valor común. Esta comunidad la formáis:

  • Colegas con los que comparto el rigor técnico.
  • Amigos que apoyan cada nueva iniciativa.
  • Seguidores que buscan una perspectiva diferente sobre cómo vivimos y nos movemos.

Entiendo esta red no como una audiencia, sino como un ecosistema de aporte mutuo. Si seguís aquí, es porque habéis encontrado algo de valor en lo que comparto, y eso es el mayor motor para seguir analizando, cuestionando y construyendo.

Gracias por ser parte de este camino.

+80,000 visitas en el último año!!!! (LinkedIn)

Alberto C.
Analista GIS, quizá colega, quizá amigo, quizá seguidor 🙂

AERIAL SURVEYOR SIMULATOR: Reimagining Aerial Photography

For three years, my office was thousands of feet in the air. As an aerial photographer, I spent my days capturing the world’s textures, layouts, and topographies from a cockpit—a masterclass in perspective that I am now transforming into a Aerial Surveyor Simulator.

https://aerialsurveyorsimulatorv2.netlify.app/
This tool bridges the gap between real-world mapping and immersive simulation, serving as a vital resource for aeronautical and cartographic training, drone mission familiarization, and the visualization of complex technical projects in a 2D environment.

Figure 1 – Aerial Surveying, back in 2000 🙂

To ensure maximum accessibility and performance, I built the entire system using a lightweight, “zero-dependency” stack:

  • Core: A single self-contained file using HTML5, Pure CSS3 (Grid/Flexbox) with native dark mode, and Vanilla JavaScript—no frameworks, no servers, no build processes.
Figure 2 – Screen captures of the whole SIMULATOR scenario

Alberto C.
Geospatial Analyst

#AerialPhotography #Cartography #VanillaJS #WebDev #Simulation #Geospatial #TechInnovation