From Overture Maps to GPKG in s: 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/

1531 followers in LinkedIn is something

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.

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.


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

Urban development in Madrid from the mid-19th century to the present day

Today I’m sharing a post from a while back about a project I worked on, looking at Madrid’s urban development since the 19th century.

All existing buildings in Madrid currently listed in the Land Registry database have their year of construction recorded. This map shows, by decade, where the bulk of that urban development took place. For example, in the 1920s it was in the Salamanca district, in the 1930s in Chamartín… shifting from development in the city centre to the outskirts.

Figure 1 – The data obtained from the land registry was used to provide content for this exhibition at the CANAL Foundation

CLIFFORD. VIEWS OF MADRID UNDER ISABEL II

The Canal Foundation presents ‘Clifford: Views of Madrid under Queen Isabella II’, an exhibition dedicated to the work of one of the pioneers of photography. It comprises a selection of almost a hundred images organised into four sections: ‘The Pleasures of Photography’, ‘Old and New Madrid’, ‘In the Service of the Monarchy’ and ‘The Construction of the Canal de Isabel II: a project worthy of the Romans’. The exhibition offers a glimpse of Madrid before the arrival of running water – a city very different from the one we know today, to which the most ambitious engineering project of its time managed, almost a century late, to open the doors to modernity.

Added Value of GIS and Geovisualization in Urban Evolutionary Analysis

  • Spatiotemporal Data Synthesis: GIS facilitates the integration of disparate cadastral datasets into a unified spatial framework, transforming longitudinal construction records into high-density visual intelligence.
  • Dynamic Morphological Tracking: The transition from static cartography to animated geovisualization (e.g., temporal GIF rendering) enables the identification of urban “pulses,” capturing the velocity and direction of sprawl, infill development, and densification cycles that remain obscured in traditional tabular formats.
  • Heatmap Analytics & Kernel Density Estimation: By utilizing spatial Interpolation and density-based clustering, GIS highlights hotspots of urban intensity, allowing for a comparative analysis between centralized vertical growth and peripheral horizontal expansion.
  • Infrastructure Correlational Mapping: Geovisual tools allow for the overlay of historical urban growth against transport network expansion, providing a technical basis for evaluating the causal relationships between infrastructure deployment and land-use evolution.
  • Interactive Chronological Filtering: Modern geovisualization interfaces support granular temporal querying, enabling researchers to isolate specific socioeconomic eras and visualize the physical footprint of policy-driven urban shifts in real-time.
Figure 2 – Exhibition leaflet
Figure 3 – GIS map using ArcGIS Desktop

Hope you guys like it!

Alberto C.
GIS Analyst

Sources:
https://www.sedecatastro.gob.es/
https://www.fundacioncanal.com/exposiciones/clifford-vistas-del-madrid-de-isabel-ii/

Bienvenido a Madrid Río, donde el carril único lo usa todo el mundo… pero solo la mitad lo usa bien

Madrid Río va a su bola. Y tengo los datos para demostrarlo. Madrid Río, una de las arterias verdes más importantes de la capital: casi 7 km de recorrido por sentido, desde el Puente de los Franceses hasta el Parque Lineal del Manzanares, sin semáforos, sin coches, y con una norma tan sencilla que casi da vergüenza tener que recordarla: ve por la derecha.

Diagrama 1 – Si no quieres atropellar o ser atropellado, ve por tu derecha 🙂

No hay señales. No hay carriles diferenciados para peatones, corredores, ciclistas, patinetes, patines o skates. Solo una convención tácita, no escrita, que funciona… cuando la gente la respeta. El problema es que mucha gente no la respeta.

Diagrama 2 – Madrid Rio esquema. Estudio Álvarez-Sala (2011)
Diagrama 3 – Si vas por tu derecha, el mundo funciona! (norma no escrita sugerida)

Y cuando eso ocurre, se genera un caos silencioso. Normalmente se resuelve en segundos con ese equilibrio dinámico tan mediterráneo que tenemos —un quiebro, una mirada, un gesto— pero no siempre. He visto bicis arrollar a personas que iban por su lado. Grupos andando en paralelo, bloqueando la vía entera. Patinetes a velocidad de vértigo a punto de llevarse por delante a un niño que, irónicamente, iba exactamente donde debía. Yo mismo llevo usando esta vía desde casi su inauguración, hace más de once años, y puedo decir que el problema no ha mejorado: ha crecido, a medida que el carril se ha llenado de nuevas tipologías de movilidad.

¿Qué quiero hacer? Quiero medir esto. Con rigor, con datos, y con herramientas actuales. Mi objetivo es retratar estadísticamente quién va en el sentido correcto y quién no, desagregando por tipo de movilidad, edad, género y comportamiento en grupo. ¿Son los ciclistas los que más incumplen, o los patinetes eléctricos? ¿La gente mayor es más respetuosa que los jóvenes? ¿Los grupos de tres o más personas son el mayor factor de caos? No lo sé aún. Pero lo voy a descubrir.

