GeoPandas para Dummies

Librerías y Entorno de Trabajo

Para el desarrollo de esta sesión utilizaremos GeoPandas, una librería de código abierto fundamental en el ecosistema de Python geoespacial. GeoPandas extiende las capacidades de pandas para permitir operaciones espaciales en tipos geométricos. Funciona integrando otras librerías esenciales: utiliza pandas para el manejo de datos tabulares, shapely para realizar las operaciones geométricas y fiona para el acceso y lectura de archivos.

Para comenzar, importamos la librería utilizando el alias estándar que nos permitirá acceder a sus funciones de manera abreviada.

import geopandas as gpd

Lectura de Fuentes de Datos Espaciales

En el análisis geográfico es común trabajar con diversos formatos de archivos. GeoPandas simplifica este proceso mediante el comando unificado read_file, el cual detecta automáticamente el formato del archivo de entrada. Trabajaremos con los formatos más estandarizados en la industria.

El formato Shapefile (.shp) es el estándar más antiguo y común. Aunque lo vemos como un solo archivo, en realidad es un conjunto de archivos que contienen la geometría (.shp), los atributos (.dbf) y la proyección (.prj). También es posible leer archivos comprimidos en formato .zip directamente sin necesidad de descomprimirlos previamente, lo cual optimiza el flujo de trabajo cuando descargamos datos de internet.

Adicionalmente, manejamos formatos modernos como GML o GeoJSON, muy utilizados en aplicaciones web y servicios de mapas.

# Lectura de un Shapefile tradicional
localidades = gpd.read_file(r'ruta/a/tu/archivo/localidades.shp')

# Lectura directa de un archivo comprimido
barrios = gpd.read_file(r'ruta/a/tu/archivo/barrios.zip')

# Lectura de formatos web (GML y GeoJSON)
centroides = gpd.read_file('ruta/archivo.gml')
datos_json = gpd.read_file('ruta/archivo.geojson')

Estructura del GeoDataFrame y Manejo de Índices

Una vez cargados los datos, obtenemos un GeoDataFrame. Estructuralmente es idéntico a un DataFrame de pandas, pero incluye una columna especial llamada geometry. Esta columna es vital porque almacena la información espacial (puntos, líneas o polígonos) y es sobre la cual se ejecutarán los cálculos geográficos.

Para facilitar el análisis, es una buena práctica redefinir el índice de la tabla. En lugar de usar números consecutivos autogenerados (0, 1, 2...), asignamos una columna con valor semántico, como el nombre de la localidad, para poder realizar búsquedas y filtros de manera más intuitiva.

# Asignar una columna de texto como índice principal
localidades_shp = localidades.set_index('Nombre')

# Exploración rápida de los primeros y últimos registros
localidades_shp.head(3)
localidades_shp.tail(4)

Atributos Geométricos Básicos: Área y Perímetro

Los objetos espaciales tienen propiedades intrínsecas. Como estamos trabajando con un sistema proyectado, podemos calcular métricas euclidianas directamente. La propiedad .area nos devuelve la superficie de cada polígono en las unidades del sistema de referencia (usualmente metros cuadrados), mientras que .boundary extrae el límite del polígono, convirtiéndolo en una línea (LineString), lo que nos sirve para calcular perímetros.

Es importante notar que estas operaciones se realizan 'fila por fila' de manera vectorizada, lo que hace el procesamiento muy eficiente.

# Cálculo del área
localidades_shp['area'] = localidades_shp.area

# Extracción del límite (perímetro)
localidades_shp['perimetro'] = localidades_shp.boundary

Geometría Constructiva: Centroides y Envolventes

Más allá de medir, podemos construir nuevas geometrías derivadas de las originales. Un concepto clave es el Centroide, que representa el centro geométrico o punto de equilibrio de un polígono. Esto es útil cuando queremos simplificar un mapa de polígonos a puntos para análisis de distribución.

Otro concepto teórico relevante es el Convex Hull (envolvente convexa). Imagina que colocas una banda elástica alrededor de los puntos que forman tu geometría; la forma resultante es el polígono convexo más pequeño que contiene a todos los puntos. Esto se usa para simplificar formas complejas o analizar dispersión.

