Featured

POI Intelligence for Urban Asset Analysis in RStudio: assetIQ

When analysing urban ASSETS, there’s true VALUE in disaggregating Points of Interest (POIs) into hyper-local layers, separating retail density, financial hubs, and hospitality ecosystems into distinct, relative indicators. This is as sharp as it sounds, my app, assetIQ was built to add value to your POIS.

https://albertogis.shinyapps.io/TOIconPOI

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_20260528 Descarga

A Geospatial solution for everybody!

Diagram 8 – Adding value to ASSET over Sydney, Australia

Diagram 9 – Adding value to ASSET over London, UK

Diagram 10 – Adding value to ASSET over Berlin, Germany

Diagram 11 – Adding value to ASSET over NY, USA

Diagram 12 – Adding value to ASSET. All 4 samples

Diagram 13 – Adding value to ASSET over Cape Town, SA

Diagram 14 – Exporting poiq attribute within building footprint for 2.5D visualization of POIs theme

ct_poiq_Service_Overture_maps_2.5D_Descarga

Diagram 15 – Correlation matrix of all poiq among each other

Diagram 17 – Highest Correlation rates. Bilbao, Spain

Diagram 18 – Highest Correlation rates. Bilbao, Spain

Alberto C.
GIS analyst

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

Start · Engage · Match. Un simulador de compatibilidad social en un bar de verdad (sin decirse ni una palabra)

1. Resumen ejecutivo

Start Engage Match es un simulador interactivo, autocontenido en un único archivo HTML/JavaScript, que modela cómo un grupo de personas con gustos, orientaciones y niveles de paciencia distintos se cruza en el espacio físico de un bar y, en determinadas condiciones, genera afinidad y —ocasionalmente— una conexión completa («match»).

A lo largo de este documento, la aplicación y el local que representa comparten un único nombre: Start · Engage · Match, que además describe con precisión sus tres fases —se pulsa Empezar, se produce el Engage al cruzarse dos personas, y algunas veces surge el Match—.

https://start-engage-match.netlify.app

El sistema combina cuatro capas de simulación que se ejecutan en tiempo real sobre un lienzo (canvas) de 20×20 metros: un modelo espacial del local (barra, pista de baile, zonas de sofás), un modelo de comportamiento individual (movimiento hacia puntos de interés, paciencia, abandono), un modelo de interacción social (el «engage» o encuentro, con una capa de compatibilidad y una capa de atracción física aleatoria) y un modelo temporal y económico que reproduce el ciclo de una noche real de bar, de 17:00 a 05:00.

El resultado se expresa como una animación en vivo, un panel de estadísticas con métricas de negocio (ingresos, coste por match) y un registro narrativo de los encuentros completos, generado automáticamente.

2. Portada de la aplicación

Al abrir el archivo, antes de acceder al simulador, se muestra una pantalla de portada con el logotipo de la aplicación y un botón «Entrar» que da paso a la sala.

El logotipo es un lockup horizontal —icono y nombre en la misma línea, como cualquier marca convencional—: a la izquierda, un pequeño icono con dos siluetas de género indefinido (sin rasgos que las identifiquen como hombre o mujer), cada una de un color distinto (ámbar y violeta, los mismos tonos que identifican la barra y la pista de baile dentro del propio simulador), con los morros encontrándose y un pequeño corazón en el punto de contacto; a la derecha, el nombre «Start · Engage · Match». El mismo icono, en una versión más pequeña, se repite junto al título dentro de la aplicación para mantener la identidad visual coherente en toda la interfaz. Es una pieza gráfica simple, generada íntegramente en SVG vectorial dentro del propio archivo, sin imágenes externas.

3. Objetivo y motivación

El proyecto nace como ejercicio de modelización de un fenómeno social —el encuentro casual entre personas— usando herramientas propias del análisis espacial y la simulación basada en agentes (agent-based modelling), aplicadas aquí a un dominio lúdico en lugar de a un dominio geográfico convencional. ¿Todo es GIS?. La respuesta es no, ¡pues eso!

Cada persona se comporta como un agente autónomo con reglas de movimiento, atributos propios y una función de decisión (el «engage») que determina si, al cruzarse físicamente con otro agente, se produce algún tipo de conexión. El objetivo explícito no es predecir comportamiento real, sino disponer de una maqueta interactiva y visual sobre la que iterar reglas de compatibilidad, probabilidad y economía, viendo el efecto agregado de cada cambio de forma inmediata.

4. Modelo del espacio físico

El local se representa como un cuadrado de 20 × 20 metros (560 × 560 píxeles en pantalla). Dentro de ese cuadrado se han definido cuatro zonas de interés, cada una con su propia probabilidad de ser elegida como destino por los agentes, y dos puertas físicamente separadas:

Barra: franja a lo largo de la pared superior (18 × 2,4 m). Probabilidad de ser destino: 30 %.
Pista de baile: zona central (10 × 7 m), la de mayor densidad esperada. Probabilidad: 40 %.
Sofás y Mesas altas: dos rincones en las esquinas inferiores (4,5 × 3,5 m cada uno). Probabilidad: 12 % cada uno.
Deambular libre: un punto aleatorio del local, sin zona asociada. Probabilidad: 6 %.

La Entrada está fijada en el centro de la pared inferior; toda persona nueva accede siempre por ese punto. La Salida está fijada en la pared izquierda, en una posición distinta de la entrada. Esta separación física busca un flujo más orgánico: quien entra atraviesa la sala hacia su primer destino, y quien se marcha recorre el local en sentido contrario hacia una puerta distinta, en lugar de desaparecer por el borde más cercano en el momento de irse.

5. Modelo de la población

Cada persona (agente) se genera con los siguientes atributos, fijados en el momento de entrar y constantes durante toda su estancia salvo que se indique lo contrario:

Nombre: asignado aleatoriamente de un listado de nombres españoles, distinto para cada género.
Género: masculino o femenino, con probabilidad 50/50.
Orientación sexual: heterosexual, homosexual o bisexual (ver apartado 7).
Preferencias personales: un vector de 5 gustos booleanos (Viajar, Cocina, Estudiar, Deporte, Música), cada uno asignado de forma independiente con probabilidad 50 %. Estas preferencias son la base de la compatibilidad (apartado 9).
Historial de intentos: el conjunto de personas con las que ya se ha intentado un encuentro (con o sin éxito), para no repetir nunca el mismo emparejamiento.
Paciencia acumulada: contador de tiempo en el local sin lograr un match completo (apartado 10).