Para ello uso mis conocimientos de análisis GIS para generar mapas de calor y patrones espaciales. La app está construida con tecnologías web estándar, sin ningún framework ni dependencia externa, como no podía ser de otra forma, me he ayudado de una IA, no la digo por no hacer publicidad pero en este caso he usado la que empieza por CL y termina por AUDE.

  • *HTML5 — estructura y contenido. Todo en un único fichero autocontenido, sin necesidad de servidor ni build process.
  • *CSS3 puro — estilos, layout con Grid y Flexbox, animaciones, variables CSS para el sistema de colores y soporte automático de dark mode con prefers-color-scheme.
  • *avaScript vanilla — sin React, sin Vue, sin jQuery, sin nada. Todo el comportamiento, la lógica de captura, el simulador, los gráficos y la gestión de datos está escrito en JS nativo del navegador.

Para los elementos específicos:

  • *Blob + URL.createObjectURL — generación y descarga del CSV directamente desde el navegador, sin servidor
  • *Canvas API (nativa del navegador) — el donut chart de OK/NOK en estadísticas
  • *localStorage — persistencia de datos entre sesiones, sin base de datos ni backend
  • *Geolocation API (nativa del navegador) — captura de coordenadas GPS en tiempo real

¿Por qué me importa? No soy escandinavo. No espero perfección. Entiendo que en Madrid las cosas funcionan a su manera y que eso tiene su gracia, pero me gustaría que mis paseos —y los de mis hijos— por esta vía fueran un poco más seguros. Y si de paso consigo que alguien lea esto, lo comparta, y la próxima vez que salga a correr piense dos segundos antes de cruzarse al carril contrario… misión cumplida.

A veces hacer las cosas un poco mejor no requiere una ordenanza municipal. Solo requiere información. Este proyecto empezó hace 11 años con un cuaderno y un bolígrafo. Hoy empieza de nuevo, con mejores herramientas. Aquí tenéis el post de hace ya 11 años!

1 Lo primero es comprender qué es ir bien y qué no, para ello he construido un simulador que me permite explicar tanto a la gente que va por la vía como a un posible ayudante que vaya a tomar datos. Registrar observaciones en tiempo real no es tan fácil como parece. Cuando hay mucha gente pasando a la vez —una bici, dos peatones y un patinete en tres segundos— es fácil perder el hilo, contar dos veces a la misma persona o simplemente no dar abasto.

Diagrama 4 – Entrenando el modelo, cada cuál va por donde debe…

Por eso antes de ir al campo, me entreno. El simulador recrea una sección transversal de la vía: dos sentidos de circulación, con sus zonas correctas marcadas en verde, una franja central de duda en naranja, y figuras que cruzan la pantalla a distintas velocidades según su tipo de movilidad. El ritmo puede ser tranquilo, normal o caótico —porque Madrid Río un domingo a las 11 de la mañana es, efectivamente, caótico.

Vídeo 1 – El simulador de Madrid Río 🙂

La mecánica es simple: clic izquierdo si va bien, clic derecho si va mal. El objetivo no es acertar —no hay respuesta correcta predefinida, tú eres el juez— sino desarrollar el reflejo de observar, clasificar y registrar sin perder de vista al siguiente que ya está entrando en pantalla.

Cuando el ritmo sube y aparecen tres iconos a la vez yendo en direcciones distintas, te das cuenta de que el trabajo de campo requiere más concentración de la que parece.

2 Pasamos a la acción. Ahora toca capturar los datos. El mayor reto de este proyecto no es el análisis —es el momento exacto en que alguien cruza la sección de medición.

Tienes aproximadamente un segundo. A veces menos, por eso la herramienta de captura está diseñada en dos tiempos. En el primero, solo hay dos botones: OK o NOK. Nada más. Sin pensar en si es hombre o mujer, sin estimar la edad, sin identificar si va en patinete o en patín eléctrico. Solo: ¿va por donde debe o no? Un toque y listo. El registro queda guardado con timestamp y coordenadas GPS.

El segundo tiempo llega cuando hay una pausa —entre una oleada y la siguiente. Entonces abres cada registro de la cola, que aparece marcado en amarillo si le falta detalle, y lo enriqueces con calma: tipo de movilidad, edad aproximada, género, si iba solo o en grupo. Campos que en el momento del cruce son imposibles de procesar, pero que con dos segundos de margen se completan sin problema.


Vídeo 2 – Extracción de datos y análisis de resultados

El resultado es un sistema que no te hace elegir entre velocidad y riqueza de datos. Primero capturas, luego describes. Y si pasan cuatro personas a la vez —cosa que pasa— no pierdes ninguna por intentar ser exhaustivo con la primera.

Diagrama 5 – Detallando el perfil del usuario de Madrid Río (1)

Este es el modelo de datos, cada persona que atraviesa el punto en el que tomo los datos ha de responder a este formulario. Lo has hecho bien?, Iba andando, corriendo?, Oba solo, acompañado?, Qué edad tenía?, Cuál es su género?.