# Generación de centroides (Puntos)
localidades_shp['centroide'] = localidades_shp.centroid

# Generación de la envolvente convexa (Polígono simplificado)
localidades_shp['envolvente'] = localidades_shp.convex_hull

Análisis de Proximidad: Distancias y Buffers

El análisis espacial muchas veces se trata de entender la relación entre objetos. Para calcular distancias, necesitamos un punto de referencia. La función .distance() calcula la distancia euclidiana más corta entre geometrías.

El concepto de Buffer (área de influencia) es fundamental en SIG. Consiste en crear un polígono que abarca toda el área dentro de una distancia específica desde una geometría dada. Esto es esencial para estudios de impacto ambiental o cobertura de servicios.

# Cálculo de distancia desde un punto fijo a todas las localidades
punto_inicial = localidades_shp['centroide'].iloc[0]
localidades_shp['distancia'] = localidades_shp['centroide'].distance(punto_inicial)

# Promedio de distancias calculadas
promedio_dist = localidades_shp['distancia'].mean()

# Creación de un buffer de 10,000 unidades
localidades_shp['buffer_10k'] = localidades_shp.buffer(10000)

Visualización y Cartografía

La visualización es la herramienta principal para la exploración de datos. GeoPandas utiliza matplotlib como motor gráfico. Podemos crear mapas de coropletas, donde el color de cada polígono representa el valor de una variable (como el área), permitiendo identificar patrones visuales rápidamente.

Para mapas más complejos que requieren superponer varias capas (por ejemplo, mostrar los polígonos de las localidades y encima sus centroides), utilizamos el concepto de axis (ejes). Primero dibujamos la capa base y guardamos su referencia en una variable ax, y luego pasamos esa referencia a las siguientes capas para que se dibujen en el mismo lienzo.

# Mapa básico coloreado por área con leyenda
localidades_shp.plot(column='area', legend=True)

# Mapa cambiando la geometría activa a centroides
localidades_shp.set_geometry('centroide').plot()

# Superposición de capas (Polígonos + Centroides)
ax = localidades_shp.plot(color='black')
localidades_shp['centroide'].plot(ax=ax, color='red')

# Mapa personalizado con tamaño de figura
localidades_shp.plot(
    column='Codigo',
    legend=True,
    figsize=(12, 20),
    legend_kwds={'loc': 'upper left'}
)

Operaciones de Conjuntos y Consultas Espaciales

Al igual que en SQL o pandas, podemos filtrar datos, pero también podemos realizar preguntas espaciales. La función .intersects() es una operación booleana que verifica si dos geometrías comparten algún espacio en común (se tocan o se solapan). Esto permite responder preguntas como '¿Qué áreas de influencia tocan una localidad específica?'.

También utilizamos .loc para selecciones por etiqueta y .query para consultas basadas en cadenas de texto, lo cual hace el código más legible.

# Selección de una localidad específica
engativa = localidades_shp.loc['Engativá']

# Selección mediante consulta SQL-like
localidad_12 = localidades_shp.query("Codigo_Loc == '12'")

# Verificación de intersección espacial
intersecciones = localidades_shp['buffer_10k'].intersects(engativa.geometry)

Sistemas de Referencia de Coordenadas (CRS)

Finalmente, un concepto crítico en cualquier análisis SIG es el CRS (Coordinate Reference System). Los datos geográficos están definidos sobre la superficie curva de la Tierra, pero para medirlos y dibujarlos en una pantalla plana, necesitamos proyectarlos.

Cada proyección tiene un código único llamado EPSG. Es vital verificar en qué sistema están nuestros datos con .crs. Si necesitamos integrar datos de diferentes fuentes o realizar cálculos en diferentes unidades (grados vs metros), debemos reproyectar los datos usando .to_crs.

# Consultar el sistema de referencia actual
print(localidades_shp.crs)

# Reproyectar a coordenadas geográficas mundiales (WGS84 - EPSG:4326)
localidades_wgs84 = localidades_shp.to_crs(epsg=4326)

