12  Formato de datos espaciales - Raster

En este capítulo vamos a aprender acerca de los dos principales formatos de datos espaciales que son útiles en la Ciencia de los datos de Tierra, ellos son el raster y el vector.

12.1 Uso de datos raster para la Ciencia de los datos de la Tierra

¿Qué son los datos raster?

Los rasters o las cuadrículas, son datos almacenados como una grilla de valores, los cuales son dibujados en mapas como pixeles. Cada valor de pixel representa un área de la superficie de la tierra formando el dato espacial. Un archivo raster está compuesto de una cuadrícula regular de celdas, las cuales poseen el mismo tamaño. Es común que ya hayamos usados datos raster, ya sea en fotografías digitales o en una imagen de Google Earth. Sin embargo, los raster que vamos a trabajar son diferentes de las fotografías, puesto que se encuentran espacialmente referenciados. Cada pixel representa un área de la superficie de la tierra. Esta área está definida por la resolución espacial del raster.

Figura 12.1: Representación y elementos de un dato raster.

Los datos raster pueden tener una o más layers

Los datos raster pueden tener una o más layers (capas o variables). Por ejemplo, un modelo de elevación, generalmente solo ofrecerá una layer que representa la elevación de la superficie de la tierra en un lugar determinado. Sin embargo, otros datos que incluyen imágenes y series temporales de datos, pueden resultar en un raster que está compuesto de varias layers. Distintos tipos de archivos pueden ser usados para acomodar diferentes tamaños y estructuras de datos raster.

Diferentes formatos de archivos raster

Existen muchos tipos diferentes de archivos que son empleados para almacenar datos raster.

Datos raster almacenados como archivos únicos

Muchos de los sets de datos que existen para la libre distribución como las imágenes Landsat están almacenadas como formato de archivos únicos. Para las imágenes Landsat, a menudo se encontrará cada banda almacenada como un archivo .tif por separado. Sin embargo, es posible también encontrar otros casos en donde todas las bandas se almacenen en sólo un archivo .tif. Los tipos de archivos más comunes que son usados para almacenar archivos únicos incluyen:

• .tif/.tiff: soporte para Tagged Image File Format. Es uno de las formas más comunes de almacenar datos raster. Algunas imágenes de satélite como Landsat, comparten sus datos de esta manera.

• .asc: soporte para* ASCII Raster Files*. Es un formato basado en texto que almacena datos raster. Este formato es usado dado que es simple de almacenar y distribuir.

Formato jerárquico de datos

Los formatos jerárquicos de datos pueden almacenar muchos tipos de datos en un solo archivo. Estos formatos son adecuados para largos sets de datos, así como si se desea trabajar con una muestra o con solo una parte de los datos a la vez. Los datos jerárquicos pueden involucrar un poco más de trabajo, pero tienden a hacer del proceso algo más eficiente. Los tipos de archivos más comunes que almacenan estos datos son:

• .hdf/.hdf5: soporte para formato de datos jerárquicos. Uno de los más comunes para almacenar datos raster. Algunas imágenes de satélite como MODIS, lo usan para compartir los datos.

• .nc (NetCDF): soporte para Network Common Data Form. Una forma común para almacenar datos de clima.

¿Qué tipo de datos son almacenados en raster?

Algunos de los datos que a menudo son suministrados en formato raster pueden ser:

• Imágenes de satélite

• Uso del suelo sobre áreas extensas

• Datos de elevación

• Datos climáticos

• Datos de batimetría

El siguiente código permite descargar datos raster que serán abiertos posteriormente.

# Importa los paquetes necesarios.
# Proporciona funciones para interactuar con el sistema operativo.
import os

# Paquetes para leer datos raster y vectoriales.
import matplotlib.pyplot as plt
import numpy as np
os.environ['USE_PYGEOS'] = '0'
import geopandas as gpd
import rasterio as rio

# Earthpy es un paquete de earthlab para trabajar con datos espaciales.
import earthpy as et
import earthpy.plot as ep

/usr/share/miniconda/envs/tallerpython/lib/python3.10/site-packages/earthpy/__init__.py:7: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
  from pkg_resources import resource_string