Diagrama 6 – Detallando el perfil del usuario de Madrid Río (y 2)

¿Y tú, por dónde vas?

Los datos que recoja en esta primera fase son míos, de un único observador, en un punto concreto de la vía. Eso tiene valor, pero tiene límites, lo que realmente haría este proyecto significativo es escala: más puntos de medición, más horas del día, más días de la semana, más condiciones meteorológicas. Madrid Río no se comporta igual un martes lluvioso de febrero que un domingo de abril con sol, si paseas, corres, vas en bici o llevas a tus hijos por Madrid Río con cierta regularidad y te apetece participar, la herramienta es libre, funciona desde el móvil sin instalar nada, y el único requisito es estar en el sitio y tener ganas de mirar.

Diagrama 7 – Exportando a CSV (ponderando rieso)
Diagrama 8 – Hotspot analysis – importando GeoJSON desde CSV

No hace falta ser analista de datos. No hace falta saber de GIS. Solo hace falta pararse cinco minutos en un punto del recorrido y registrar lo que pasa delante de ti.

Si te interesa colaborar, tienes dudas sobre la metodología, o simplemente quieres debatir si esto tiene solución o estamos condenados a esquivarnos eternamente, déjalo en los comentarios o escríbeme directamente.

Y la próxima vez que salgas a Madrid Río… ya sabes. Ve por la derecha, por ti, por el que viene de frente, y por el niño que va exactamente donde debe.

Diagrama 8 – Interfaz de la Aplicación

Espero que os haya parecido interesante,

Alberto C.
Analista GIS

https://eas.es/proyectos/madrid-rio/
https://habitat.aq.upm.es/dubai/14/bp-52.html
https://es.wikipedia.org/wiki/Madrid_R%C3%ADo
revistes.ub.edu/index.php/ScriptaNova/article/view/42419
https://redalyc.org/journal/3692/369242871009/html/

Diagrama 9 – Propuesta señalética

Population Estimation through Dynamic LULC-Based Settlement Validation

This methodology involves the deployment of a centralized processing interface within the Google Earth Engine (GEE) environment. The provided visualization captures the core interface of the custom GEE application, which serves as the hub for the multi-sensor LULC validation pipeline. Within this dashboard, users can define a specific Area of Interest (AOI)—highlighted here over Spain—and configure key parameters, including temporal ranges for the acquisition of sentinel-derived products. Crucially, the interface is designed to load and compare two primary datasets simultaneously: 10m Dynamic World (near real-time, probability-based LULC) and ESA WorldCover (10m resolution structured LULC). The contrasting classification schemes are represented by the legends on the left and right sides of the map view, which illustrate the varying definitions of ‘Built-up’ and urban areas between the two products. Establishing this visual and statistical comparison at the application level is the prerequisite for calculating the spatial disagreement threshold, or delta, that guides the subsequent merging and population estimation phases.

Figure 1 – Dynamic World and World Cover. LULC for pop estimation

In this stage of the workflow, the application performs a localized comparative analysis within the selected AOI (drawing a rectangle or freehand) cluster. The visualization highlights the spatial intersection between the Dynamic World and WorldCover datasets, where the red-shaded areas represent the high-confidence built-up zones confirmed by both sources.

A critical component of this dashboard is the automated calculation of the spatial disagreement (delta); as seen in the right-hand panel, the system evaluates the consistency between the LULC assets. If this delta remains below the 20% threshold, the urban pixels are merged to create a validated settlement mask. This mask then acts as the primary constraint for the population disaggregation phase. The interface displays the final output in real-time—in this instance, estimating a population of 9,643 individuals—by clipping the GHSL 2023 grid to the validated footprint. This ensures that the population count is geographically anchored only to verified human structures, providing a far more precise metric than the generalized figures typically found in administrative boundary datasets.

Figure 2 – Dynamic World extraction. If delta challenging them both is smaller than 20%, we go ahead

This third capture demonstrates the robustness of the methodology when applied to larger, more complex urban fabrics. In this scenario, the dashboard processes a more extensive metropolitan area, showcasing the comparative histogram between ESA WorldCover and Dynamic World (highlighted in the right-hand panel). These bar charts represent the pixel-wise distribution across various LULC classes, specifically comparing the extent of “Built-up” or “Urban” areas between the two assets. By analyzing these frequency distributions, we can quantitatively assess the spatial agreement before merging.

The visualization confirms that even at this broader scale, the logic remains consistent: the system validates the urban footprint by identifying where the two global products intersect geographically. As shown in the population indicator, the disaggregated estimate for this specific zone reaches 83,522 individuals. This illustrates the scalability of the GEE application, moving seamlessly from small clusters to large urban centers while maintaining the same rigorous validation threshold. The ability to visualize these distributions through real-time charting allows for immediate auditing of the spatial data quality, ensuring that the final population figures are derived from a high-confidence settlement mask rather than a single, potentially biased source.