Esta guía de estudio detalla el flujo de trabajo para realizar un análisis de densidad poblacional utilizando la librería GeoPandas en Python, tomando como caso de estudio la ciudad de Bogotá. En esta clase se explicará cómo procesar datos espaciales desde su carga hasta su visualización final, integrando conceptos de unión espacial y estadística descriptiva.

Configuración del Entorno e Importación de Librerías

El primer paso fundamental en este flujo de trabajo es la configuración del entorno y la importación de las librerías necesarias. Se utilizará GeoPandas, una extensión de Pandas diseñada para trabajar con datos geoespaciales. La convención estándar para importar esta librería es mediante el siguiente comando. Adicionalmente, para realizar operaciones de cruce espacial, se requerirá importar una herramienta específica contenida dentro de las utilidades de GeoPandas. Es esencial contar con un entorno previamente configurado, como un ambiente en Anaconda, donde GeoPandas ya se encuentre instalado.

import geopandas as gpd
from geopandas.tools import sjoin

Carga de Datos Espaciales

Una vez preparado el entorno, se procederá a la carga de los datos espaciales. En este ejercicio se utilizarán dos conjuntos de datos principales en formato GeoJSON. El primero corresponde a una "grilla" o cuadrícula hexagonal que divide la zona de estudio, la cual posee una geometría de tipo polígono. El segundo conjunto de datos contiene los "centroides de manz anas", que son geometrías de tipo punto representando el centro geométrico de las manzanas catastrales, y que contienen la información demográfica clave, como la densidad poblacional.

Para cargar estos archivos se utiliza el método .read_file(), asignando la ruta del archivo correspondiente y almacenando el resultado en una variable que actuará como un GeoDataFrame. Es recomendable incluir el prefijo r antes de la ruta del archivo para evitar conflictos con caracteres especiales en la dirección.

# Carga de la grilla hexagonal
grilla = gpd.read_file(r'ruta/a/tu/archivo/grilla.geojson')

# Carga de los centroides de manzanas
centroides = gpd.read_file(r'ruta/a/tu/archivo/centroides_manzanas.geojson')

Exploración Visual de los Datos

La exploración visual de los datos es crucial antes de cualquier análisis. Se utilizará el método .plot() sobre cada GeoDataFrame para inspeccionar la geometría. Para mejorar la visualización, se puede ajustar el tamaño de la figura utilizando el parámetro figsize=(ancho, alto). Esto permite verificar que la grilla efectivamente cubre el área de interés y que los puntos de los centroides se distribuyen espacialmente de manera coherente sobre dicha grilla.

# Visualización de la grilla
grilla.plot(figsize=(20, 12))

# Visualización de los centroides
centroides.plot(figsize=(20, 12))

Unión Espacial (Spatial Join)

El núcleo del análisis reside en la operación de unión espacial o "Spatial Join". Dado que la información de densidad está en los puntos (manzanas) y la estructura de análisis es la grilla (hexágonos), se deben relacionar ambas capas basándose en su ubicación espacial.

La función sjoin toma dos GeoDataFrames y transfiere los atributos de uno al otro basándose en una relación de contenencia espacial (si un punto está dentro de un polígono). El comando genera un nuevo GeoDataFrame que combina la información de ambas fuentes. Esto permite saber a qué hexágono de la grilla pertenece cada manzana censal.

# Unión espacial entre centroides y grilla
union_espacial = sjoin(centroides, grilla)

Limpieza y Filtrado de Datos

Posteriormente, se realizará una limpieza y filtrado de datos mediante el método .query(). Es común encontrar registros con valores nulos o cero que pueden sesgar el análisis. Se aplicará un filtro para conservar únicamente los registros donde la densidad poblacional sea mayor a cero, utilizando una sintaxis similar a SQL. Esto reduce el tamaño del conjunto de datos y enfoca el análisis en zonas habitadas.

# Filtrar registros con densidad mayor a cero
datos_filtrados = union_espacial.query('Densidad_P > 0')

Agrupación y Estadísticas por Zona

