Radares Meteorológicos
Resumen
En este cuadernillo (Notebook) aprenderemos:
Breve introducción a los radares meterológicos
Red de Nacional de radares meteorológicos de Colombia
Consulta de datos usando el bucket de AWS de radares
Generación de gráficas de reflectividad
Generación de campos de lluvia usando la ecuación de Marshall & Gunn (1953)
Prerequisitos
Conceptos |
Importancia |
Notas |
---|---|---|
Necesario |
Lectura de datos multidimensionales |
|
Necesario |
Lectura de datos de radar |
|
Necesario |
Lectura de datos de radar |
|
Necesario |
Lectura de datos de radar |
|
Útil |
Fundametos básicos en radares meteorológicos |
|
Útil |
Entender la metadata de los datos |
Tiempo de aprendizaje: 30 minutos
1. Radares meteorológicos
Los radares meteorológicos son sensores activos que emiten pulsos de energía a la atmósfera y miden la energía retrodispersa de todos los objetos a lo largo del haz de medición como se puede observar en la siguiente animación del curso Fundamentos en radares meterológicos (The Commet program - meted.ucar.edu)
Tradicionalmente, los radares meteorológicos operan en modo de Indicador de Posición Plana (PPI por sus siglas en inglés de Plan Position Indicator). Este modo primero fija la elevacion de la antena a un ángulo determinado con respecto a la horizontal y gira 360 grados con respecto al norte. Luego repite esta tarea míltiples veces a diferentes elevaciones hasta completar un Patrón de Cobertura Volumétrico (VCP por sus siglas en inglés de Volume Coverage Pattern) como se puede obervar en la siguiente animación del curso Fundamentos en radares meterológicos (The Commet program - meted.ucar.edu)
2. Red Nacional de Radares de Colombia
Colombia cuenta actualmente con una red de 12 radares meterológicos que hacen vigilancia continua a las condiciones de precipitación en el país como se puede observar en la siguiente tabla:
Entidad |
Número de radares |
banda |
---|---|---|
IDEAM |
4 |
C |
Aerocivil |
6 |
C |
SIATA |
1 |
C |
IDIGER |
1 |
X |
Librerías
A continuación vamos a importar las librerías que utilizaremos en este cuadernillo
import boto3
import botocore
import cartopy.crs as ccrs
import cmweather
import folium
import fsspec
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xradar as xd
from folium import plugins
0.3.0
2.1 Ubicación de los radares meteorológicos
Podemos usar folium
para crear un mapa dinámico donde visualizamos la ubicación de los radares meteorológicos en Colombia
min_lon, max_lon, min_lat, max_lat = -90, -72, -1, 14
map_ = folium.Map(
location=[8, -76],
zoom_start=6,
min_lat=min_lat,
max_lat=max_lat,
min_lon=min_lon,
max_lon=max_lon,
zoom_control=False,
control_scale=True,
scrollWheelZoom=False,
width=1000,
height=700,
)
folium.TileLayer("openstreetmap").add_to(map_)
folium.TileLayer("cartodbpositron").add_to(map_)
folium.TileLayer("cartodbdark_matter").add_to(map_)
folium.LayerControl().add_to(map_)
minimap = folium.plugins.MiniMap()
map_ = map_.add_child(minimap)
# map_
Ahora podemos agregar la ubicación de cada radar en el mapa
def locate_radar(row):
df = row.to_frame()
html = df.iloc[:-3].to_html(
classes="table table-striped table-hover table-condensed table-responsive"
)
popup = folium.Popup(html, max_width=2650)
folium.Marker(location=[row.lat, row.lon], popup=popup).add_to(map_)
rad = float(row["rad"])
circle = folium.vector_layers.Circle(
location=[row.lat, row.lon],
color=row.color,
fill=True,
fill_color=f"{row['color']}",
radius=rad * 2500,
fill_opacity=0.1,
)
circle.add_to(map_)
df_radares = pd.read_csv("../data/radar_locations.csv")
Aplicamos la función locate_radar
a cada fila del dataframe usando el método pd.apply
df_radares.apply(locate_radar, axis=1)
map_
2.2 Radares meteorológicos del IDEAM disponibles en AWS
El Instituto de Hidrología, Meteorología y Estudios Ambientales - IDEAM ha dispuesto de manera pública los datos del radar meteorológico.
La estructura del bucket es s3://s3-radaresideam/l2_data/YYYY/MM/DD/Radar_name/RRRAAMMDDTTTTTT.RAWXXXX dónde:
YYYY es el año de 4 dígitos
MM es el mes de 2 dígitos
DD es el mes de 2 dígitos día
Radar_name nombre del radar. Las opciones son Guaviare, Munchique, Barrancabermja y Carimagua
RRRAAMMDDTTTTTT es el nombre del archivo del radar con lo siguiente:
RRR las tres primeras letras del nombre del radar (p. ej., GUA para el radar Guaviare)
YY es el año de 2 dígitos
MM es el mes de 2 dígitos
DD es el día de 2 dígitos
TTTTTT es la hora en la que se realizó el escaneo (GTM)
RAWXXXX Formato de archivo Sigmet y código único proporcionado por el software IRIS
Veamos un ejemplo de cómo luce el repositorio usando AWS cli y el siguiente commando en la terminal aws s3 ls --no-sign-request s3://s3-radaresideam/
!aws s3 ls --no-sign-request s3://s3-radaresideam/l2_data/2021/09/19/
<botocore.awsrequest.AWSRequest object at 0x7fad91df9a10>
3. Acceso a los datos de radar usando Python
Podemos acceder a los datos de radar usando Python
y librerías como Py-Art
o Xradar
. Para leer los datos de radar primero necesitamos generar una conexión entre el repositorio de datos s3
de Amazon y nuestro computador. Esto lo podemos lograr mediante la librería boto3
# direccion del bucket
str_bucket = "s3://s3-radaresideam/"
# creacion de un objeto s3
s3 = boto3.resource(
"s3",
config=botocore.client.Config(
signature_version=botocore.UNSIGNED, user_agent_extra="Resource"
),
)
# Creamos un objeto s3.Bucket
bucket = s3.Bucket("s3-radaresideam")
# Generamos una consulta tipo
query = "l2_data/2022/08/09/Carimagua/CAR22080919"
# Listamos los archivos dentro de nuestro bucket
radar_files = [f"{str_bucket}{i.key}" for i in bucket.objects.filter(Prefix=f"{query}")]
radar_files[:5]
['s3://s3-radaresideam/l2_data/2022/08/09/Carimagua/CAR220809190003.RAWDSVV',
's3://s3-radaresideam/l2_data/2022/08/09/Carimagua/CAR220809190315.RAWDSW0',
's3://s3-radaresideam/l2_data/2022/08/09/Carimagua/CAR220809190401.RAWDSW3',
's3://s3-radaresideam/l2_data/2022/08/09/Carimagua/CAR220809190505.RAWDSW8',
's3://s3-radaresideam/l2_data/2022/08/09/Carimagua/CAR220809191003.RAWDSWM']
Esta forma es un poco complicada. Podemos utilizar fsspec
que permite hacer un “montaje” del bucket en nuestro computador y es mucho más simplificado.
fs = fsspec.filesystem("s3", anon=True)
files = sorted(fs.glob("s3-radaresideam/l2_data/2022/08/09/Carimagua/CAR22080919*"))
files[:5]
/usr/share/miniconda3/envs/atmoscol2023/lib/python3.11/site-packages/fsspec/registry.py:272: UserWarning: Your installed version of s3fs is very old and known to cause
severe performance issues, see also https://github.com/dask/dask/issues/10276
To fix, you should specify a lower version bound on s3fs, or
update the current installation.
warnings.warn(s3_msg)
['s3-radaresideam/l2_data/2022/08/09/Carimagua/CAR220809190003.RAWDSVV',
's3-radaresideam/l2_data/2022/08/09/Carimagua/CAR220809190315.RAWDSW0',
's3-radaresideam/l2_data/2022/08/09/Carimagua/CAR220809190401.RAWDSW3',
's3-radaresideam/l2_data/2022/08/09/Carimagua/CAR220809190505.RAWDSW8',
's3-radaresideam/l2_data/2022/08/09/Carimagua/CAR220809191003.RAWDSWM']
Una vez que tengamos nuestros archivos de interés, podemos cargar uno en memoria usando Xradar
o Py-Art
. Por ejemplo, estamos interesados en el octavo archivo (index=7
).
3.1 Lectura de datos usando Xradar
Los datos de radar de IDEAM también se puden leer usando xradar.io.open_iris_datatree
. Esta librería, a diferencia de Py-Art
usa como base Xarray
retornando un xarray.Dataset
en caso de seleccionar un barrido en particular o un xarray.datatree
en caso de leer todas las elevaciones contenidas dentro del mismo archivo.
# preparamos un archivo para abrir de manera local
task_file = fsspec.open_local(
f"simplecache::s3://{files[7]}", s3={"anon": True}, filecache={"cache_storage": "."}
)
# leemos nuestro archivo usando xradar
radar = xd.io.open_iris_datatree(task_file)
display(radar)
<xarray.DatasetView> Dimensions: () Data variables: volume_number int64 0 platform_type <U5 'fixed' instrument_type <U5 'radar' time_coverage_start <U20 '2022-08-09T19:15:05Z' time_coverage_end <U20 '2022-08-09T19:16:21Z' longitude float64 -71.33 altitude float64 206.0 latitude float64 4.564 Attributes: Conventions: None version: None title: None institution: None references: None source: None history: None comment: im/exported using xradar instrument_name: None
Este datatree
solo contine una elevación. Podemos acceder usando notación similar a la de un diccionario
display(radar["sweep_0"])
<xarray.DatasetView> Dimensions: (azimuth: 720, range: 994) Coordinates: elevation (azimuth) float32 ... time (azimuth) datetime64[ns] 2022-08-09T19:15:56.891000 ..... * range (range) float32 1e+03 1.3e+03 ... 2.986e+05 2.989e+05 longitude float64 ... latitude float64 ... altitude float64 ... * azimuth (azimuth) float32 0.03571 0.5795 1.019 ... 359.0 359.6 Data variables: (12/17) DBTH (azimuth, range) float32 ... DBZH (azimuth, range) float32 ... VRADH (azimuth, range) float32 ... WRADH (azimuth, range) float32 ... ZDR (azimuth, range) float32 ... KDP (azimuth, range) float32 ... ... ... DB_DBZE8 (azimuth, range) float32 ... sweep_mode <U20 ... sweep_number int64 ... prt_mode <U7 ... follow_mode <U7 ... sweep_fixed_angle float64 ...
3.2 Gráfico de reflectividad
Podemos usar la funcionalidad xarray.plot
que se encuentra incorporada en xarray
.
radar["sweep_0"]["DBZH"].plot(cmap="ChaseSpectral", vmin=-10, vmax=60)
<matplotlib.collections.QuadMesh at 0x7f7733468410>
Como podemos observar en el gráfico anterior los datos se encuentran ordenados por las dimensiones y coordenadas azimuth
y range
. Podemos usar la xr.georeference
para generar las coordenadas polares a nuestro dataset
radar = radar.xradar.georeference()
display(radar["sweep_0"])
<xarray.DatasetView> Dimensions: (azimuth: 720, range: 994) Coordinates: elevation (azimuth) float64 0.4779 0.4779 0.4779 ... 0.4779 0.4779 time (azimuth) datetime64[ns] 2022-08-09T19:15:56.891000 ..... * range (range) float32 1e+03 1.3e+03 ... 2.986e+05 2.989e+05 longitude float64 -71.33 latitude float64 4.564 altitude float64 206.0 * azimuth (azimuth) float32 0.03571 0.5795 1.019 ... 359.0 359.6 crs_wkt int64 0 x (azimuth, range) float64 0.6231 0.8101 ... -2.191e+03 y (azimuth, range) float64 999.9 1.3e+03 ... 2.987e+05 z (azimuth, range) float64 214.4 216.9 ... 7.948e+03 Data variables: (12/17) DBTH (azimuth, range) float32 ... DBZH (azimuth, range) float32 ... VRADH (azimuth, range) float32 ... WRADH (azimuth, range) float32 ... ZDR (azimuth, range) float32 ... KDP (azimuth, range) float32 ... ... ... DB_DBZE8 (azimuth, range) float32 ... sweep_mode <U20 ... sweep_number int64 ... prt_mode <U7 ... follow_mode <U7 ... sweep_fixed_angle float64 ...
Como podemos ver x
, y
y z
fueron agregados como nuevos campos de coordenadas. Ahora podemos proceder a generar el gráfico pero esta vez de manera georeferenciada
radar["sweep_0"]["DBZH"].plot(x="x", y="y", cmap="ChaseSpectral", vmin=-10, vmax=60)
<matplotlib.collections.QuadMesh at 0x7f773229fdd0>
Para generar gráficas de radar en un sistema coordenado de referencia podemos usar cartopy
. Extraemos el sistema de referencia desde el objeto radar usando xd.georeference.get_crs
proj_crs = xd.georeference.get_crs(radar["sweep_0"].ds)
cart_crs = ccrs.Projection(proj_crs)
Ahora sí procedemos a generar el gráfico de reflectividad en un sistema coordenado de referencia
fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
radar["sweep_0"]["DBZH"].plot(
x="x",
y="y",
cmap="ChaseSpectral",
transform=cart_crs,
cbar_kwargs=dict(pad=0.075, shrink=0.75),
vmin=-10,
vmax=60,
)
ax.coastlines()
ax.gridlines(draw_labels=True, ls="--", lw=0.5)
/usr/share/miniconda3/envs/atmoscol2023/lib/python3.11/site-packages/cartopy/mpl/geoaxes.py:1781: UserWarning: The input coordinates to pcolormesh are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolormesh.
result = super().pcolormesh(*args, **kwargs)
<cartopy.mpl.gridliner.Gridliner at 0x7f7732c2c910>
3.3 Gráficos interactivos usando hvplot
Los gráficos estáticos pueden ser usados principalmente para publicaciones donde la interactividad no es necesaria. Sin embargo, en cursos como el taller de programación de AtmosCol 2023 podemos utilizar las ventajas de la interactividad para aprovechar mejor la información disponible.
import hvplot.xarray
ref = radar["sweep_0"].DBZH.hvplot.quadmesh(
x="x",
y="y",
cmap="ChaseSpectral",
clabel="Horizontal Reflectivity (dBZ)",
title=f'Horizontal Reflectivity \n {radar.attrs["instrument_name"]} Radar',
clim=(-20, 60),
height=500,
rasterize=True,
width=600,
)
ref
Ahora podemos interacturar con más de un plot a la vez
zdr = radar["sweep_0"].ZDR.hvplot.quadmesh(
x="x",
y="y",
cmap="jet",
clabel="Differential Reflectivity (dB)",
title=f'Horizontal Reflectivity \n {radar.attrs["instrument_name"]} Radar',
clim=(-1, 6),
height=500,
rasterize=True,
width=600,
)
(ref + zdr).cols(1)
4. Estimados de lluvia utilizando la ecuación de Marshall and Gunn (1953)
Uno de los productos de mayor interés para la comunidad científica, y en general para todos los usuarios de la información de radares meteorológicos, es la estimación de lluvia. Desde 1948 los científicos han tratado de convertir las mediciones de reflectividad (\(Z\)) en mediciones de intensidad de precipitación (\(R\)). Para esto se han desarrollado múltiples ecuaciones empíricas \(Z-R\). En 1953, Marshall and Gunn encontraron una de las relaciones más comunes e incluso utilizadas en la actualidad, que sigue la forma:
Despejando para R tenemos:
Cuidado!!!
La reflectividad del radar esta dada en unidades decibélicas es decir
\(Z [dBZ] = 10 * log_{10}(Z[mm^{6}m^{3}])\)
Por lo tanto debemos convertir a unidades lineales
\(Z[mm^{6}m^{3}] = 10 ^{\left(\frac{Z[dBZ]}{10}\right)}\)
# Reflectividad en unidades decibelicas
ref_log = radar["sweep_0"].DBZH
# Reflectividad en unidades lineales
ref_lin = 10 ** (ref_log / 10)
# Estimado de lluvia usando Marshall and Gunn (1953)
rr = (1 / 200) ** (1 / 1.6) * ref_lin ** (1 / 1.6)
Ahora podemos visualizar el campo de lluvia generado
fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
rr.plot(
x="x",
y="y",
cmap="jet",
transform=cart_crs,
cbar_kwargs=dict(
pad=0.075, shrink=0.75, label=r"$Intensidad \ de \ lluvia \ [mm hr^{-1}]$"
),
vmin=0,
vmax=50,
)
ax.set_title(r"$Estimado \ de \ lluvia$")
ax.coastlines()
ax.gridlines(draw_labels=True, ls="--", lw=0.5)
/usr/share/miniconda3/envs/atmoscol2023/lib/python3.11/site-packages/cartopy/mpl/geoaxes.py:1781: UserWarning: The input coordinates to pcolormesh are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolormesh.
result = super().pcolormesh(*args, **kwargs)
<cartopy.mpl.gridliner.Gridliner at 0x7f7715f03b90>
Conclusiones
En este cuadernillo aprendimos sobre la red nacional de radares meteorológicos de Colombia. Logramos acceder a los datos del repositorio de AWS usando Python. Utilizamos Py-Art y Xradar para acceder, visualizar los datos de manera intectiva usando hvplot y generar estimados de precipitación utilzando la ecuación de Marshall and Gunn (1953).
Resources and references
Radar Cookbook [https://projectpythia.org/radar-cookbook/README.html] DOI[https://doi.org/10.5281/zenodo.8075855]
Rose, B. E. J., Kent, J., Tyle, K., Clyne, J., Banihirwe, A., Camron, D., May, R., Grover, M., Ford, R. R., Paul, K., Morley, J., Eroglu, O., Kailyn, L., & Zacharias, A. (2023). Pythia Foundations (Version v2023.05.01) https://doi.org/10.5281/zenodo.7884572
Marshall, J. S., and W. M. K. Palmer, 1948: THE DISTRIBUTION OF RAINDROPS WITH SIZE. J. Atmos. Sci., 5, 165–166, https://doi.org/10.1175/1520-0469(1948)005<0165:TDORWS>2.0.CO;2.
Helmus, J.J. & Collis, S.M., (2016). The Python ARM Radar Toolkit (Py-ART), a Library for Working with Weather Radar Data in the Python Programming Language. Journal of Open Research Software. 4(1), p.e25. DOI: http://doi.org/10.5334/jors.119
Xradar Python library [https://docs.openradarscience.org/projects/xradar/en/stable/#] DOI[https://docs.openradarscience.org/projects/xradar/en/stable/#]