Figure 3 – World Cover. LULC for pop estimation

This fourth visualization provides a critical look at the methodology’s sensitivity to agricultural-urban fringes and smaller, isolated settlements. In this specific Area of Interest (AOI), the application processes a more fragmented settlement pattern, where the population estimate is calculated at 3,009 individuals.

Figure 4 – GHSL23 pop estimation

The side-panel charting once again highlights the comparative analysis of LULC classes, emphasizing the “Urban” category (marked in the red circle on the histogram). This demonstrates the system’s ability to maintain a rigorous validation threshold even in peri-urban environments where spectral mixing between built-up areas and bare soil or crops is common. By requiring cross-product agreement before merging, the model avoids over-estimating population in rural-urban transition zones—a frequent error in static administrative datasets. The consistency of the population disaggregation, as seen here, proves that the workflow is robust across various settlement morphologies, ensuring that only structurally verified pixels are used to “anchor” the human presence data provided by GHSL 2023.

Figure 5 – Dynamic World and Worl Cover. LULC for pop estimation

In summary, the transition from rigid administrative boundaries to dynamic, LULC-validated population modeling represents a significant advancement in urban spatial analytics. By implementing a dual-source validation pipeline—leveraging the temporal frequency of Dynamic World and the structural precision of ESA WorldCover—we create a refined “Human Presence” mask that is far more representative of physical reality. This methodology effectively addresses the Modifiable Areal Unit Problem (MAUP) by anchoring population counts from GHSL 2023 to verified built-up pixels rather than generalized polygons. For researchers and urban planners, this granular approach is not merely a technical refinement; it is a vital tool for identifying socio-spatial vulnerabilities and ensuring that resource allocation is based on high-resolution, empirical evidence of where communities actually thrive. As urban landscapes continue to evolve at an unprecedented pace, these satellite-derived workflows provide the scalability and precision necessary to map the future of our cities accurately.

Hope you like it,

Alberto C.
Geospatial Analyst

Agricultura de precisión (II). APP para integración con Catastro rural en España

La convergencia entre los las aplicaciones geoespaciales actuales y la administración pública ofrece una oportunidad sin precedentes para la optimización agronómica. La capacidad de procesamiento de Google Earth Engine (GEE), vinculada a la cartografía vectorial del Catastro rural, permite transformar las series temporales de misiones como Sentinel-2 en herramientas de diagnóstico directo sobre la parcela. Este enfoque desplaza el análisis de una observación puramente visual a una monitorización cuantitativa basada en la respuesta espectral de los cultivos.

El núcleo de esta aplicación reside en la intersección geométrica de las parcelas catastrales con colecciones de imágenes multiespectrales. Mediante el uso de la API de JavaScript en GEE, se automatiza el cálculo de indicadores biofísicos críticos como el NDVI (Índice de Vegetación de Diferencia Normalizada), el NDWI (Índice de Agua de Diferencia Normalizada), el EVI (Índice de Vegetación Mejorado) y el SAVI (Índice de Vegetación Ajustado al Suelo). Estos índices no solo reflejan el vigor fotosintético, sino que permiten identificar anomalías de crecimiento, estrés hídrico o variaciones en la densidad foliar que son invisibles al ojo humano en las fases tempranas del ciclo fenológico.

La integración técnica busca simplificar la toma de decisiones complejas en el campo. Al normalizar los datos del Catastro, la herramienta facilita que cualquier usuario técnico pueda extraer un informe en formato PDF con la evolución histórica y actual de su explotación. Este reporte actúa como un cuaderno de campo digital, proporcionando la base científica necesaria para determinar con precisión las ventanas temporales óptimas para la siembra, la aplicación selectiva de fitosanitarios para el control de plagas o la programación de la cosecha según el estado de madurez y los niveles de humedad detectados en el terreno.

A satellite map view of the area around Cadalso de los Vidrios, featuring geographical details such as roads, vegetation indices, and cadastral reference information.
Figura 1 – Interfaz de la aplicación en GEE

El funcionamiento de la aplicación se articula mediante la consulta de la base de datos catastral (cargada en GEE como un ASSET), permitiendo al usuario introducir una REFERENCIA CATASTRAL específica o buscar por POLÍGONO/PARCELA. Una vez reconocida la geometría de la parcela, el sistema actúa como un motor de filtrado espacial y temporal que realiza una llamada a las colecciones de imágenes satelitales disponibles en el repositorio de Google Earth Engine. Esta arquitectura permite que, de forma automatizada, se extraigan los valores de los píxeles contenidos exclusivamente dentro de los límites de la propiedad, eliminando el ruido de las parcelas colindantes o de elementos infraestructurales que podrían sesgar los resultados del análisis agronómico.