# Obtiene los datos y setea un directorio de trabajo.
et.data.get_data("colorado-flood")
os.chdir(os.path.join(et.io.HOME, 'earth-analytics', 'data'))

# Para abrir el archivo primero crear la ruta a él.
lidar_dtm_path = os.path.join("colorado-flood",
                              "spatial",
                              "boulder-leehill-rd",
                              "pre-flood",
                              "lidar",
                              "pre_DTM.tif")
print(lidar_dtm_path)

'colorado-flood\spatial\boulder-leehill-rd\pre-flood\lidar\pre_DTM.tif'

Luego crear la conexión al archivo. La expresión with rio.open() crea lo que se llama un administrador de contexto para abrir archivos. Esto permite crear una conexión con el archivo, sin modificar el archivo propiamente tal.

# Luego crear una conexión con el archivo.
with rio.open(lidar_dtm_path) as src:
    # Leer el dato de entrada y asigna a una variable (lidar_dtm).
    lidar_dtm = src.read(1)

La siguiente expresión permite observar el contenido del archivo recientemente abierto.

lidar_dtm

array([[-3.4028235e+38, -3.4028235e+38, -3.4028235e+38, ...,
         1.6956300e+03,  1.6954199e+03,  1.6954299e+03],
       [-3.4028235e+38, -3.4028235e+38, -3.4028235e+38, ...,
         1.6956000e+03,  1.6955399e+03,  1.6953600e+03],
       [-3.4028235e+38, -3.4028235e+38, -3.4028235e+38, ...,
         1.6953800e+03,  1.6954399e+03,  1.6953700e+03],
       ...,
       [-3.4028235e+38, -3.4028235e+38, -3.4028235e+38, ...,
         1.6814500e+03,  1.6813900e+03,  1.6812500e+03],
       [-3.4028235e+38, -3.4028235e+38, -3.4028235e+38, ...,
         1.6817200e+03,  1.6815699e+03,  1.6815599e+03],
       [-3.4028235e+38, -3.4028235e+38, -3.4028235e+38, ...,
         1.6818900e+03,  1.6818099e+03,  1.6817400e+03]],
      shape=(2000, 4000), dtype=float32)

Si queremos enmascarar los valores que se emplean para indicar los datos no válidos tenemos que hacer una pequeña modificación.

# Luego crear una conexión con el archivo.
with rio.open(lidar_dtm_path) as src:
    # Leer el dato de entrada y asigna a una variable (lidar_dtm).
    lidar_dtm = src.read(1, masked=True)

lidar_dtm

masked_array(
  data=[[--, --, --, ..., 1695.6300048828125, 1695.419921875,
         1695.429931640625],
        [--, --, --, ..., 1695.5999755859375, 1695.5399169921875,
         1695.3599853515625],
        [--, --, --, ..., 1695.3800048828125, 1695.43994140625,
         1695.3699951171875],
        ...,
        [--, --, --, ..., 1681.449951171875, 1681.3900146484375, 1681.25],
        [--, --, --, ..., 1681.719970703125, 1681.5699462890625,
         1681.5599365234375],
        [--, --, --, ..., 1681.8900146484375, 1681.8099365234375,
         1681.739990234375]],
  mask=[[ True,  True,  True, ..., False, False, False],
        [ True,  True,  True, ..., False, False, False],
        [ True,  True,  True, ..., False, False, False],
        ...,
        [ True,  True,  True, ..., False, False, False],
        [ True,  True,  True, ..., False, False, False],
        [ True,  True,  True, ..., False, False, False]],
  fill_value=-3.4028235e+38,
  dtype=float32)

En términos simples el código hace:

1. rio.open(): rio es el alias para rasterio (ojo a la sección de importaciones).
2. open(): crea una conexión con el archivo en nuestro computador.
3. src.read(): trae el dato dentro de Python, de manera que lo podamos usar en el código.
4. masked=True: enmascara todos los valores nodata en el arreglo. Esto implica que esos valores no serán ploteados y tampoco incluidos en los cálculos matemáticos que se hagan con Python.

El dato que se está abriendo cuyo nombre de archivo es pre_DTM.tif es un modelo digital de elevaciones derivado de datos LiDAR. Es fácil de notar que el código empleado para abrir este tipo de datos es un poco más complicado que el usado para abrir los datos tabulares con Panda.