6. Modelo de movimiento

En lugar de un paseo aleatorio uniforme por todo el local —que en la práctica apenas generaría cruces entre pocas personas repartidas en 400 m²—, cada agente elige un punto de destino concreto dentro de una de las zonas de interés (ver apartado 4) y camina hacia él en línea recta con velocidad propia (0,5–1,0 m/s equivalentes).

Al llegar a su destino, el agente permanece inmóvil un intervalo aleatorio (entre 0,7 y 2,2 segundos simulados aproximadamente) —simulando que pide en la barra o baila— y a continuación elige un nuevo destino, repitiendo el ciclo. Esta atracción hacia puntos calientes concretos (sobre todo la pista de baile) es lo que produce cruces frecuentes y realistas incluso cuando la población activa es reducida.

7. Modelo de atracción y orientación sexual

Antes de intentar cualquier encuentro, el sistema comprueba si existe atracción mutua real entre las dos personas que se han cruzado físicamente, en función de su género y orientación:

Heterosexual: atracción únicamente hacia el género opuesto.
Homosexual: atracción únicamente hacia su mismo género.
Bisexual: atracción hacia ambos géneros.

La atracción debe ser mutua: si una persona no está interesada en el género de la otra, no se produce ningún intento de encuentro, aunque se hayan cruzado físicamente; simplemente siguen su camino, y ese cruce queda registrado como «sin interés» para no volver a evaluarse entre esas dos mismas personas. Esto significa, por ejemplo, que dos hombres heterosexuales que se crucen nunca iniciarán un encuentro entre sí, mientras que dos mujeres bisexuales sí podrán hacerlo.

El porcentaje de población no heterosexual y el reparto entre bisexualidad y homosexualidad dentro de ese porcentaje son parámetros configurables (apartado 16).

8. El encuentro: «engage» y «challenge» visual

Cuando dos personas con atracción mutua se cruzan físicamente (sus círculos se solapan en el lienzo) y no se han probado antes entre sí, ambas quedan congeladas en el sitio durante aproximadamente 1,8 segundos: es el «challenge». Durante ese tiempo se dibuja una línea de conexión entre ambas y, de forma escalonada, van apareciendo cinco casillas de color —una por cada una de las cinco preferencias—, iluminadas si ambas personas coinciden en esa preferencia (ya sea porque a las dos les gusta o porque a ninguna le gusta) y apagadas si no coinciden.

Transcurrido ese tiempo se resuelve el encuentro (apartado 9) y ambas personas quedan marcadas como ya intentadas entre sí, de modo que nunca vuelven a evaluarse la una a la otra.

Periodo de gracia de entrada: durante los primeros 5 segundos desde que una persona ha entrado por la puerta, no puede iniciar ningún encuentro. Sin esta pausa, dado que todo el mundo entra por el mismo punto, la mayoría de los encuentros se producirían de forma artificial justo en el umbral de la puerta, antes de que la gente tuviera ocasión de dispersarse hacia la barra, la pista o los rincones.

9. Compatibilidad y «chispa»: el modelo de dos niveles

La resolución del encuentro se calcula en dos pasos independientes, lo que da lugar a tres desenlaces posibles:

Paso 1 — Compatibilidad de intereses: se cuenta en cuántas de las 5 preferencias coinciden ambas personas. Si ese número iguala o supera un umbral configurable (por defecto 4 de 5, ajustable por franja horaria; ver apartado 11), se consideran «compatibles».

Paso 2 — La chispa: solo si son compatibles, se lanza una probabilidad —la «chispa»— configurable (50 % por defecto). Solo si esta también se cumple se considera un match completo.

Los tres desenlaces posibles son, por tanto:

Match completo (compatibles + chispa): aparece un destello dorado con un corazón, ambas personas brillan, se muestra su nombre y se genera una entrada en el registro narrativo (apartado 13). Salen juntas del local.
Compatible sin chispa (compatibles, pero sin suerte): destello azul, no hay match; ambas personas siguen su noche con normalidad, contando ese intento como uno de los tres disponibles.
Sin compatibilidad: destello rojo; mismo efecto práctico que el caso anterior.

Esta separación entre «compatibilidad objetiva» y «chispa aleatoria» evita que el modelo sea puramente determinista: dos personas con gustos idénticos no están garantizadas de conectar, igual que en la vida real la afinidad de intereses no siempre se traduce en atracción.

10. Condiciones de abandono del local

Cada persona puede abandonar el bar por cuatro vías distintas, todas ellas evaluadas continuamente:

Tres intentos agotados: si una persona acumula 3 encuentros sin lograr un match completo, ya no puede conseguirlo esa noche. Se le concede un breve intervalo de «terminando la copa» (entre 0,8 y 2,2 segundos) antes de dirigirse a la salida.
Paciencia agotada: si pasa un tiempo prudencial (configurable, 25 segundos reales por defecto) sin lograr ningún match completo, se cansa de esperar y abandona, con el mismo intervalo de cortesía anterior.
Sin candidatos disponibles: si a una persona no le queda nadie activo en el local con quien no haya probado ya (o con quien no exista atracción mutua), abandona de inmediato, sin ningún intervalo de cortesía —no tiene sentido que siga esperando si no hay nadie más con quien intentarlo—. Este es el caso típico de las dos últimas personas del local que ya se probaron entre sí sin éxito.
Cierre del local: a las 05:00 (apartado 11) se detiene la entrada de gente nueva; quienes ya estén dentro siguen su curso normal hasta resolverse, momento en el que la simulación concluye.

11. Modelo temporal: el ciclo de una noche

La simulación reproduce el horario de un bar real, de 17:00 a 05:00 (12 horas), comprimido en un intervalo de tiempo real configurable (5 minutos por defecto). El reloj interno se muestra en pantalla y determina tres franjas horarias con reglas propias:

17:00 – 00:00 («ambiente tranquilo»): se aplican el umbral de compatibilidad y la probabilidad de chispa configurados por el usuario, sin modificación.
00:00 – 04:00 («late», «la cosa se anima»): el umbral de compatibilidad baja automáticamente a un máximo de 3 de 5 (es decir, hace falta coincidir en menos cosas para ser compatible), y la probabilidad de chispa aumenta de forma gradual y lineal desde el valor configurado hasta un 90 %, a medida que se acerca la última hora.
04:00 – 05:00 («última hora»): el umbral se mantiene en 3 de 5 y la chispa se produce siempre (100 %) si hay compatibilidad, reproduciendo la conocida sensación de «cuanto más tarde, menos exigente se es».

Visualmente, el local se tiñe de un violeta muy sutil durante la franja intermedia y de un tono cálido/rojizo durante la última hora, a modo de ambientación de «última copa».

12. Modelo económico

La entrada al local tiene un coste fijo por persona: 5 € hasta medianoche, y el doble — 10 € — desde las 00:00 en adelante. Una vez dentro, la persona no vuelve a pagar y puede permanecer el tiempo que desee, sujeto únicamente a las condiciones de abandono del apartado 10.

El panel de estadísticas (la «cuenta del bar») acumula en tiempo real los ingresos totales realmente cobrados —no una simple multiplicación, sino la suma de lo que pagó cada persona según la hora exacta en la que entró— y calcula un coste medio por match completo (ingresos totales entre número de matches completos), como indicador simplificado de la rentabilidad del local en relación con su función social.

13. Generación narrativa del registro de la noche

Cada vez que se produce un match completo, el sistema compone automáticamente una breve crónica combinando los nombres reales de ambas personas, una frase de apertura elegida al azar de un banco de diez variantes (miradas, silencios, una sonrisa de más), una frase de cierre que menciona los gustos concretos que ambos comparten —o, si no comparten ningún gusto en positivo, que coinciden en lo que no les va—, una mención a la orientación de cada persona cuando no es heterosexual (por ejemplo, «(ambos son gais)» o «Marta es bisexual»), y finalmente un epílogo elegido al azar de un banco de seis variantes con distintos grados de insinuación, siempre sugerente y nunca explícito.

Estas crónicas se muestran en una sección con desplazamiento propio dentro del panel de estadísticas, conservando las últimas 25 entradas.

14. Ambientación sonora

El simulador incorpora una pista de música electrónica generada en tiempo real mediante la Web Audio API del navegador —no se reproduce ningún archivo de audio externo, evitando así cualquier cuestión de derechos—. El motor sintetiza un bombo (kick) con caída de frecuencia, un bajo en patrón repetitivo de cuatro notas y un hi-hat de ruido filtrado, todo ello a 126 pulsaciones por minuto en bucle continuo, con un interruptor para activarla o silenciarla.

15. Panel de métricas

El panel «Cuenta del bar» muestra en todo momento: entradas totales, personas en sala, número de encuentros intentados, cruces sin interés por incompatibilidad de orientación, encuentros compatibles, matches completos, personas que abandonaron sin match, la tasa de match completo por encuentro, el precio de entrada vigente, los ingresos totales acumulados y el coste medio por match completo.

16. Parámetros configurables

Todos los parámetros del modelo son ajustables en tiempo real mediante controles deslizantes, lo que permite explorar el efecto de cada regla sin modificar el código:

Parámetro	Rango / valor por defecto	Efecto
Personas	6 – 60 (24 por defecto)	Tamaño de la población simulada
Umbral match	2 – 5 sobre 5 (4 por defecto)	Coincidencias necesarias para ser «compatibles» (antes de medianoche)
Chispa	10 % – 90 % (50 % por defecto)	Probabilidad de match completo una vez hay compatibilidad
Paciencia	10 – 60 s (25 s por defecto)	Tiempo sin match completo antes de abandonar por cansancio
% no heterosexual	0 – 30 % (12 % por defecto)	Proporción de la población gay, lesbiana o bisexual
De ese %, % bi	0 – 100 % (50 % por defecto)	Reparto entre bisexualidad y homosexualidad exclusiva
Noche dura	3 – 20 min (5 min por defecto)	Minutos reales que representan las 12 horas de apertura
Reponer gente	activado / desactivado	Si se genera gente nueva para sustituir a quien se va
Música	activada / desactivada	Pista electrónica de ambientación

17. Notas de implementación técnica

El simulador es un único archivo HTML autocontenido (HTML, CSS y JavaScript en el mismo documento), sin dependencias externas ni conexión a internet, por lo que funciona igual abierto como archivo local que servido desde un navegador. El motor de animación usa Canvas 2D con un bucle basado en requestAnimationFrame; la física de movimiento, las colisiones y la resolución de encuentros se recalculan en cada fotograma sobre una lista de agentes en memoria, sin frameworks ni librerías externas.

Esta sencillez deliberada facilita que cualquier regla del modelo —umbrales, probabilidades, zonas, horarios— pueda ajustarse directamente en el código o mediante los controles de la interfaz, sin necesidad de recompilar ni de dependencias de terceros.

18. Limitaciones y posibles líneas futuras

El modelo de compatibilidad es simétrico y binario por preferencia (gusta / no gusta); no pondera la intensidad de cada gusto ni introduce afinidades parciales.
La orientación bisexual se modela como atracción indiferenciada a ambos géneros, sin matices adicionales.
El coste de entrada y los ingresos son indicadores simplificados; no se ha modelado gasto en consumiciones ni otros ingresos del local.
Posibles ampliaciones: perfiles de personalidad más ricos (más de 5 preferencias, pesos distintos), grupos de amigos que entran y se mueven juntos, eventos especiales (barra libre, DJ invitado) que alteren temporalmente las reglas, o exportación de los resultados de una noche simulada a un informe descargable.

Ahora, si os gusta, a a contarlo,

Alberto Concejal
GIS Analyst

KALMAN RADAR TRACKER: SEGUIMIENTO DE BLANCOS AÉREOS