Para obtener estadísticas por zona, se introducirá el concepto de agrupación con el método .groupby(). Dado que múltiples manzanas pueden caer dentro de un mismo hexágono de la grilla, es necesario agrupar los datos por el identificador único de la grilla (ID) y sumar los valores de densidad.

La instrucción se estructura seleccionando la columna de interés y aplicando una función de agregación. Esto genera una serie de datos que muestra la densidad total acumulada por cada celda de la grilla, permitiendo identificar, mediante ordenamiento, cuáles son las zonas con mayor concentración de población.

# Agrupar por ID de grilla y sumar densidad
densidad_por_zona = datos_filtrados.groupby('id')['Densidad_P'].sum()

# Ver las zonas con mayor densidad
densidad_por_zona.sort_values(ascending=False).head(10)

Generación del Mapa Temático (Coropletas)

Finalmente, se procederá a la generación del mapa temático o de coropletas. Se utilizará nuevamente el método .plot() sobre el GeoDataFrame resultante de la unión espacial. Para representar la densidad, se asignará la columna de densidad al parámetro column. Para facilitar la interpretación visual de la distribución de datos, se aplicará un esquema de clasificación basado en cuantiles mediante el parámetro scheme='quantiles' y se definirá el número de clases con k. Adicionalmente, se activará la leyenda con legend=True.

El resultado será un mapa donde los colores más intensos indicarán las zonas de mayor densidad poblacional, que en el caso de Bogotá se observan generalmente hacia la periferia, mientras que las zonas centrales tienden a presentar densidades residenciales más bajas.

# Mapa de coropletas con cuantiles
union_espacial.plot(
    column='Densidad_P',
    scheme='quantiles',
    k=5,
    legend=True,
    figsize=(20, 12)
)

En esta clase se aborda el uso de la librería GeoPandas para realizar un análisis de datos abiertos, específicamente enfocado en el reporte de hurtos en los municipios de Colombia. El objetivo principal es conectar una fuente de datos tabular extraída de una API con información espacial para generar visualizaciones cartográficas.

Obtención y carga de datos desde fuentes externas

Para iniciar el ejercicio, es necesario importar la librería principal bajo el alias estándar. A diferencia de ejercicios anteriores donde se cargan archivos locales, en esta sesión se obtienen los datos directamente desde el portal de Datos Abiertos de Colombia (datos.gov.co). Se selecciona el conjunto de datos de reportes de hurtos y se utiliza la URL proporcionada por la API del portal (formato GeoJSON o JSON) para leer la información directamente en el entorno de Python. El comando utilizado para esta lectura remota es gpd.read_file(url), asignando el resultado a una variable que contendrá el GeoDataFrame.

import geopandas as gpd

# Cargar datos desde la API de Datos Abiertos
url = 'https://www.datos.gov.co/api/...'
hurtos = gpd.read_file(url)

Es importante notar que, aunque se utiliza GeoPandas, los datos iniciales de hurtos carecen de una columna de geometría activa, presentándose como una tabla alfanumérica con información como departamento, municipio, tipo de hurto y, crucialmente, el código DANE.

Preprocesamiento y Agrupación de Datos (Groupby)

Una vez cargados los datos, se procede a cambiar el índice del DataFrame utilizando el código DANE mediante el método set_index. Posteriormente, se identifica un problema conceptual para el mapeo: los datos consisten en registros individuales de delitos (filas repetidas por municipio), pero para el mapa se requiere el total de hurtos por entidad territorial.

Para solucionar esto, se explica el concepto de agrupación con el método groupby. Se agrupan los datos por la columna del código DANE y se aplica una función de conteo (count) para obtener un nuevo DataFrame que resume la cantidad total de incidentes por municipio.

# Establecer el código DANE como índice
hurtos = hurtos.set_index('codigo_dane')

# Agrupar por municipio y contar hurtos
hurtos_agrupados = hurtos.groupby('codigo_dane')['id'].count()
hurtos_agrupados = hurtos_agrupados.reset_index(name='count')

Vinculación de Datos Espaciales (Merge vs Spatial Join)

