Vamos a aplicar algunas de las herramientas que hemos visto hasta ahora para obtener algunas visualizaciones más sofisticadas. Estas incluirán el uso de datos vectoriales de un GeoDataFrame de GeoPandas, se construirán gráficos estáticos y dinámicos a partir de un arreglo 3D y se combinarán datos vectoriales y datos ráster.
Como contexto, los archivos que examinaremos se basan en el incendio McKinney que ocurrió en el 2022, en el Bosque Nacional Klamath (al oeste del condado de Siskiyou, California). Los datos vectoriales representan una instantánea del límite de un incendio forestal. Los datos ráster corresponden a la alteración que se observó en la superficie de la vegetación (esto se explicará con mayor detalle más adelante).
Importación preliminar y direcciones de los archivos¶
Para empezar, se necesitan algunas importaciones típicas de paquetes. También definiremos algunas direcciones a archivos locales que contienen datos geoespaciales relevantes.
from warnings import filterwarnings
filterwarnings('ignore')
from pathlib import Path
import numpy as np, pandas as pd, xarray as xr
import geopandas as gpd
import rioxarray as rio# Imports for plotting
import hvplot.pandas, hvplot.xarray
import geoviews as gv
from geoviews import opts
gv.extension('bokeh')FILE_STEM = Path.cwd().parent if 'book' == Path.cwd().parent.stem else 'book'
ASSET_PATH = Path(FILE_STEM, 'assets')
SHAPE_FILE = ASSET_PATH / 'shapefiles' / 'mckinney' / 'McKinney_NIFC.shp'
RASTER_FILES = list(ASSET_PATH.glob('OPERA*VEG*.tif'))
RASTER_FILE = RASTER_FILES[0]Graficar datos vectoriales desde un GeoDataFrame¶
El paquete GeoPandas contiene muchas herramientas útiles para trabajar con datos geoespaciales vectoriales. En particular, el GeoDataFrame de GeoPandas es una subclase del DataFrame de Pandas que está específicamente diseñado, por ejemplo, para datos vectoriales almacenados en shapefiles.
Como ejemplo, vamos a cargar algunos datos vectoriales desde la ruta local SHAPEFILE (tal como se definió anteriormente).
shape_df = gpd.read_file(SHAPE_FILE)
shape_dfEl objeto shape_df es un GeoDataFrame que tiene más de dos docenas de columnas de metadatos en una sola fila. La columna principal que nos interesa es la de geometry. Esta columna almacena las coordenadas de los vértices de un MULTIPOLYGON, es decir, un conjunto de polígonos.
shape_df.geometryLos vértices de los polígonos parecen expresarse como pares de (longitude, latitude). Podemos determinar qué sistema de coordenadas específico se utiliza examinando el atributo crs del GeoDataFrame.
print(shape_df.crs)Utilicemos hvplot para generar un gráfico de este conjunto de datos vectoriales.
shape_df.hvplot()La proyección en este gráfico es un poco extraña. La documentación de HoloViz incluye una discusión sobre las consideraciones cuando se visualizan datos geográficos. El punto importante a recordar en este contexto es que la opción geo=True es útil.
Vamos a crear dos diccionarios de Python, shapeplot_opts y layout_opts, y construiremos una visualización.
shapeplot_opts= dict(geo=True)
layout_opts = dict(xlabel='Longitude', ylabel="Latitude")
print(f"{shapeplot_opts=}")
print(f"{layout_opts=}")
shape_df.hvplot(**shapeplot_opts).opts(**layout_opts)La documentación de HoloViz incluye una discusión sobre las consideraciones cuando se grafican datos geográficos. El punto importante que se debe recordar en el contexto inmediato es que la opción geo=True es útil.
- Definir
color=Nonesignifica que los polígonos no se rellenarán. - Al especificar
line_widthyline_colorse modifica la apariencia de los límites de los polígonos.
shapeplot_opts.update( color=None ,
line_width=2,
line_color='red')
print(shapeplot_opts)
shape_df.hvplot(**shapeplot_opts).opts(**layout_opts)Rellenemos los polígonos con color naranja utilizando la opción color=orange y hagamos que el área rellenada sea parcialmente transparente especificando un valor para alpha entre cero y uno.
shapeplot_opts.update(color='orange', alpha=0.25)
print(shapeplot_opts)
shape_df.hvplot(**shapeplot_opts).opts(**layout_opts)Agregado de un mapa base¶
A continuación, proporcionemos un mapa base y superpongamos los polígonos trazados utilizando el operador *. Observa el uso de paréntesis para aplicar el método opts a la salida del operador *.
basemap = gv.tile_sources.ESRI(height=500, width=500, padding=0.1)
(shape_df.hvplot(**shapeplot_opts) * basemap).opts(**layout_opts)El mapa base no necesita ser superpuesto como un objeto separado, puede especificarse utilizando la palabra clave tiles. Por ejemplo, al establecer tiles='OSM' se utilizan los mosaicos de OpenStreetMap. Observa que las dimensiones de la imagen renderizada probablemente no sean las mismas que las de la imagen anterior con los mosaicos ESRI, ya que antes especificamos explícitamente la altura de 500, con la opción height=500, y el ancho de 500, con la opción width=500, en la definición del basemap.
shapeplot_opts.update(tiles='OSM')
shape_df.hvplot(**shapeplot_opts).opts(**layout_opts)Vamos a eliminar la opción tiles de shapeplot_opts y a vincular el objeto del gráfico resultante al identificador shapeplot.
del shapeplot_opts['tiles']
print(shapeplot_opts)
shapeplot = shape_df.hvplot(**shapeplot_opts)
shapeplotCombinación de datos vectoriales con datos ráster en un gráfico estático¶
Combinemos estos datos vectoriales con algunos datos ráster. Cargaremos los datos raster desde un archivo local utilizando la función rioxarray.open_rasterio. Por practicidad, usaremos el encadenamiento de métodos para reetiquetar las coordenadas del DataArray cargado y usaremos el método squeeze para convertir el arreglo tridimensional que se cargó en uno bidimensional.
raster = (
rio.open_rasterio(RASTER_FILE)
.rename(dict(x='longitude', y='latitude'))
.squeeze()
)
rasterUtilizaremos un diccionario de Python image_opts para almacenar argumentos útiles y pasarlos a hvplot.image.
image_opts = dict(
x='longitude',
y='latitude',
rasterize=True,
dynamic=True,
cmap='hot_r',
clim=(0, 100),
alpha=0.8,
project=True,
)
raster.hvplot.image(**image_opts)Podemos superponer shapeplot con los datos ráster graficados utilizando el operador *. Podemos utilizar las herramientas Pan, Wheel Zoom y Box Zoom en la barra de herramientas de Bokeh a la derecha del gráfico para hacer zoom y comprobar que los datos vectoriales fueron trazados encima de los datos ráster.
raster.hvplot.image(**image_opts) * shapeplotAdemás, podemos superponer los datos vectoriales y ráster en mosaicos ESRI utilizando basemap.
raster.hvplot.image(**image_opts) * shapeplot * basemapObserva que muchos de los píxeles blancos oscurecen el gráfico. Resulta que estos píxeles corresponden a ceros en los datos ráster, los cuales pueden ser ignorados sin problema. Podemos filtrarlos utilizando el método where.
raster = raster.where((raster!=0))
layout_opts.update(title="McKinney 2022 Wildfires")
(raster.hvplot.image(**image_opts) * shapeplot * basemap).opts(**layout_opts)Construcción de gráficos estáticos a partir de un arreglo 3D¶
Vamos a cargar una secuencia de archivos ráster en un arreglo tridimensional y generaremos algunos gráficos. Para empezar, construiremos un bucle para leer DataArrays de los archivos RASTER_FILES y utilizaremos xarray.concat para generar un único arreglo de rásters tridimensional (es decir, tres rásters de apilados verticalmente). Aprenderemos las interpretaciones específicas asociadas con el conjunto de datos ráster en un cuaderno computacional posterior. Por el momento, los trataremos como datos sin procesar con los que experimentaremos.
stack = []
for path in RASTER_FILES:
data = rio.open_rasterio(path).rename(dict(x='longitude', y='latitude'))
band_name = path.stem.split('_')[-1]
data.coords.update({'band': [band_name]})
data.attrs = dict(description=f"OPERA DIST product", units=None)
stack.append(data)
stack = xr.concat(stack, dim='band')
stack = stack.where(stack!=0)Renombramos los ejes longitude y latitude y filtramos los píxeles con valor cero para simplificar la visualización posterior.
stackUna vez que el stack de DataArray esté construido, podemos enfocarnos en la visualización.
Si queremos generar un gráfico estático con varias imágenes, podemos utilizar hvplot.image junto con el método .layout. Para ver cómo funciona, empecemos por redefinir los diccionarios image_opts y layout_opts.
image_opts = dict( x='longitude',
y='latitude',
rasterize=True,
dynamic=True,
cmap='hot_r',
crs=stack.rio.crs,
alpha=0.8,
project=True,
aspect='equal',
shared_axes=False,
colorbar=True,
tiles='ESRI',
padding=0.1)
layout_opts = dict(xlabel='Longitude', ylabel="Latitude")Para acelerar el renderizado, construiremos inicialmente un objeto view que seleccione subconjuntos de píxeles. Definimos inicialmente el parámetro steps=200 para restringir la vista a cada 200 píxeles en cualquier dirección. Si reducimos los steps, se tarda más en renderizar. Establecer steps=1 o steps=None equivale a seleccionar todos los píxeles.
steps = 200
subset = slice(0, None, steps)
layout_opts.update(frame_width=250, frame_height=250)
view = stack.isel(latitude=subset, longitude=subset)
view.hvplot.image(**image_opts).opts(**layout_opts).layout()El método layout traza de manera predeterminada cada uno de los tres rásteres seleccionados a lo largo del eje band horizontalmente.
Construcción de una vista dinámica a partir de un arreglo 3D¶
Otra forma de visualizar un arreglo tridimensional es asociar un widget de selección a uno de los ejes. Si llamamos a hvplot.image sin agregar el método .layout el resultado es un mapa dinámico. En este caso, el widget de selección nos permite elegir cortes a lo largo del eje band.
Una vez más, aumentar el valor del parámetro steps reduce el tiempo de renderizado. Reducirlo a 1 o a None renderiza a resolución completa.
steps = 200
subset = slice(0, None, steps)
layout_opts.update(frame_height=400, frame_width=400)
view = stack.isel(latitude=subset, longitude=subset)
view.hvplot.image(**image_opts).opts(**layout_opts)Más adelante apilaremos muchos rásters con distintas marcas temporales a lo largo de un eje time. Cuando haya muchos cortes, el widget de selección se renderizará como un deslizador en vez de como un menú desplegable. Podemos controlar la ubicación del widget utilizando una opción de palabra clave widget_location.
view.hvplot.image(widget_location='bottom_left', **image_opts, **layout_opts)Observa que agregar la opción widget_location modifica ligeramente la secuencia en la que se especifican las opciones. Es decir, si invocamos algo como
view.hvplot.image(widget_location='top_left', **image_opts).opts(**layout_opts)se genera una excepción (de ahí la invocación en la celda de código anterior). Hay algunas dificultades sutiles en la elaboración de la secuencia correcta para aplicar opciones al personalizar objetos HoloViz/Hvplot (principalmente debido a las formas en que los espacios de los nombres se superponen con Bokeh u otros backends de renderizado). La mejor estrategia es empezar con el menor número de opciones posible y experimentar añadiendo opciones de forma incremental hasta que obtengamos la visualización final que deseamos.
Combinación de datos vectoriales con datos ráster en una vista dinámica.¶
Por último, vamos a superponer nuestros datos vectoriales anteriores, los límites del incendio forestal, con la vista dinámica de este arreglo tridimensional de rásteres. Podemos utilizar el operador * para combinar la salida de hvplot.image con shapeplot, la vista renderizada de los datos vectoriales.
steps = 200
subset = slice(0, None, steps)
view = stack.isel(latitude=subset, longitude=subset)
image_opts.update(colorbar=False)
layout_opts.update(frame_height=500, frame_width=500)
(view.hvplot.image(**image_opts) * shapeplot).opts(**layout_opts)Una vez más, la correcta especificación de las opciones puede requerir un poco de experimentación. Las ideas importantes que se deben tomar en cuenta aquí son:
- cómo cargar conjuntos de datos relevantes con
geopandasandrioxarray, y - cómo utilizar
hvplotde forma interactiva e incremental para generar visualizaciones atractivas.