La potencia de esta integración reside en la generación dinámica de gráficos de series temporales y leyendas explicativas que traducen los valores espectrales en estados fenológicos comprensibles. Al seleccionar una parcela, la interfaz despliega una comparativa de índices (NDVI para vigor, NDWI para humedad y SAVI/EVI para densidades foliares complejas) que facilitan la detección de patrones de estrés hídrico o anomalías de crecimiento. Esta capa de información técnica es la que finalmente se consolida en el informe PDF, ofreciendo un registro histórico y actual del comportamiento del cultivo. De este modo, la herramienta no solo visualiza datos espaciales, sino que proporciona un respaldo estadístico para programar tareas críticas como la fertilización, el tratamiento de plagas o la determinación de la ventana de cosecha más eficiente según la madurez hídrica del terreno.

A satellite image of a land parcel with a yellow outlined boundary. Attribute information is displayed on the right, including details like area, coordinates, and various identifiers.
Figura 2 – Visualización de capas SHAPE descargadas desde https://www.sedecatastro.gob.es/. Global Mapper

La interfaz se ha estructurado para priorizar la visualización simultánea de la cartografía y la estadística multiespectral. En el panel lateral, el usuario gestiona la entrada de datos mediante la referencia catastral, lo que activa de forma inmediata el renderizado de la parcela seleccionada y la carga de sus atributos espaciales. El visor central ofrece una representación de alta resolución del terreno, mientras que el área de análisis despliega gráficos interactivos que segmentan las series temporales por índices. Esta disposición permite correlacionar directamente las variaciones visuales sobre el mapa con las fluctuaciones en las curvas de vigor y humedad, facilitando una interpretación técnica rápida de la evolución fenológica del cultivo sin salir del entorno de trabajo.

Aerial map view displaying a defined agricultural plot outlined in red, with various vegetation types and nearby structures. Accompanying graphs show indices for vegetation health, including SAVI and EVI, over a specified time period.
Figura 3 – Resultado de la búsqueda por “Referencia Catastral”

El componente analítico se materializa en la generación de series temporales multianuales, donde la superposición de los índices permite una lectura cruzada de la salud del cultivo. Al graficar conjuntamente métricas como el NDVI y el NDWI, es posible identificar con rigor técnico el desfase entre el estrés hídrico y la pérdida de vigor fotosintético. Esta serie histórica no solo sirve para el diagnóstico actual, sino que permite establecer líneas base de comportamiento fenológico, fundamentales para predecir anomalías climáticas o evaluar la efectividad de tratamientos agronómicos previos a escala de detalle.

Line graph comparing NDVI, NDWI, SAVI, and EVI over time from May 2020 to January 2026, displaying fluctuations in values.
Figura 4 – Resultado de la búsqueda por “Referencia Catastral”. Detalle de extracción de todos los índices

La aplicación culmina con la generación de un informe PDF que traduce el Big Data en decisiones de campo inmediatas. Al automatizar la integración entre GEE y Catastro, se elimina la barrera técnica, permitiendo que cualquier usuario interprete los ciclos críticos de su explotación sin necesidad de procesamientos complejos. Es, en definitiva, convertir el rigor científico en una herramienta operativa, legible y directamente aplicable para maximizar la eficiencia en cada parcela.

Espero que lo encuentres interesante, si es así, manda un saludo al menos! La verdad es que hay un millón de cosas que se pueden hacer usando Javascript en GEE, construyendo APP operativas desde el minuto cero, orientadas a Smart Cities o en este caso más concretamente a Smart Towns 🙂

Alberto C.
Analista Geoespacial

catastro.gob.es @catastro_es #MinisterioTransformaciónDigital #GEE #AgriculturaDePrecisión #GIS #RemoteSensing #Teledetección #BigData #Sentinel2

Super-résolution 1 m a Cadalso de los Vidrios, Madrid avec Sentinel 2 (10m). Magique !

Passer d’une résolution de 10 mètres à 1 mètre change radicalement la perspective du suivi agricole : on ne regarde plus une parcelle dans sa globalité, on observe ce qui se passe à l’intérieur même des rangs de culture. Ce saut qualitatif est possible grâce à l’algorithme S2DR3, un modèle de Deep Learning qui ne se contente pas d’agrandir les pixels, mais reconstruit l’information manquante. En s’appuyant sur les corrélations entre les différentes bandes spectrales de Sentinel-2 et en s’entraînant sur des images de très haute résolution, l’IA parvient à synthétiser une image à 1 m/pixel d’une précision étonnante.

Figura 1 – Super-résolution Sentinel-2 à 1 m – Cadalso de los Vidrios, Madrid

Le processus est entièrement automatisé : le code interroge le catalogue Copernicus pour extraire la scène la plus récente sur Cadalso de los Vidrios, en écartant systématiquement les passages nuageux pour garantir une donnée pure. Pour un utilisateur sur geovisualization.net, cela signifie obtenir un fond de plan ultra-précis sans dépendre de vols de drones coûteux ou d’imagerie commerciale privée.