Para dotar de componente espacial a la tabla de estadísticas de hurtos, se carga un segundo archivo que corresponde a la capa vectorial (Shapefile) de los municipios de Colombia, la cual sí contiene la geometría (polígonos).

En este punto se hace una distinción teórica fundamental: dado que la tabla de hurtos no tiene geometría, no es posible realizar una unión espacial (sjoin). En su lugar, se debe realizar una unión atributiva o alfanumérica utilizando el método merge. Este proceso vincula ambas tablas utilizando una columna común, que es el código DANE.

# Cargar la capa de municipios con geometría
municipios = gpd.read_file(r'ruta/municipios.shp')

# Unión alfanumérica por código DANE
resultado = municipios.merge(
    hurtos_agrupados,
    left_on='MPIO_CDPMP',
    right_on='codigo_dane'
)

Análisis de Valores y Ordenamiento

Con los datos ya unificados, se procede a identificar cuáles son los municipios con mayores índices de delincuencia. Para ello se utiliza el método sort_values, ordenando la información basándose en la columna de conteo generada en la agrupación. Es necesario configurar el parámetro ascending=False para visualizar los datos de mayor a menor.

# Ordenar de mayor a menor
resultado.sort_values('count', ascending=False).head(10)

Este análisis preliminar permite identificar cuantitativamente las ciudades con mayor problemática (Bogotá, Cali, Medellín) antes de pasar a la representación gráfica.

Visualización Cartográfica y Simbología

Finalmente, se utiliza el método .plot() para generar el mapa temático. Debido a que la unión (merge) conserva solo los municipios que tienen registros de hurtos, se decide graficar dos capas: una capa base con todos los municipios de Colombia (para dar contexto geográfico) y, sobre ella, la capa de los municipios con hurtos.

Para mejorar la interpretación, se aplica una simbología clasificada basada en la columna de conteo. Se explora el uso de esquemas de clasificación como quantiles definiendo el número de clases con el parámetro k.

# Capa base con todos los municipios
ax = municipios.plot(color='lightgray', edgecolor='white')

# Capa temática con hurtos
resultado.plot(
    ax=ax,
    column='count',
    scheme='quantiles',
    k=5,
    legend=True,
    figsize=(15, 15)
)

Durante este proceso, se discute la distribución de los datos, notando que existen valores atípicos extremos (municipios con 1 hurto frente a otros con cientos), lo cual puede dificultar la visualización y generar advertencias en la generación de las clases, sugiriendo ajustar el número de cuantiles para una mejor representación visual.

En esta clase utilizaremos la librería GeoPandas para realizar un análisis espacial enfocado en determinar qué colegios del departamento de Antioquia, Colombia, se verían afectados ante el supuesto desbordamiento de ríos. Para ello, aplicaremos conceptos de superposición espacial, consultas de atributos y generación de áreas de influencia o buffers.

Preparación del Entorno y Carga de Datos

Inicialmente, prepararemos el entorno de trabajo importando la librería necesaria bajo el alias común gpd. La carga de datos se realizará utilizando la función gpd.read_file(), la cual nos permite leer formatos de archivos espaciales como Shapefiles. En este ejercicio, cargaremos tres capas vectoriales distintas: los departamentos de Colombia, la red de drenajes (ríos) y los establecimientos educativos.

import geopandas as gpd

# Carga de capas vectoriales
departamentos = gpd.read_file(r'ruta/departamentos.shp')
rios = gpd.read_file(r'ruta/drenajes.shp')
colegios = gpd.read_file(r'ruta/establecimientos_educativos.shp')

# Verificación de datos
departamentos.head()
departamentos.plot()

Es fundamental recordar que al definir las rutas de los archivos en Windows, se antepone una r antes de las comillas para evitar errores de codificación en la lectura del path.

Filtrado del Área de Estudio

Procederemos a filtrar la información geográfica para acotar nuestra área de estudio. Dado que las capas cubren todo el territorio nacional, utilizaremos el método .query() para extraer únicamente el polígono correspondiente al departamento de Antioquia. Este paso es vital para optimizar el procesamiento, evitando realizar cálculos sobre zonas que no son de interés.