Cuando alguien me pregunta sobre radar, pienso sobre todo en radares montados en satélites (sesgo geospacial) pero en realidad hay mucho más, hoy voy a hablaros de de radares aeroportados, de filtros de Kalman y seguimiento de blancos aéreos en movimiento… ¡Qué interesante!

Imagen 1- KALMAN RADAR TRACKER – El vuelo del blanco

Lo primero que pienso no es en el radar en sí, sino en el problema que resuelve, porque ese problema lo llevo resolviendo de otra forma desde hace años sin llamarlo por su nombre técnico. Un radar mide la posición de un avión con ruido. Un GPS mide la posición de un coche con ruido. Un sensor SAR mide el desplazamiento del terreno con ruido. En los tres casos hay una señal real escondida detrás de mediciones que saltan, que tiemblan, que nunca coinciden exactamente con la trayectoria verdadera. Y en los tres casos la respuesta es la misma matemática: combinar lo que predice el modelo físico con lo que dice el sensor, ponderando cada fuente según cuánto te fías de ella.

Eso es un filtro de Kalman, despojado de jerga. Llevo construyendo herramientas geoespaciales que rondan esta misma idea sin que nadie me lo pidiera explícitamente. Cuando trabajo con series temporales de NDVI en Earth Engine y aplico un suavizado para separar la tendencia real del ruido atmosférico de cada imagen Landsat, estoy haciendo una versión simplificada de lo mismo. Cuando proceso cambios en backscatter de Sentinel-1 sobre una mina y tengo que decidir qué variación es señal geológica y qué es ruido del sensor, vuelvo a estar en el mismo terreno conceptual. La estadística no cambia, cambia el dominio de aplicación.

Imagen 2- KALMAN RADAR TRACKER – Seguimiento de blancos aéreos

Lo interesante de pasar unos días metido en simulación de trayectorias radar fue confirmar hasta qué punto esto es transferible. Construí un simulador que genera el vuelo de un avión con maniobras aleatorias, le añade ruido de medición como si fuera un radar real, y luego aplica el filtro para reconstruir la trayectoria verdadera a partir de esas mediciones imperfectas. La mejora respecto a quedarte solo con el dato bruto del sensor ronda el cuarenta o cincuenta por ciento de reducción de error, dependiendo de cuánto ruido metas y cuánto maniobre el blanco. Verlo animado, ping a ping, mientras la línea cian del filtro converge sobre la trayectoria real mientras la línea naranja del radar bruto sigue saltando erráticamente, es de las pocas veces que una ecuación de álgebra matricial se vuelve intuitiva con solo mirarla.

Imagen 3 – KALMAN RADAR TRACKER – Seguimiento de blancos aéreos

Esto me lleva a algo que pienso desde hace tiempo sobre el sector geoespacial y por qué cada vez se parece más a otros sectores que en apariencia no tienen nada que ver. La frontera entre GIS, teledetección, radar de defensa y ciencia de datos se está disolviendo, no porque las aplicaciones converjan, sino porque la base matemática siempre fue la misma y durante años cada comunidad la vistió con su propio vocabulario. Un analista GIS que entiende bien la incertidumbre espacial entiende sin mucho esfuerzo el tracking radar.

Imagen 5 – KALMAN RADAR TRACKER – Exportando a KML todas las trayectorias y datos puntuales

Un ingeniero de radar que entiende el filtrado de señal entiende sin mucho esfuerzo el procesado SAR. La diferencia real no está en la herramienta, está en el dominio de aplicación y en el contexto operativo, que sí importan mucho, pero no son la barrera que parecen desde fuera.

Imagen 6 – KALMAN RADAR TRACKER – Exportando a KML todas las trayectorias y datos puntuales

Esto también explica por qué el geointeligencia y la observación terrestre están viviendo un momento tan interesante ahora mismo. Hay una demanda creciente de perfiles que sepan moverse entre estos mundos, que entiendan tanto la física de la señal como la lógica espacial del problema, y que no se queden bloqueados cuando el vocabulario cambia de “ground truth” a “blanco” o de “pixel” a “celda de resolución radar”. La tecnología emergente en este espacio no va de inventar matemática nueva, va de aplicar matemática ya madura a dominios que históricamente estuvieron separados por silos institucionales más que por silos técnicos.

Voy a seguir publicando sobre esto, porque cuanto más meto las manos en código de tracking y en estadística de procesos espacio-temporales, más confirmo que el geógrafo que sabe programar y el ingeniero que sabe leer un mapa están resolviendo, en el fondo, el mismo tipo de incertidumbre.

Descarga un KML de ejemplo. Renombra el ZIP en KML y ya lo tienes!Descarga

Alberto Concejal
GIS Analyst

Sources:
https://kalmanfilter.net/

Imagen 7 – KALMAN RADAR TRACKER – Output: GIF y KML

VENEZUELA EARTHQUAKE RESPONSE using DuckDB, Overture Maps and R

Just wanted to update on the usage of the tool I developed (OVERTURE MAPS EXTRACTOR) for extraction of Open data from Overture Maps for a quick hands on.

Overture Maps Extractor developed by Alberto Concejal (Open Interface AOI1)

Summary:

Latest release 2026/06/17
3 km buffer over Caracas downtown: 5 minutes
Buildings 52,006 items
Roads 5,482 items
POIS 4,711 items
LULC 398 items
LAND 441 items
Admin Bounds 77 items
Infrastructure 2,442 items

Exported to a 17 MB geopackage GPKG file

Please let me know if this interests anybody, free use of course,

Overture Maps Extractor developed by Alberto Concejal (Open Interface AOI2)

Overture Maps Extractor developed by Alberto Concejal (Global Mapper)

Hugs to all my friends from Venezuela!

Alberto Concejal
Geospatial Analyst

Overture Maps Extractor developed by Alberto Concejal (buildings’ output)

Sources:
https://developmentseed.org/stac-map/?href=https://vantor-opendata.s3.amazonaws.com/events/Venezuela-Earthquake-Jun-2026/B1400011000BDF10.json
https://radiantearth.github.io/stac-browser/#/external/vantor-opendata.s3.amazonaws.com/events/Venezuela-Earthquake-Jun-2026/collection.json?.language=es
https://x.com/vantortech/status/2070496092131569833