Une fois cette base de 1 mètre générée, l’exploitation en agriculture de précision devient chirurgicale. En calculant des indices comme le NDVI pour la vigueur végétative ou le NDWI pour le stress hydrique, on peut détecter des anomalies de croissance ou des besoins en irrigation sur des micro-zones spécifiques. L’utilisation d’indices plus fins comme le SAVI (qui ajuste l’influence du sol) ou l’EVI (plus sensible en zone de forte biomasse) permet de piloter chaque étape du cycle cultural.

Figura 2 – Extracción de índices de vegetación desde Sentinel -2
Figura 3 – S2DR3 Python Super-résolution Sentinel-2 à 1 m

On peut ainsi décider du moment exact de la plantation selon l’humidité résiduelle du sol, moduler l’apport d’engrais zone par zone, ou identifier les foyers de maladies avant qu’ils ne se propagent à toute l’exploitation. C’est le complément technique indispensable pour transformer la donnée satellite en une véritable aide à la décision sur le terrain.

# =========================================================
# S2DR3 - Super-résolution Sentinel-2 à 1 m
# Zone : Cadalso de los Vidrios, Madrid
# =========================================================
# --- 1. Monter Google Drive ---
from google.colab import drive
import os
from datetime import datetime
drive.mount('/content/drive')
# Nouveau dossier spécifique pour Cadalso
output_path = '/content/drive/MyDrive/Sentinel2_Cadalso1m'
!mkdir -p {output_path}
# Nettoyage et création du lien symbolique pour le moteur S2DR3
!rm -rf /content/output
!ln -s {output_path} /content/output
# --- 2. Installer le paquet S2DR3 ---
# Utilisation du wheel que vous avez validé
!pip -q install https://storage.googleapis.com/0x7ff601307fa5/s2dr3-20250905.1-cp312-cp312-linux_x86_64.whl
# --- 3. Importer le module principal ---
import s2dr3.inferutils
# --- 4. Définir la zone d'intérêt (Cadalso de los Vidrios) ---
# Coordonnées du centre du village.
# Le moteur S2DR3 traite généralement une zone autour de ce point (buffer inclus).
lonlat = (-4.44, 40.30)
# --- 5. Définir la date Sentinel-2 ---
# On utilise la date d'aujourd'hui pour forcer la recherche de la scène la plus récente
date_actuelle = datetime.now().strftime('%Y-%m-%d')
# --- 6. Lancer le traitement ---
# La fonction 'test' gère automatiquement le téléchargement de la scène
# la moins nuageuse proche de la date indiquée.
print(f"Lancement de la super-résolution pour Cadalso (Date cible : {date_actuelle})...")
s2dr3.inferutils.test(lonlat, date_actuelle)
print(f"✅ Traitement terminé. Les fichiers sont dans : {output_path}")

Le processus se déroule entièrement dans l’environnement Google Colab à l’aide d’un script Python.

Figura 4 – Sample 2 – S2DR3 Python Super-résolution Sentinel-2 à 1 m (Durham, UK)
Figura 5 – Sample 3 – S2DR3 Python Super-résolution Sentinel-2 à 1 m (Madrid, Spain)
Figura 6 – Sample 4 – S2DR3 Python Super-résolution Sentinel-2 à 1 m (Central Park, New York, USA)


J’espère que cela vous plaira,

Cordialement,
Alberto C.
Analyste géospatial

#AgroTech #S2DR3 #Télédétection #SuperResolution #AI #AgTech #Geovisualization #Sentinel2 #InteligenciaArtificia #Sentinel2 #Cadalso

https://www.sigterritoires.fr/index.php/es/uso-de-s2dr3-en-google-colab-para-el-estudio-de-los-corales-en-mauricio/
https://colab.research.google.com/
https://medium.com/@ya_71389/sentinel-2-deep-resolution-3-0-c71a601a2253
https://www.youtube.com/watch?v=zNX4vp1hGWI

Setting up Mapterhorn terrain in RStudio

¿Alguna vez has querido visualizar el relieve de un territorio en 3D directamente desde RStudio, sin depender de software GIS externo? Mapterhorn es un proyecto open source que distribuye modelos digitales de elevación (MDT) de alta resolución — hasta 2 metros en España — empaquetados en formato PMTiles, un estándar moderno que permite servir datos geoespaciales sin necesidad de un servidor propio.

En este post veremos cómo configurar Mapterhorn en R usando el paquete mapgl en Rstudio, que nos permite crear mapas interactivos con terreno 3D en pocas líneas de código. El resultado: visualizaciones como la que ves abajo, con sombreado de relieve (hillshade) generado directamente desde los datos de elevación del IGN.

No necesitas experiencia previa en cartografía — si sabes instalar paquetes en R, puedes seguir este tutorial.

Figura 1 – Empezamos a visualizar Mapterhorn

El siguiente paso es combinar el terreno de Mapterhorn con un estilo cartográfico que aporte contexto geográfico. En este ejemplo vemos la Sierra Norte de Madrid — con el embalse de Bustarviejo al fondo y municipios como Torrelaguna en primer plano — renderizada con una inclinación de cámara (pitch) de 60° y exageración de relieve moderada.