# Filtrar el departamento de Antioquia
antioquia = departamentos.query("DeNombre == 'Antioquia'")

Operación de Superposición (Overlay)

Con el área de estudio delimitada, realizaremos una operación espacial conocida como Overlay o superposición. El objetivo es recortar la capa de ríos para obtener solo aquellos que cruzan el departamento de Antioquia. Teóricamente, esta operación identifica la intersección geométrica entre las capas y devuelve únicamente las geometrías que comparten espacio en ambas, transfiriendo también sus atributos.

# Recortar ríos por Antioquia
rios_antioquia = gpd.overlay(rios, antioquia, how='intersection')

Verificación del Sistema de Referencia de Coordenadas (CRS)

Antes de proceder con el análisis de inundación, es obligatorio verificar el Sistema de Referencia de Coordenadas (CRS) de nuestras capas mediante el atributo .crs. Para realizar cálculos de distancias métricas, como un buffer, los datos deben estar proyectados en un sistema de coordenadas planas (en metros), como lo es Magna Sirgas para Colombia.

# Verificar el CRS
print(rios_antioquia.crs)

# Reproyectar si es necesario (de grados a metros)
rios_antioquia = rios_antioquia.to_crs(epsg=3116)

Si los datos estuvieran en coordenadas geográficas (grados latitud/longitud), sería necesario reproyectarlos con el método .to_crs() antes de continuar, ya que un buffer calculado en grados generaría resultados erróneos.

Generación de Zonas de Inundación (Buffer)

El núcleo de nuestro análisis consiste en simular el desbordamiento de los ríos. Para esto, aplicaremos el método .buffer() sobre la capa de ríos recortada, asignando un valor de 500 metros. Esto generará una nueva geometría de tipo polígono alrededor de cada línea de río.

Aquí introducimos un concepto avanzado de GeoPandas: un GeoDataFrame puede tener múltiples columnas geométricas, pero solo una es la geometría activa sobre la cual se realizan las operaciones espaciales. Por defecto, la geometría activa son las líneas de los ríos, por lo que usaremos el comando .set_geometry('buffer') para indicar que ahora queremos trabajar con los polígonos de inundación.

# Crear buffer de 500 metros alrededor de los ríos
rios_antioquia['buffer'] = rios_antioquia.buffer(500)

# Cambiar la geometría activa al buffer
rios_antioquia = rios_antioquia.set_geometry('buffer')

Identificación de Colegios en Riesgo

Finalmente, determinaremos qué colegios caen dentro de estas zonas de inundación. Realizaremos un nuevo overlay de tipo intersección, esta vez entre la capa de colegios filtrada para Antioquia y la capa de ríos (que ahora tiene activa la geometría del buffer).

# Filtrar colegios de Antioquia
colegios_antioquia = gpd.overlay(colegios, antioquia, how='intersection')

# Identificar colegios en zona de inundación
colegios_en_riesgo = gpd.overlay(colegios_antioquia, rios_antioquia, how='intersection')

# Cantidad de colegios afectados
total_colegios = colegios_antioquia.shape[0]
colegios_afectados = colegios_en_riesgo.shape[0]
porcentaje = (colegios_afectados / total_colegios) * 100

print(f'Colegios en riesgo: {colegios_afectados} de {total_colegios} ({porcentaje:.2f}%)')

El resultado será un nuevo conjunto de datos que contiene solo los colegios en riesgo. Para cuantificar el impacto, utilizaremos el atributo .shape, que nos devolverá la cantidad de filas (colegios afectados) y columnas, permitiéndonos presentar una estadística concreta del análisis.