Bingham Canyon: El deslizamiento más grande de la historia minera moderna analizado con radar SENTINEL-1

Análisis de cambios con SAR (Radar de Apertura Sintética) usando Sentinel-1 sobre la mina Bingham Canyon en Utah, donde ocurrió uno de los mayores deslizamientos de tierra de la historia minera el 10 de abril de 2013. El análisis no es PRE/POST del evento en sentido estricto. Es una detección de cambios entre dos períodos posteriores al deslizamiento:

*PRE en el script = oct 2014 – jun 2015 (primera referencia disponible)
POST en el script = jul 2015 – mar 2016 (un año después)

Lo que se detecta no es “antes vs después del colapso de 2013”, sino la evolución de la cicatriz entre 2015 y 2016: reconfiguración de taludes, movimiento de material, estabilización o actividad residual de la zona afectada.

Las cicatrices geomorfológicas de un deslizamiento de esa magnitud no desaparecen en meses. La rugosidad anómala, los depósitos de escombros y la geometría alterada del pit siguen siendo detectables por el radar años después. Lo que el análisis captura es la dinámica post-colapso, no el colapso en sí. el radar detecta los cambios superficiales ocurridos entre 2015 y 2016 sobre la zona afectada por el deslizamiento de 2013. Es un matiz importante que señalo de entrada para no generar confusión.

Imagen 1 – Contexto (imagen reciente 2025)

La clave, el backscatter, pero ¿qué es?. El satélite emite un pulso de microondas hacia la Tierra, ese pulso golpea la superficie y se dispersa en todas direcciones. El backscatter es la fracción que regresa exactamente hacia el sensor. Se mide en decibelios (dB), donde valores más negativos = menos energía devuelta.

Qué determina cuando vuelve?. Tres factores principales: 1 Humedad/dieléctrico: suelo húmedo o vegetación densa absorben y devuelven más energía 2 Rugosidad superficial: una roca fragmentada devuelve muchísimo más que una superficie lisa 3 Geometría: esquinas y estructuras verticales crean “corner reflectors” que disparan el backscatter

Imagen 2 – Los datos crudos: así “ve” el radar

Antes del colapso*: la mina tiene taludes estables, geometría conocida, señal SAR consistente.

Después: millones de toneladas de roca fragmentada y removida cambian radicalmente la rugosidad y geometría de la superficie. El radar lo ve como un cambio brusco de backscatter, aunque haya nubes, aunque sea de noche.

Ahí está la clave: el radar no necesita luz solar ni cielo despejado. Ve a través de todo.

Valor Cercano a 0 dB Casi toda la energía vuelve (metal, agua agitada, roca desnuda)

Valor −10 a −15 dBVegetación, suelo moderado

Valor < −20 dBAgua calma, superficies muy lisas (casi nada vuelve)

En este análisis, un Δ de +3 dB o más indica que esa zona devuelve el doble de energía que antes, señal inequívoca de que algo cambió físicamente en la superficie.

Como S1 no se lanzó hasta 2014, el script usa las primeras imágenes disponibles como referencia base en lugar de tener un “antes” real del evento. El flujo tiene tres grandes bloques:

Preparación de datos

Filtra la colección COPERNICUS/S1_GRD en modo Interferometric Wide (IW), manteniendo solo imágenes con ambas polarizaciones VV y VH. Define dos ventanas temporales: PRE (oct 2014 – jun 2015) y POST (jul 2015 – mar 2016). De cada ventana extrae una mediana, que por sí sola ya reduce bastante el speckle. Encima aplica un filtro boxcar 3×3 adicional por banda.

Imagen 4 – La mina en 3D usando datos LIDAR del USGS de 2019

Detección de cambios

Calcula la diferencia POST − PRE en decibelios para VV, VH y su media. Un Δ positivo indica aumento de backscatter (material acumulado, mayor rugosidad superficial), un Δ negativo indica pérdida o remoción. Aplica un umbral de ±3 dB para distinguir cambio significativo del ruido de fondo.

Salidas y análisis

Visualiza cuatro capas SAR (PRE/POST × VV/VH), los tres mapas de diferencia y un RGB multitemporal donde el canal R lleva POST-VV, G lleva PRE-VV y B lleva la diferencia, produciendo tonos magenta donde hay ganancia y verdes donde hay pérdida. Calcula el área afectada en hectáreas con pixelArea y genera una serie temporal de backscatter medio en el AOI entre 2014 y 2016. Finalmente exporta a Drive los rásteres de diferencia VV, VH y la máscara de cambio binarizada.