12.2 Explorar los valores del dato raster y su estructura

A continuación, vamos a dar una mirada a los datos. Note que la estructura del dato devuelta por type() es un objeto Python retornado por rasterio, el cual es un array del tipo numpy. Los arrays de numpy son maneras eficientes de almacenar y trabajar con datos raster en Python.

# Conocer el tipo de objeto.
print(type(lidar_dtm))

# Listar los valores minimos y maximos del array.
print(lidar_dtm.min(), lidar_dtm.max())

# Listar las dimensiones de array (filas, columnas).
print(lidar_dtm.shape)

<class 'numpy.ma.MaskedArray'>
1676.21 2087.43
(2000, 4000)

Si queremos observar cómo se ven nuestros datos en el territorio los podemos plotear. Para eso usamos una función del paquete earthpy con los parámetros de visualización.

# Ploteamos el modelo digital de elevación.

f, ax = plt.subplots()

ep.plot_bands(lidar_dtm,
              scale=False,
              cmap='Greys',
              title="Modelo Digital de Elevación LiDAR",
              ax=ax)
plt.show()

12.3 Las Imágenes, otro tipo de datos raster

Otro tipo de datos raster que se pueden observar son las imágenes. Estas serán reconocidas por ustedes si han usado Google Earth u otro tipo de herramienta que tenga un layer de imágenes. En Python es posible abrir y plotear datos de imágenes. En el código que se muestra a continuación vamos a descargar algunos datos de imágenes que fueron colectados antes de un incendio ocurrido cerca de Colorado, EE.UU.

# Descarga datos de usando una dirección Internet.
et.data.get_data(url="https://ndownloader.figshare.com/files/23070791");

Downloading from https://ndownloader.figshare.com/files/23070791
Extracted output to /home/runner/earth-analytics/data/earthpy-downloads/naip-before-after

Crea una ruta apuntando al archivo descargado.

# Crea una ruta apuntando al archivo con los datos. Note que es un archivo .tif.
naip_pre_fire_path = os.path.join("earthpy-downloads",
                                  "naip-before-after",
                                  "pre-fire",
                                  "crop",
                                  "m_3910505_nw_13_1_20150919_crop.tif")

print(naip_pre_fire_path)

'earthpy-downloads\naip-before-after\pre-fire\crop\m_3910505_nw_13_1_20150919_crop.tif'

Abre el archivo.

# Abre el archivo usando rasterio.

with rio.open(naip_pre_fire_path) as naip_prefire_src:
    naip_pre_fire = naip_prefire_src.read()

naip_pre_fire

array([[[113, 117, 137, ...,  54,  51,  74],
        [113, 117, 131, ...,  63,  54,  54],
        [111, 117, 120, ...,  78,  76,  52],
        ...,
        [191, 192, 193, ...,  58,  69,  76],
        [192, 192, 193, ...,  53,  62,  71],
        [193, 193, 193, ...,  51,  59,  66]],

       [[114, 114, 126, ...,  58,  54,  72],
        [114, 112, 120, ...,  70,  60,  58],
        [111, 114, 115, ...,  85,  87,  58],
        ...,
        [183, 184, 185, ...,  61,  75,  84],
        [184, 185, 185, ...,  56,  66,  78],
        [186, 186, 186, ...,  52,  58,  65]],

       [[ 80,  87,  95, ...,  55,  54,  63],
        [ 79,  83,  90, ...,  57,  55,  55],
        [ 81,  84,  87, ...,  62,  65,  55],
        ...,
        [161, 161, 163, ...,  54,  58,  64],
        [162, 164, 165, ...,  53,  58,  62],
        [165, 166, 166, ...,  51,  54,  57]],

       [[145, 143, 139, ...,  74,  47,  65],
        [145, 146, 139, ...,  98,  59,  57],
        [142, 144, 144, ..., 119, 107,  54],
        ...,
        [167, 169, 172, ...,  76, 105, 119],
        [170, 171, 171, ...,  60,  81, 105],
        [173, 173, 174, ...,  44,  55,  67]]],
      shape=(4, 2312, 4377), dtype=int16)