En esta clase se abordará la integración de datos raster y vectoriales utilizando Python. Si bien GeoPandas es el estándar para manejar datos vectoriales (puntos, líneas y polígonos), no tiene capacidad nativa para manejar imágenes o modelos digitales de elevación. Para suplir esto, utilizaremos Rasterio, una librería fundamental en el ecosistema geoespacial de código abierto que permite leer y escribir formatos raster. Rasterio maneja la información espacial como arreglos N-dimensionales basados en NumPy, lo que facilita el cálculo matemático eficiente sobre las matrices de píxeles. El objetivo práctico será extraer valores de elevación de un Modelo Digital de Elevación (DEM) en ubicaciones específicas definidas por una capa de puntos shapefile (ciudades y nevados), para finalmente enriquecer la tabla de atributos de dichos puntos con la nueva información altimétrica.

Configuración del Entorno de Trabajo

Para replicar el ejercicio, se trabajará sobre un entorno de Anaconda previamente configurado (llamado por ejemplo "geo_env"). Es indispensable tener instaladas las librerías geopandas, rasterio y matplotlib. Una vez en Jupyter Notebook, se inicia el script importando los módulos necesarios.

import geopandas as gpd
import rasterio
import matplotlib.pyplot as plt
from rasterio.plot import show

Lectura y Visualización de Datos Raster

El primer paso lógico consiste en cargar la información raster. Se utiliza el método rasterio.open('ruta_del_archivo') asignándolo a una variable que contendrá el objeto raster. Este método es homólogo al read_file de GeoPandas pero optimizado para matrices. Rasterio soporta formatos como GeoTIFF (.tif) o HGT.

# Cargar el modelo de elevación
dem = rasterio.open(r'ruta/modelo_elevacion.tif')

# Visualización básica
show(dem)

# Visualización con curvas de nivel
show(dem, contour=True)

Con el comando show(objeto_raster), se despliega una visualización básica del modelo de elevación. Es posible refinar esta visualización agregando parámetros como contour=True para generar curvas de nivel automáticamente.

Carga de Datos Vectoriales y Visualización Conjunta

Posteriormente se cargan los puntos de muestreo utilizando GeoPandas. Al tener ambas fuentes de datos, es común querer visualizarlas en un mismo gráfico para validar la superposición espacial. Dado que las funciones de ploteo básicas pueden ser limitadas para superposiciones complejas, utilizaremos Matplotlib.

# Cargar puntos de muestreo
puntos = gpd.read_file(r'ruta/ciudades_nevados.shp')

# Crear figura y ejes
fig, ax = plt.subplots(figsize=(12, 12))

# Graficar raster y superponer puntos
show(dem, ax=ax)
puntos.plot(ax=ax, color='red', markersize=50)

Lógica de Extracción de Coordenadas

Para extraer los valores de los píxeles subyacentes a cada punto, Rasterio requiere una lista de coordenadas (pares X, Y). Dado que GeoPandas almacena la geometría en la columna geometry, se debe iterar sobre ella para separar las componentes. Se utiliza la función zip(), la cual combina las listas de coordenadas X y Y en pares ordenados (tuplas).

# Extraer coordenadas X, Y de cada punto
lista_coordenadas = [
    (x, y) for x, y in zip(puntos.geometry.x, puntos.geometry.y)
]

print(lista_coordenadas)

Muestreo y Asignación de Valores al DataFrame

El proceso central de extracción se realiza mediante el método .sample() del objeto raster. La instrucción raster.sample(lista_coordenadas) recorre el archivo raster y devuelve un generador con los valores de los píxeles en las ubicaciones indicadas.

# Muestrear valores de elevación
elevaciones = [valor[0] for valor in dem.sample(lista_coordenadas)]

# Agregar columna al GeoDataFrame
puntos['elevacion'] = elevaciones

# Visualizar resultados
puntos[['nombre', 'elevacion']]

Al visualizar nuevamente la tabla de atributos, cada punto (ciudad o nevado) tendrá asociada su altura exacta en metros sobre el nivel del mar derivada del DEM. Este flujo de trabajo es replicable para extraer cualquier tipo de información raster, como tipos de cobertura de suelo, índices de vegetación o variables climáticas, hacia una capa vectorial de puntos.

En esta clase se utilizará la librería Rasterio para el manejo de información en formato raster, enfocándonos principalmente en imágenes satelitales. Es fundamental comprender que un archivo raster (como un GeoTIFF) funciona como una matriz de valores numéricos, por lo que Rasterio se apoya en la estructura de NumPy (Numerical Python) para procesar estos datos como objetos n-dimensionales.