Imagen 5 – RGB Multitemporal (R=POST, G=PRE, B=Diff

Imagen 6 – El cambio, clasificado ¡a que mola!

Diagrama 1 – Sumarios de ganancias y pérdidas

Un detalle metodológico relevante: el script no fuerza una órbita concreta (ascending o descending) porque la mediana temporal compensa la mezcla de geometrías de adquisición, aunque para un análisis más riguroso lo ideal sería separar órbitas.

Conclusión

El deslizamiento de Bingham Canyon ocurrió en abril de 2013. Sentinel-1 no existía todavía. Y aun así, el radar fue capaz de leer las cicatrices que dejó en el terreno más de un año después.

Eso dice algo importante sobre la naturaleza del SAR: no necesita estar en el momento exacto para detectar que algo cambió. La geometría del terreno, la rugosidad de la roca fragmentada, la reconfiguración de los taludes… todo queda impreso en el backscatter durante meses, incluso años.

Lo que este script demuestra no es solo que GEE puede procesar imágenes SAR en la nube sin descargar un solo píxel. Demuestra que con datos abiertos, código reproducible y una ventana temporal bien elegida, es posible reconstruir la huella espacial de un evento catastrófico desde cero.

El siguiente paso natural es la tercera dimensión. El radar nos dice dónde cambió la superficie. El LiDAR USGS 3DEP (Imagen 4) nos dice cuánto volumen se desplazó. La combinación de ambas fuentes: SAR multitemporal + MDT de alta resolución, es exactamente el tipo de análisis que los equipos de gestión de riesgos geológicos necesitan, y que hoy es accesible sin infraestructura dedicada, solo conociendo las fuentes.

Imagen 7 – Usando LIDAR de alta resolución para medir ganancias y pérdidas de la zona de análisis

Espero que te haya gustado!

Alberto Concejal
GIS Analyst
https://storymaps.arcgis.com/stories/a5398ca1f6ea4b05b90272e38281c903
https://apps.nationalmap.gov/lidar-explorer

SolarScope: cuando el catastro, el LiDAR y el sol se sientan a la misma mesa

Llevo unos días dándole vueltas a una idea que, en el fondo, es bastante sencilla: si tenemos la huella de cada edificio, su altura y un modelo digital de superficies de alta resolución, ¿por qué seguimos viendo estudios de potencial solar que tratan los tejados como manchas homogéneas sobre un mapa? De esa pregunta, y de unas cuantas sesiones intensas de R, ha salido SolarScope, una aplicación Shiny que estoy desarrollando para hacer scoring de potencial fotovoltaico tejado a tejado, con datos abiertos y un flujo que se puede reproducir tanto en España como en cualquier sitio del mundo donde pueda conectarme a un DSM de alta resolución!

Imagen 1 – Solar Scope, midiendo el potencial fotovoltaico de TODOS los edificios de España

La motivación es doble. Por un lado, profesional: vengo de quince años moviéndome entre geografía, teledetección y GIS aplicado, y cada vez que he tocado proyectos de energía solar he visto el mismo cuello de botella. Los modelos de potencial suelen apoyarse en rásteres de baja resolución (SRTM, Copernicus DEM a 30 m) que son perfectamente válidos para planificación territorial a gran escala, pero que se quedan cortos en cuanto entras en el detalle de una nave industrial o un polígono residencial: ahí lo que importa es la sombra que proyecta el edificio de al lado, la orientación real de la cubierta y cuántos metros cuadrados aprovechables tiene cada tejado después de descontar lucernarios, antenas y pasillos técnicos. Por otro lado, hay una motivación más simple: tenía ganas de construir algo vistoso, con mapas 3D, que sirviera como pieza de presentación ante empresas del sector —y de paso, demostrar que con herramientas open source se puede llegar muy lejos sin depender de licencias de software propietario.

Cómo funciona SolarScope (vídeo YouTube)

De dónde vienen los datos

El corazón de SolarScope es la combinación de tres fuentes. La primera es Overture Maps, el proyecto colaborativo (Meta, Microsoft, Amazon, TomTom y otros) que publica footprints de edificios a escala global con atributos de altura y número de plantas, distribuidos como Parquet sobre S3. Aquí es donde DuckDB se convierte en la pieza más elegante del stack: con la extensión httpfs y spatial, puedo lanzar una consulta SQL directamente contra los ficheros remotos, filtrar por bounding box usando los campos bbox.xmin/xmax/ymin/ymax y traerme solo los edificios que caen dentro del área de interés, sin descargar nada de más. Cuando la altura no está disponible —que pasa más a menudo de lo que gustaría— caigo en una estimación a partir del número de plantas, y si tampoco hay eso, asumo una nave de una planta. No es perfecto, pero es razonable para naves industriales, que es justo el tipo de cubierta que más interesa a un desarrollador solar.

https://albertogis.shinyapps.io/SolarScope

Imagen 2 – Solar Scope, midiendo el potencial fotovoltaico de TODOS los edificios de España

La segunda fuente es el Modelo Digital de Superficies (MDS) del vuelo LiDAR del PNOA, servido por el CNIG a través de un servicio WCS. Aquí sí que tuve que arremangarme: el endpoint correcto no es el de siempre (el del MDT del IGN), sino wcs-mds.idee.es/mds, y el identificador de cobertura tampoco sigue la convención que uno esperaría —nada de Elevacion25830_5, sino simplemente mds05—. Una vez resuelto eso (y la consabida reproyección de EPSG:4326 a EPSG:25830, porque el servicio trabaja en metros UTM, no en grados), el DSM permite calcular, para cada edificio, si hay construcciones vecinas más altas que le proyecten sombra, y con eso derivar un factor de sombreado real en lugar de uno inventado. Para Estados Unidos, el equivalente es USGS 3DEP, que ofrece DSM de hasta 1 metro de resolución allá donde hay cobertura LiDAR —el conector está escrito y solo pendiente de pruebas con datos reales, pero la arquitectura ya contempla ambos países sin tocar el resto de la aplicación.

La tercera pieza es, simplemente, la geometría de cada footprint: a partir de la relación de aspecto del bounding box estimo un factor de orientación, asumiendo que la mayoría de cubiertas industriales son planas (donde la orientación pesa poco) pero penalizando ligeramente las naves muy alargadas en sentido norte-sur frente a las que se extienden este-oeste, que en el hemisferio norte “miran” más al sur.

El cálculo: del tejado al kWp

Con esos tres ingredientes —área, altura, orientación y sombreado—, el motor de scoring aplica un modelo bastante directo inspirado en los valores de referencia de PVGIS: irradiación global horizontal anual (en torno a 1.650 kWh/m²/año para Madrid), un porcentaje de área útil tras descontar elementos técnicos, una eficiencia de sistema fotovoltaico del 20% y una densidad de potencia instalable de 0,18 kWp por metro cuadrado útil. El resultado es, para cada edificio, una potencia instalable en kWp, una producción anual estimada en MWh y un Solar POI Score normalizado de 0 a 100 que combina producción total con superficie disponible —porque un tejado grande y mediocre puede ser más interesante para un desarrollador que uno pequeño y perfecto—. Los tejados se clasifican en tres niveles (bajo, medio, alto), lo que permite, de un vistazo, identificar qué activos merecen una visita técnica.

Imagen 3 – Solar Scope. Interfaz apaisada, búsqueda de usabilidad al máximo

La aplicación: de la idea al mapa interactivo

Todo esto vive dentro de una app Shiny con un diseño oscuro deliberadamente “de producto”, construido sobre bslib con tipografías Inter y Space Grotesk, y una paleta que va del morado profundo al amarillo solar —la misma que reaparece en cada gráfico, en la leyenda del mapa y, por qué no, en la portada de la propia aplicación, donde un skyline de edificios generado en SVG (con sus correspondientes “sombreros” de colores en cada tejado) hace de fondo desenfocado tras el logo.

Imagen 4 – Solar Scope. Vista en 3D. Dentro de poco meteré los nuevos Landmarks de MAPBOX!!!!!

La interacción es deliberadamente simple: el usuario hace clic en el mapa para fijar un punto central, ajusta con un slider el radio de análisis —entre 100 y 1.000 metros— y pulsa “Analizar zona”. En ese momento, la app reconstruye el bounding box, relanza las consultas a Overture y al MDS del CNIG, recalcula todas las métricas y redibuja tanto el mapa 2D (con CartoDB Dark Matter como base) como una vista 3D con extrusión de edificios por altura y coloreada por score, esta última construida con mapdeck, que en el fondo es un envoltorio de deck.gl sobre Mapbox GL. Si el conector de datos reales falla por cualquier motivo —sin conexión, cambio de release de Overture, servicio del CNIG caído—, la app cae automáticamente a un generador de datos sintéticos con la misma estructura, y lo indica con un discreto badge de color para que nunca te quedes con una pantalla en blanco en mitad de una demo. Para quien necesite explotar los resultados fuera de la app, hay un botón de exportación a GeoPackage con todos los atributos calculados, y otro que genera un informe ejecutivo en HTML autocontenido —mapa, KPIs y ranking incluidos— listo para enviar por correo.

Como guiño visual adicional, y porque a veces lo más vistoso no tiene por qué ser lo más complejo, incorporé también una capa opcional de imagen de alta resolución (Esri World Imagery) con control de opacidad sobre el mapa base oscuro: sirve únicamente para contextualizar visualmente la zona —el cálculo de sombreado sigue apoyándose en el DSM del CNIG—, pero al venir servida como teselas por CDN carga muchísimo más rápido que el WMS del PNOA, que para un uso puramente visual resultaba innecesariamente pesado.

Imagen 5 – Solar Scope. Acceso pormenorizado a todos los edificios del AOI

El stack, en una frase

Si tuviera que resumir la pila tecnológica en un párrafo: R y Shiny como columna vertebral, bslib y CSS personalizado para la interfaz, sf y terra para todo lo geoespacial vectorial y ráster, DuckDB con httpfs/spatial como motor de consulta sobre datos en la nube sin descargas intermedias, leaflet para el mapa 2D y mapdeck para la vista 3D, DT para las tablas interactivas y rmarkdown para los informes ejecutivos. Todo open source, todo reproducible, y todo pensado para que el mismo esqueleto sirva tanto para un polígono industrial en Vicálvaro como, cambiando el conector de DSM, para un parque empresarial en Arizona.

Para quién es esto

El público natural son desarrolladores y operadores de energía solar que necesitan priorizar carteras de tejados —ya sea para autoconsumo industrial, comunidades energéticas o grandes cubiertas logísticas— sin tener que encargar un estudio LiDAR específico cada vez que aparece una oportunidad. También tiene sentido para consultoras de sostenibilidad que necesiten estimar potencial fotovoltaico como parte de informes ESG, para administraciones locales que quieran mapear el potencial solar de sus polígonos industriales, o simplemente para cualquier estudio de GIS que quiera mostrar que el análisis espacial de alta resolución no tiene por qué vivir solo dentro de un escritorio ArcGIS. Y sí, en estos días concretos la estoy preparando como pieza de demostración para una conversación con una empresa del sector solar —si sale adelante, ya contaré más por aquí.

Imagen 6 – Solar Scope, midiendo el potencial fotovoltaico de TODOS los edificios de España

Como suele pasar en este oficio, la parte más “glamurosa” —el mapa 3D girando con la cámara, los colores del score, la portada con el skyline desenfocado— es la que menos tiempo me ha llevado. Lo que de verdad ha consumido las horas ha sido, cómo no, encontrar el COVERAGEID correcto de un servicio WCS y pelearme con la codificación de caracteres en una consola de R en Windows: si alguna vez os preguntáis por qué los geógrafos envejecemos mal, no es por el sol de las salidas de campo, es por los acentos UTF-8 dentro de backticks de R. Reproject responsibly.

Imagen 7 – Solar Scope. Exportando geometrías a otro software, en este caso Global Mapper

Iré documentando en este blog los siguientes pasos: cerrar el conector de Estados Unidos con USGS 3DEP, refinar la orientación de cubierta con geometría de fachadas reales, y —si el tiempo lo permite— un modelo de sombreado hora a hora basado en la posición solar real a lo largo del año.

Si todo va bien, la próxima entrada debería poder escribirse en mucho menos tiempo que esta: la idea es que cada iteración de SolarScope deje el terreno un poco más allanado para la siguiente.

Espero que os haya interesado!

Alberto Concejal
Geospatial Analyst

Time Series Analyzer: Análisis multitemporal de índices espectrales en Google Earth Engine

He desarrollado una aplicación interactiva en Google Earth Engine para la extracción y análisis estadístico automático de series temporales de cinco índices espectrales (NDVI, EVI, SAVI, NDWI y NBR) sobre cualquier geometría definida por el usuario en cualquier sitio del mundo. El objetivo es pasar de una imagen satélite puntual a una comprensión temporal del territorio: qué ha pasado, qué patrón subyace, y qué cabe esperar.

Imagen 1 – Time Series analyzer, the interface

La app está escrita íntegramente en JavaScript usando la GEE UI API, sin dependencias externas. El procesamiento corre en la nube de GEE sobre las colecciones oficiales de Landsat (L5/L7/L8/L9, Collection 2 SR, desde 1984) y Sentinel-2 SR Harmonized (desde 2017). El cálculo de la ACF se realiza en el cliente una vez extraída la serie como FeatureCollection. La interfaz se organiza en tres paneles: controles a la izquierda (selección del índice, plataforma, fechas, nubosidad límite), mapa y serie temporal en el centro, y análisis estadístico a la derecha con sistema de pestañas para separar gráfico, resultado e interpretación. Aquí abajo puedes acceder a la aplicación:

https://analysis202601.projects.earthengine.app/view/timeseriesanalyzer

Imagen 2 – the original index data over the selected dates

Metodología

El flujo parte del filtrado de imágenes por fecha, geometría y porcentaje de nubosidad máximo, con enmascarado de nubes por QA_PIXEL (Landsat) o QA60 (Sentinel-2). Sobre esa colección se calculan los índices espectrales aplicando los factores de escala de Collection 2 sobre reflectancia de superficie. Los cinco índices disponibles son NDVI, EVI, SAVI, NDWI y NBR, cada uno orientado a una pregunta diferente sobre la cobertura: estado general de la vegetación, respuesta en zonas densas, corrección del efecto suelo, contenido hídrico o severidad de incendios.

La serie temporal resultante se analiza estadísticamente mediante cuatro módulos. El más relevante metodológicamente es la función de autocorrelación (ACF), que mide la correlación de la serie consigo misma en distintos desfases temporales. Un pico significativo en lag 1 indica que la vegetación tiene inercia: lo que ocurre hoy predice lo que ocurrirá en la siguiente adquisición. Un pico en lag 6 o lag 23 —dependiendo de la cadencia del satélite— revela un ciclo semestral o anual. El umbral de significancia se calcula según el criterio de Bartlett al 95% (±1.96/√n), lo que permite distinguir entre estructura real y ruido estadístico.

Imagen 3 – Función de Autocorrelación y sus intervalos de confianza

La interpretación de la ACF se estructura automáticamente en cuatro bloques: memoria a corto plazo, ciclo estacional con estimación del período en meses, longitud de memoria total y recomendación de modelo. Si la autocorrelación permanece significativa hasta lags altos, la serie tiene memoria larga y probablemente una tendencia no estacionaria que requiere diferenciación antes de modelar. Si hay un pico estacional claro, la estructura es compatible con un modelo SARIMA(p,d,q)(P,D,Q)[s] donde s es el período detectado. La app complementa esto con un análisis de tendencia lineal que aporta la pendiente y el R², y una curva suavizada que permite comparar la amplitud del ciclo estacional entre años.

Comprensión de patrones y proyección

Lo que hace útil el análisis temporal no es la fotografía de un año sino la estructura que emerge a lo largo de varios ciclos. La ACF permite responder preguntas concretas: ¿es el ciclo estacional de este pastizal estable año a año o su amplitud está decreciendo? ¿La caída de NDVI en 2022 fue un evento puntual o marcó un cambio de régimen? ¿Cuántos meses tarda la vegetación en recuperar sus valores tras un incendio?

Identificada la estructura temporal, el paso natural será el modelado predictivo (para un próximo post si me lo pedís:-)). Una serie con estacionalidad anual clara y memoria moderada es directamente modelable con SARIMA, cuyo ajuste óptimo puede automatizarse con auto.arima() en R o SARIMAX en Python sobre el CSV exportado por la app. El modelo resultante permite proyectar valores futuros del índice, construir intervalos de confianza estacionales y detectar anomalías como desviaciones significativas respecto al patrón esperado, lo que equivale a una señal de alerta temprana ante degradación, sequía o cambio de uso del suelo.