install.packages(c("mapgl", "terra", "elevatr", "sf", "tidyterra", "ggplot2"))
---------------------------------
library(mapgl)
library(terra)
library(elevatr)
library(sf)
maplibre(
style = carto_style("voyager"),
center = c(-3.703, 40.416),
zoom = 14,
pitch = 40
) |>
add_raster_dem_source(
id = "mapterhorn_pro",
url = "pmtiles://https://download.mapterhorn.com/planet.pmtiles",
tileSize = 512,
encoding = "terrarium"
) |>
set_terrain(source = "mapterhorn_pro", exaggeration = 0.1) |>
add_layer(
id = "sombras",
type = "hillshade",
source = "mapterhorn_pro",
paint = list(
"hillshade-exaggeration" = 0.5,
"hillshade-shadow-color" = "#333333",
"hillshade-illumination-direction" = 315
)
)

El resultado es un mapa donde el relieve deja de ser un dato abstracto y se convierte en algo inmediatamente legible: puedes ver de un vistazo la diferencia de altitud entre el fondo del valle del Jarama y las cumbres de la sierra, algo que un mapa plano convencional no transmite con la misma fuerza.

Figura 2 – Aportando contexto geográfico con un basemap (Google)
library(mapgl)
maplibre(
style = carto_style("voyager"),
center = c(-3.60, 40.75),
zoom = 11,
pitch = 60
) |>
add_raster_dem_source(
id = "mapterhorn",
url = "pmtiles://https://download.mapterhorn.com/planet.pmtiles",
tileSize = 512,
encoding = "terrarium"
) |>
set_terrain(source = "mapterhorn", exaggeration = 1.5)

Análisis de visibilidad: ¿se ven Cibeles y la Puerta de Alcalá entre sí?

Más allá de la visualización estética, los datos de elevación de Mapterhorn permiten hacer análisis geoespaciales reales. Un ejemplo clásico es el análisis de línea de visión (Line of Sight, LOS): dado un observador en un punto A, ¿puede ver el punto B sin que el terreno lo obstaculice?

En este caso trazamos el perfil de elevación entre la Fuente de Cibeles y la Puerta de Alcalá — apenas 500 metros de distancia en el centro de Madrid. La línea marrón representa el terreno real; la línea roja discontinua es la línea de visión teórica desde los ojos del observador (a 1,70 m de altura, ten en cuenta que la escala del eje Y está magnificada para mejor comprensión).

El resultado es claro: el terreno queda por encima de la línea de visión en todo el trayecto, lo que indica que el suelo sube progresivamente desde Cibeles hacia Alcalá. Curiosamente, aunque ambos monumentos son perfectamente visibles entre sí en la realidad (la calle es recta y despejada), el propio desnivel del terreno hace que la línea de visión rasante quede por debajo del suelo si no se tienen en cuenta los edificios ni la altura real de los monumentos.

Ahora comparemos nuestro Mapterhorn 2m con el Copernicus 30m, a pesar de que esta comparación es inconsistente más allá de una interpretación puramente visual.

Figura 4 – Visualizamos el DSM-ish 30m de Copernicus y el DTM-ish 2m del IGN (alerta de inconsistencia QC)
library(elevatr)
library(terra)
library(ggplot2)
# 1. Zona y descarga de ambos MDTs
puntos <- data.frame(x = c(-4.0, -3.8), y = c(40.7, 40.9))
dem_low <- rast(get_elev_raster(puntos, z = 8, prj = "EPSG:4326", clip = "bbox"))
dem_high <- rast(get_elev_raster(puntos, z = 14, prj = "EPSG:4326", clip = "bbox"))
# 2. Perfil de elevación (path profile) a lo largo de una transecta
pasos <- 200
lon_seq <- seq(-4.0, -3.8, length.out = pasos)
lat_seq <- seq(40.7, 40.9, length.out = pasos)
coords <- cbind(lon_seq, lat_seq)
z_low <- extract(dem_low, coords)[, 1]
z_high <- extract(dem_high, coords)[, 1]
dist_m <- seq(0, 25000, length.out = pasos) # ~25km de transecta
# 3. Path profile comparativo
perfil_df <- data.frame(
distancia = rep(dist_m, 2),
elevacion = c(z_low, z_high),
fuente = rep(c("Copernicus GLO-30 (~30m)", "Mapterhorn IGN (~2m)"), each = pasos)
)
ggplot(perfil_df, aes(x = distancia / 1000, y = elevacion, color = fuente)) +
geom_line(linewidth = 0.8, alpha = 0.85) +
scale_color_manual(values = c("Copernicus GLO-30 (~30m)" = "#e74c3c",
"Mapterhorn IGN (~2m)" = "#2c3e50")) +
labs(
title = "Perfil de elevación: Copernicus vs Mapterhorn",
subtitle = "Transecta SW-NE sobre la Sierra de Guadarrama",
x = "Distancia (km)",
y = "Altitud (m)",
color = NULL
) +
theme_minimal(base_size = 13) +
theme(legend.position = "top")
# 4. RMSE — remuestrea el MDT de alta res a la resolución del bajo
dem_high_resampled <- resample(dem_high, dem_low, method = "bilinear")
z_low_full <- values(dem_low)
z_high_full <- values(dem_high_resampled)
# Elimina NAs
idx <- complete.cases(z_low_full, z_high_full)
rmse <- sqrt(mean((z_high_full[idx] - z_low_full[idx])^2))
cat(sprintf("RMSE entre Copernicus GLO-30 y Mapterhorn IGN: %.2f metros\n", rmse))
# 5. Diferencia espacial (mapa de error)
diff_raster <- dem_high_resampled - dem_low
plot(diff_raster,
main = "Diferencia de elevación: Mapterhorn - Copernicus (m)",
col = hcl.colors(100, "RdBu", rev = TRUE))
Figura 5 – RMSE entre Copernicus GLO-30 y Mapterhorn IGN: 9.34 metros (alerta de inconsistencia QC)