import rasterio
import numpy as np
import matplotlib.pyplot as plt
from rasterio.plot import show, show_hist

Carga y Exploración de Metadatos

El primer paso práctico consiste en la lectura de los datos. Se empleará el método rasterio.open('ruta_archivo.tif') para cargar la imagen en una variable. Una vez cargado el objeto, es crucial explorar sus metadatos antes de cualquier análisis visual.

# Cargar imagen satelital
imagen = rasterio.open(r'ruta/imagen_landsat.tif')

# Explorar metadatos
print(f'Bandas: {imagen.count}')
print(f'Resolución: {imagen.width} x {imagen.height} píxeles')
print(f'Extensión: {imagen.bounds}')
print(f'CRS: {imagen.crs}')
print(f'Índices disponibles: {imagen.indexes}')

Visualización Básica y Simbología

Para visualizar la información geoespacial, se hará uso de la función show(objeto_raster). Dado que las imágenes satelitales crudas pueden no tener el contraste adecuado, se explicará cómo visualizar bandas individuales aplicando mapas de colores (colormaps).

En términos teóricos, se abordará el caso de uso de una imagen Landsat 5, donde las bandas 1, 2 y 3 corresponden al espectro visible (azul, verde, rojo respectivamente) y la banda 4 al infrarrojo cercano.

# Visualizar banda con mapa de colores
show(imagen.read(4), cmap='Spectral')

# Visualizar bandas RGB con colormaps intuitivos
show(imagen.read(3), cmap='Reds')    # Banda roja
show(imagen.read(2), cmap='Greens')  # Banda verde
show(imagen.read(1), cmap='Blues')   # Banda azul

Manipulación de Arreglos con NumPy y Estadísticas

Una parte central de la guía es la extracción de valores. Se utilizará el método .read(indice_banda) para extraer una banda específica y almacenarla en una nueva variable. Al hacer esto, la variable se convierte automáticamente en un array de NumPy, lo que habilita el uso de funciones estadísticas directas.

# Extraer banda como array NumPy
banda_nir = imagen.read(4)

# Estadísticas básicas
print(f'Valor mínimo: {banda_nir.min()}')
print(f'Valor máximo: {banda_nir.max()}')
print(f'Promedio: {banda_nir.mean():.2f}')

Análisis de Histogramas

Para interpretar la distribución de los niveles digitales de la imagen, se generarán histogramas utilizando la función show_hist. Se configurarán parámetros como el título, el número de barras y la transparencia para mejorar la legibilidad.

# Generar histograma de la imagen
show_hist(
    imagen,
    title='Distribución de valores por banda',
    bins=50,
    alpha=0.7
)

Composición de Bandas (Falso Color)

Finalmente, se explicará cómo realizar composiciones de bandas (Composite) para análisis multiespectral avanzado. Dado que Rasterio maneja los datos por separado, se utilizará la función np.dstack() de NumPy para apilar tres bandas individuales en un solo arreglo tridimensional. Se generarán combinaciones teóricas clásicas, como el Falso Color Estándar (Infrarrojo Cercano, Rojo, Verde) para resaltar vegetación.

# Extraer bandas para composición
banda_r = imagen.read(3)  # Rojo
banda_g = imagen.read(2)  # Verde
banda_nir = imagen.read(4)  # Infrarrojo cercano

# Falso color: NIR, Rojo, Verde
falso_color = np.dstack((banda_nir, banda_r, banda_g))

# Visualizar composición
fig, ax = plt.subplots(figsize=(12, 12))
plt.imshow(falso_color)
plt.title('Composición Falso Color (NIR-R-G)')
plt.show()

Debido a que las imágenes pueden visualizarse pequeñas por defecto, se ajustará el tamaño del gráfico definiendo un entorno de figura con plt.subplots(figsize=(12, 12)), permitiendo una interpretación detallada de coberturas como nieve, nubes o vegetación.