Aplicaciones

La combinación de Landsat desde 1984 con series temporales largas abre un abanico amplio de usos. En el contexto de pastizales y agricultura, permite cuantificar la respuesta de la cobertura vegetal a años secos o húmedos y comparar la dinámica entre parcelas con distinto manejo. En el ámbito forestal, la caracterización fenológica mediante el ciclo anual del NDVI o EVI permite detectar adelantos o retrasos en la brotación asociados a cambios climáticos. Con NBR, el seguimiento post-incendio revela tanto la severidad inicial como la trayectoria de recuperación. En zonas de expansión urbana, la tendencia negativa sostenida en NDVI a lo largo de décadas es un indicador directo de sellado del suelo y pérdida de cobertura vegetal.

En investigación, la exportación CSV permite integrar los resultados en flujos de trabajo más complejos: correlación con datos meteorológicos, validación con trabajo de campo o calibración de modelos de simulación de ecosistemas. En docencia, la app funciona como entorno de demostración interactivo para explicar autocorrelación, estacionalidad y modelado de series temporales sobre datos reales sin necesidad de instalar nada.

Aquí va el resultado del análisis de autocorrelation en una zona muy cercana a Zaragoza, España, usando Landsat, un periodo que comprendía los tres últimos años.

Imagen 6 – Análisis extra cerca de Ejea de los Caballeros en Zaragoza, España