Pero como avanzaba, Copernicus GLO-30 y Mapterhorn IGN no son dos mediciones del mismo fenómeno a distinta resolución — son dos productos con orígenes, metodologías y fechas de captura completamente diferentes. Copernicus GLO-30 se deriva del radar de apertura sintética del satélite TanDEM-X, que mide la superficie de la Tierra desde el espacio con una resolución nominal de 30 metros. Al ser una medición radar, captura la parte superior de la vegetación y los edificios, no el suelo desnudo — lo que en cartografía se llama un Modelo Digital de Superficie (MDS), no un Modelo Digital del Terreno (MDT) estricto.

Mapterhorn distribuye los datos del IGN España, obtenidos mediante vuelos LiDAR que permiten filtrar la vegetación y los edificios para obtener el suelo real. El resultado es un MDT de 2 metros de resolución que representa la topografía con una precisión altimétrica muy superior.

Por tanto, el RMSE que calculamos no mide únicamente el error de Copernicus — mide la diferencia acumulada entre dos modelos distintos, que incluye diferencias reales de resolución, diferencias metodológicas y el efecto de la vegetación sobre el radar. En zonas boscosas como la sierra de Guadarrama, estas diferencias pueden ser especialmente pronunciadas, llegando a varios metros en las áreas más densamente arboladas.

Figura 6 – Mapterhorn en Rstudio

Lo que hemos visto en este post es solo el punto de entrada. Una vez que tienes los datos de elevación conectados, las posibilidades se multiplican rápidamente.

En el ámbito del análisis de terreno puedes calcular pendientes y orientaciones para estudios de riesgo de erosión, planificación de infraestructuras o modelado de radiación solar. A partir del mismo MDT puedes delimitar cuencas hidrográficas automáticamente, trazar la red de drenaje teórica o calcular el área de captación de cualquier punto del territorio. Para proyectos de energía renovable, cruzar las orientaciones de ladera con datos de irradiación permite identificar zonas óptimas para instalaciones solares con pocas líneas de código.

En visualización avanzada puedes combinar el terreno 3D con cualquier capa vectorial propia — rutas de senderismo, límites administrativos, localizaciones de campo — y exportar el resultado como HTML interactivo para publicar directamente en la web. El paquete mapgl permite animar la cámara con fly_to() para crear recorridos virtuales sobre el terreno, algo especialmente útil para presentaciones o divulgación. También puedes generar mapas de sombreado artístico cambiando el ángulo de iluminación para conseguir efectos visuales dignos de cartografía profesional.

Para aplicaciones más especializadas, los datos de Mapterhorn son la base ideal para análisis de intervisibilidad a escala regional — por ejemplo, calcular desde cuántos puntos del territorio es visible una antena, un parque eólico o una torre de vigilancia de incendios. En urbanismo y ordenación del territorio permiten modelar el impacto visual de nuevas construcciones o calcular perfiles de ruido considerando la topografía real.

Todo esto sin salir de RStudio, sin licencias de software GIS y con datos de la máxima resolución disponible para España.

Mapterhorn ha eliminado uno de los principales obstáculos que tenían los usuarios de R para trabajar con datos de elevación de alta resolución: la complejidad de acceder, descargar y procesar MDTs de fuentes institucionales. Con una sola URL y tres líneas de código tienes acceso a los mejores datos topográficos disponibles para España y buena parte de Europa, listos para visualizar o analizar.

Espero que te haya interesado,

Alberto C.
GIS Analyst

#RStats #Mapterhorn #GIS #Cartografía #MapGL #TerrainAnalysis #OpenData #Geovisualización #RStudio #LiDAR #IGN #SpatialData #Mapping #DataViz #Topografía #RemoteSensing #OpenSource #GeospatialR #MDT #Hillshade

https://mapterhorn.com/
https://protomaps.com/blog/mapterhorn-terrain/
https://github.com/mapterhorn
https://posit.co/download/rstudio-desktop/
https://centrodedescargas.cnig.es/