El plotear una imagen es un poco distinto, ya que una imagen suele estar conformada por múltiples bandas. No obstante, es posible plotearlas de forma individual usando plot_bands(), o se puede plotear una composición en color usando plot_rgb().

# Plotea cada banda de la imagen de forma separada.
ep.plot_bands(naip_pre_fire, figsize=(10, 5))

array([[<Axes: title={'center': 'Band 1'}>,
        <Axes: title={'center': 'Band 2'}>,
        <Axes: title={'center': 'Band 3'}>],
       [<Axes: title={'center': 'Band 4'}>, <Axes: >, <Axes: >]],
      dtype=object)

# Plotea una composición de color.
f, ax = plt.subplots()

ep.plot_rgb(naip_pre_fire, title="NAIP data pre-incendio", ax=ax)

plt.show()

12.4 Visualización Raster

La visualización de datos ráster es fundamental en el análisis de información geográfica porque permite interpretar y comunicar patrones espaciales complejos de manera intuitiva. Los datos ráster, como imágenes satelitales, modelos de elevación digital o mapas de uso del suelo, representan fenómenos continuos en el espacio, y su visualización facilita:

1. Exploración y detección de patrones: Ayuda a identificar tendencias, anomalías o relaciones espaciales que serían difíciles de discernir a partir de tablas o valores numéricos.

2. Toma de decisiones informada: Permite que los responsables de políticas y usuarios finales comprendan mejor los datos geográficos para tomar decisiones precisas, como planificación territorial o manejo de recursos naturales.

3. Análisis multiescala: Proporciona una forma visual de evaluar fenómenos tanto a nivel local como global, integrando múltiples capas de información.

4. ** Comunicación efectiva**: Transforma datos técnicos en representaciones gráficas accesibles, facilitando el entendimiento entre audiencias diversas.

Para ejemplificar esto, cargaremos la imagen img-rgb.tif disponible en el siguiente link. Esta es una imagen compuesta por tres bandas, rojo, verde y azul. Visualicémosla:

<Figure size 672x480 with 0 Axes>

raster = rio.open("img-rgb.tif")

raster_data = raster.read()

plt.figure()
ep.plot_rgb(raster_data)
plt.show()

Como se habrá notado, la imagen se visualiza con un tono oscuro que dificulta la identificación de sus componentes. El primer paso consiste en analizar la distribución de los valores en cada banda, en este caso, el orden de las bandas es rojo, verde y azul.

plt.figure()
ep.hist(raster_data,
        colors = ['r', 'g', 'b'],
        title = ['Red', 'Green', 'Blue'],
        cols = 3,
        figsize = (8, 4))
        
plt.show()

<Figure size 672x480 with 0 Axes>

En el histograma se puede observar que los datos se encuentran mayormente hacia el lado izquierdo y que hay pocos valores en el lado derecho, de esta forma al visualizar la imagen en una composición, esta resultará con poco contraste.

Dentro de la función plot_rgb() podemos definir dos parámetros que ayudarán para mejorar la visualización:

1. Estiramiento: El estiramiento (stretch = True) modifica el rango de valores de píxeles de una imagen para que se extienda a lo largo de todo el rango de visualización que entrega el monitor. En este sentido, se cambian los valores digitales para usar el rango disponible, en definitiva, se cambia el histograma. Esto incrementa el contraste de los objetos y mejora la diferenciación de sus entornos, las áreas con tonos claros aparecen más brillantes y las con tonos oscuros aparecen más oscuras, haciendo que la interpretación visual sea más fácil.

La fórmula que se aplica sobre los ND de la imagen es:

\[ DN_{st}=255*\frac{(DN-DN_{min})}{(DN_{max}-DN_{min})} \tag{12.1}\]

2. Corta de colas de histogramas: str_clip=2 permite especificar qué parte de las colas de los datos desea recortar, esto es para evitar los números extremos. Cuanto mayor sea el número, más se estirarán o aclararán los datos.

f, ax = plt.subplots(figsize=(8,6))

ep.plot_rgb(raster_data,
            rgb=[0,1,2],
            stretch = True,
            str_clip = 1,
            ax = ax)