Imagen 7 – Función de Autocorrelación y sus intervalos de confianza

WHAT THIS ANALYSIS DOES

The ACF measures how similar the series is to a time-shifted version of itself. Each bar is a lag: lag 1 = correlation between observation i and the next; lag 2 = two steps back, etc. Dashed lines mark the 95% significance threshold (±1.96/√n). Bars exceeding that threshold indicate a real pattern, not noise.

AUTOMATIC RESULT

1. SHORT-TERM MEMORY (lag 1 = 0.89) ► Significant (very strong). Consecutive values are correlated. A high index today predicts a high index in the next image.

2. SEASONAL CYCLE (strongest lag = 2, r = 0.75) ► Significant peak → clear periodicity detected. Landsat 8-9 cadence: cycle = 8-16 days. Sentinel cadence: cycle = 3-5 days. This is consistent with a short intra-seasonal fluctuation.

3. MEMORY LENGTH (last significant lag = 24) ► Long memory — ACF remains significant up to lag 24. This suggests a non-stationary trend in the data. Recommendation: apply first-order differencing before fitting ARIMA/SARIMA.

4. MODEL SUGGESTION ► Non-stationary series. Difference first, then fit ARIMA(p,1,q).

READING GUIDE

Reading guide: → High lag 1 (>0.5): vegetation has strong “memory”, changes slowly. → Peak at lag 6 or 12 (~16-day data): semi-annual or annual cycle — typical grassland seasonality. → ACF that decays slowly into a “tail”: possible non-stationary trend in the series. → ACF oscillating between positive and negative: clear cyclic pattern, ideal for SARIMA modeling. → All bars within the lines: series with no temporal structure — possibly noise or heterogeneous geometry.

Espero que os haya resultado interesante,

Alberto Concejal
https://geovisualization.net/

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.

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.

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.

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