plt.show()

Figura 12.2: Representación del cambio de distribución de valores en el histograma.

Otro método para mejorar la visualización de imágenes es manipular directamente los datos (aunque siempre hay que tener cuidado de no alterar los datos originales), para ello nos podemos ayudar de la librería Skimage, diseñada para el procesamiento de imágenes.

from skimage import exposure

Podemos ajustar la corrección gamma que se utiliza para compensar ciertas propiedades de la percepción humana de la luz y el color, con el objetivo de maximizar el uso eficiente del ancho de banda o la distribución de bits en imágenes digitales. Sin la corrección gamma, se asignan demasiados bits o ancho de banda a los valores de intensidad más brillantes, que el ojo humano tiene dificultad para diferenciar, y muy pocos bits a las intensidades más oscuras, a las que somos más sensibles y que requieren mayor resolución para mantener la misma calidad visual.

\[ Valor_{out}=gain*(Valor_{in}^{gamma}) \tag{12.2}\]

Además del ajuste gamma, el parámetro gain (ganancia) se utiliza para escalar la intensidad general de la imagen después de aplicar la corrección gamma. Un valor de gain mayor que 1 aumenta el brillo global de la imagen, mientras que un valor menor que 1 lo reduce. La combinación de gamma y gain permite un control más preciso, ajustando tanto el contraste no lineal como el brillo general.

En términos prácticos:

• Para una gamma superior a 1, el histograma de la imagen se desplazará hacia la izquierda, oscureciendo la imagen de salida en comparación con la entrada.

• Para valores de gamma menores que 1, el histograma se desplazará hacia la derecha, haciendo que la imagen de salida sea más brillante que la de entrada.

• El gain puede compensar estos efectos, ajustando el brillo general sin alterar la curva de contraste determinada por el gamma.

gamma_aj = exposure.adjust_gamma(raster_data, 
                                 gamma = 1.5, 
                                 gain = 6)

plt.figure()

ep.hist(gamma_aj,
        colors = ['r', 'g', 'b'],
        title = ['Red', 'Green', 'Blue'],
        cols = 3,
        figsize = (8, 4),
        bins = 20)

plt.show()

<Figure size 672x480 with 0 Axes>

f, ax = plt.subplots(figsize = (8,6))

ep.plot_rgb(gamma_aj,
            rgb = [0,1,2],
            stretch = False,
            ax = ax)
plt.show()

Ecualización de histograma

Este método aumenta el contraste global de las imágenes, utilizando el rango completo de valores posibles, resultando en un histograma más uniforme respecto al original. Es especialmente útil cuando el histograma original contiene la mayor parte de los datos en un rango pequeño, provocando problemas en el contraste.

eq = exposure.equalize_hist(raster_data)

ep.hist(eq*255, # Podemos volver los datos al rango 0-255.
        colors=['r', 'g', 'b'],
        title=['Red', 'Green', 'Blue'],
        cols=3,
        figsize=(8, 4),
        bins=12);

plt.show()

Así, obtenemos la siguiente composición:

f, ax = plt.subplots(figsize=(8,6))

ep.plot_rgb(eq,
            rgb=[0,1,2],
            ax=ax)
            
plt.show()

A pesar de la mejora respecto a la original, la imagen se ve un poco sobreexpuesta, si deseamos recortar los valores extremos del histograma, podemos usar la variante de esta función a la cual se le puede introducir un parámetro para el corte de las colas. Pruebe con varios valores para observar cómo se desplazan los valores en histograma.

eq_ad = exposure.equalize_adapthist(raster_data,
                                    clip_limit = 0.03)
                                    
ep.hist(eq_ad,
        colors=['r', 'g', 'b'],
        title=['Red', 'Green', 'Blue'],
        cols=3,
        figsize=(8, 4),
        bins=12);

plt.show()

f, ax = plt.subplots(figsize=(8,6))

ep.plot_rgb(eq_ad,
            rgb=[0,1,2],
            ax=ax)
            
plt.show()

Estas funciones generan una versión procesada de los datos originales como un array NumPy, sin modificar el archivo original. Si se desea guardar el resultado como un archivo raster georreferenciado, se puede exportar utilizando rasterio.