Clase 11: Manejo de archivos FITS#
FITS (Flexible Image Transport System) es el formato de archivo estándar más utilizado en astronomía para el almacenamiento, la transmisión y la manipulación de datos científicos. Diseñado específicamente para ser independiente de la plataforma y duradero a largo plazo, un archivo FITS consta de una o más secciones que incluyen un encabezado (header) en texto legible con metadatos en formato key-value, seguido de datos binarios que pueden representar imágenes multidimensionales, espectros o tablas de datos. Este formato permite incluir información sobre las condiciones de observación, las coordenadas celestiales y la calibración.
Visualizando imágenes en archivos FITS#
Imagen procesada del disco de escombros de Fomalhaut basada en imágenes tomadas por el Hubble Space Telescope. https://github.com/saint-germain/cazandoplanetas
from astropy.io import fits
fomalhaut=fits.open("https://github.com/saint-germain/cazandoplanetas/raw/66fb7ce82feb5e983e465b7cfd1a638d8bc3f038/Intermedio/PCA_FOMAL_full_centroids.fits")
fomalhaut.info()
Filename: /Users/germanchaparro/.astropy/cache/download/url/babdce75c785d3e44afc57fb60e3e07d/contents
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 1250 (2000, 2000) float32
hst_img=fomalhaut[0].data
hst_img
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]], shape=(2000, 2000), dtype='>f4')
hst_img.max()
np.float32(0.9672259)
plt.imshow(hst_img)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 1
----> 1 plt.imshow(hst_img)
NameError: name 'plt' is not defined
plt.figure(figsize=(10,10))
plt.imshow(hst_img*(hst_img>-0.0045)*(hst_img<0.004),cmap="gist_heat")
plt.scatter(1173,1205,s=1000,facecolors='none', edgecolors='b',lw=2)
plt.xlim(900,1250)
plt.ylim(1100,1400)
Explorando un cubo de datos en un archivo FITS#
Datos de espectros de CO galáctico (CO J = 1–0 line, vrest = 115.27 GHz) tomados por el radiotelescopio MINI (Crédito: Leonardo Bronfman, U. de Chile).
!wget https://github.com/saint-germain/astropy/raw/refs/heads/master/southgal_fixbadc.fits
# esta línea da un warning debido a que es un fits no estándar
cubo = fits.open("southgal_fixbadc.fits") #abrir objeto cubo de datos
cubo.info()
cubo[0].data
cubo[0].header
def values(h,j):
# construye los ejes en velocidad, longitud y latitud con la información en el header
N=h['NAXIS'+str(j)];
val=np.zeros(N);
for i in range(0,N):
val[i] = (i+1-float(h['CRPIX'+str(j)]))*float(h['CDELT'+str(j)]) + float(h['CRVAL'+str(j)]);
return val;
data = cubo[0].data #extraer matriz de datos
header = cubo[0].header #extraer el header del archivo fits
#print header
#Estos seran los tres arreglos con los valores reales de los tres ejes del cubo
velocidad=values(header,1)
longitud=values(header,2)
latitud=values(header,3)
# Escogemos algún valor arbitrario para longitud, latitud
# para extraer el espectro para ese apuntamiento
i_l=-1
i_b=-1
T = data[i_b][i_l][:]
lon=longitud[i_l]
lat=latitud[i_b]
plt.plot(velocidad,T)
plt.title(r'Espectro de la emisión de CO $J=1-0$ para ($l,b$)=(%i,%i)'%(lon,lat))
maxeix=np.argmax(data,axis=2)
varr=np.array([velocidad[i] for i in maxeix.ravel()]).reshape(maxeix.shape)
plt.figure(figsize=(20,5))
plt.imshow(varr,cmap='bwr',extent=[longitud[0],longitud[-1],latitud[0],latitud[-1]])
También existe una librería llamada spectral_cube que puede procesar automáticamente este tipo de cubos de datos en formato fits (estandarizados): https://learn.astropy.org/tutorials/FITS-cubes.html
Tablas (Votable) y WCS#
VOTABLE es un estándar de formato basado en XML diseñado específicamente para el intercambio de datos tabulares dentro de la comunidad astronómica, bajo la supervisión de la IVOA (International Virtual Observatory Alliance). Un archivo VOTABLE permite incluir metadatos enriquecidos que describen las unidades, los sistemas de coordenadas y la procedencia de la información.
Exploración de datos xml en formato Votable: el caso de las estrellas en el grupo móvil Lower Centaurus Crux (LCC).
Ir a SIMBAD, buscar “lcc”, children objects, descargar Votable
Enlace: http://simbad.u-strasbg.fr/simbad/sim-id?Ident=NAME%20LCC&NbIdent=query_hlinks&Coord=12%2019-57.1&parents=1&children=1938&submit=children&siblings=19&hlinksdisplay=h_all
!wget https://github.com/saint-germain/astropy/raw/refs/heads/master/lcc_simbad_votable.xml
from astropy.io.votable import parse_single_table
votable = parse_single_table("lcc_simbad_votable.xml").to_table()
votable
fig = plt.figure(figsize=(10, 10))
ax = plt.subplot()
plt.scatter(votable['RA_d'],votable['DEC_d'])
plt.xlabel(r'RA')
plt.ylabel(r'Dec')
Las coordenadas mundiales (World Coordinates) sirven para ubicar una medida en un espacio de parámetros multidimensional. Un Sistema de Coordenadas Mundiales (WCS en inglés) especifica las coordenadas físicas, o mundiales (WC en inglés), que se adjuntan a cada píxel o spaxel de una imagen o matriz de \(N\) dimensiones. Se ha desarrollado un elaborado conjunto de normas y convenciones para el formato FITS (Flexible Image Transport System) (Wells et al. 1981). Un ejemplo típico de WCS es la especificación de la Ascensión Recta (AR) y la Declinación (Dec) en el cielo asociada a una determinada ubicación de píxel o vóxel en una imagen celeste bidimensional (Greisen y Calabretta 2002; Calabretta y Greisen 2002).
Hay dos formas principales de inicializar un objeto WCS: con un diccionario de Python (o un objeto tipo diccionario, como la cabecera de un archivo FITS) o con listas de Python. En este ejemplo, inicializaremos un objeto astropy.wcs.WCS con dos dimensiones, como sería necesario para representar una imagen.
Adaptado y traducido de https://learn.astropy.org/tutorials/celestial_coords1.html
from astropy.wcs import WCS
# world coordinate system https://learn.astropy.org/tutorials/celestial_coords1.html
wcs_input_dict = {
'CTYPE1': 'RA',
'CUNIT1': 'deg',
'CTYPE2': 'DEC',
'CUNIT2': 'deg'
}
wcs_helix = WCS(wcs_input_dict)
wcs_helix
fig = plt.figure(figsize=(10, 10))
ax = plt.subplot(projection=wcs_helix)
plt.scatter(votable['RA_d'],votable['DEC_d'])
plt.xlabel(r'RA')
plt.ylabel(r'Dec')
overlay = ax.get_coords_overlay('icrs')
overlay.grid(color='black', ls='dotted')
fig = plt.figure(figsize=(10, 10))
ax = plt.subplot(projection=wcs_helix)
plt.scatter(votable['RA_d'],votable['DEC_d'])
plt.xlabel(r'RA')
plt.ylabel(r'Dec')
overlay = ax.get_coords_overlay('galactic')
overlay.grid(color='black', ls='dotted')
Haciendo un mapa del gas y el polvo en la Nube Menor de Magallanes#
Polvo: Encontrar y Descargar una Imagen de Herschel#
Vamos a ver cómo aparece la Nube Menor de Magallanes (SMC) con observaciones de emisión a 350 micras de Herschel para rastrear la emisión térmica (en continuo del polvo). Esto se puede hacer con astroquery.
Podemos consultar los datos por misión, dar un vistazo rápido a la tabla de resultados y descargar los datos después de seleccionar una longitud de onda o filtro específico.
Como estamos buscando datos de Herschel provenientes de una misión de la ESA, utilizaremos la clase astroquery.ESASky.
Específicamente, el método ESASKY.query_region_maps() nos permite buscar una región específica del cielo utilizando ya sea un objeto SkyCoord de Astropy o una cadena de texto que especifique el nombre de un objeto. En este caso, podemos simplemente buscar la SMC. También se puede especificar un radio de búsqueda alrededor del objeto.
Tutorial original por Dhanesh Krishnarao (DK), Shravan Shetty, Diego Gonzalez-Casanova, Audra Hernandez, Kris Stern, Kelle Cruz, Stephanie Douglas.
!pip install astroquery
from astroquery.esasky import ESASky
from astropy import units as u
import numpy as np
from astroquery.utils import TableList
import matplotlib.pyplot as plt
%matplotlib inline
maps_list = ESASky.list_maps()
print(maps_list)
# Query for Herschel data in a 1 degree radius around the SMC
result = ESASky.query_region_maps("SMC", radius=1 * u.deg, missions="Herschel")
print(result)
# Query for Herschel data in a 1 degree radius around the SMC
result = ESASky.query_region_maps("SMC", radius=1 * u.deg, missions="Herschel")
print(result)
Aquí, el resultado es un TableList que contiene 24 productos de datos de Herschel que pueden descargarse. Podemos ver qué información está disponible en este TableList examinando las keys de la tabla de Herschel.
result["HERSCHEL"].keys()
Queremos encontrar una imagen de 350 micras, así que necesitamos examinar con más detalle los filtros utilizados para estas observaciones.
result["HERSCHEL"]["filter"]
Afortunadamente para nosotros, hay una observación realizada con tres filtros: 250, 350 y 500 micras. Este es el objeto que queremos descargar. Una forma de hacerlo es creando una máscara booleana para seleccionar la entrada de la tabla correspondiente al filtro deseado. Luego, el método ESASky.get_maps() descargará nuestros datos proporcionando un argumento de tipo TableList.
La siguiente descarga puede tomar varios minutos.
filters = result["HERSCHEL"]["filter"].astype(
str
) # Convert the list of filters from the query to a string
# Construct a boolean mask, searching for only the desired filters
mask = np.array(["250, 350, 500" == s for s in filters], dtype="bool")
# Re-construct a new TableList object containing only our desired query entry
target_obs = TableList(
{"HERSCHEL": result["HERSCHEL"][mask]}
) # This will be passed into ESASky.get_maps()
IR_images = ESASky.get_maps(target_obs) # Download the images
IR_images["HERSCHEL"][0][
"350"
].info() # Display some information about the 350 micron image
Vamos a extraer únicamente la información de WCS y datos de la imagen de 350 micrones. Para obtener un objeto WCS simplemente pasamos a la función WCS() el objeto FITS como argumento.
herschel_header = IR_images["HERSCHEL"][0]["350"]["image"].header
herschel_wcs = WCS(IR_images["HERSCHEL"][0]["350"]["image"]) # Extract WCS information
herschel_imagehdu = IR_images["HERSCHEL"][0]["350"]["image"] # Extract Image data
print(herschel_wcs)
Con esto, podemos mostrar esta imagen utilizando matplotlib con WCSAxes y el objeto LogNorm() para aplicar una escala logarítmica a nuestra imagen.
from matplotlib.colors import LogNorm
# Set Nans to zero
himage_nan_locs = np.isnan(herschel_imagehdu.data)
herschel_data_nonans = herschel_imagehdu.data
herschel_data_nonans[himage_nan_locs] = 0
# Initiate a figure and axis object with WCS projection information
fig = plt.figure(figsize=(18, 12))
ax = fig.add_subplot(111, projection=herschel_wcs)
# Display the moment map image
im = ax.imshow(herschel_data_nonans, cmap="viridis", norm=LogNorm(vmin=2, vmax=50))
# ax.invert_yaxis() # Flips the Y axis
# Add axes labels
ax.set_xlabel("Right Ascension", fontsize=16)
ax.set_ylabel("Declination", fontsize=16)
ax.grid(color="white", ls="dotted", lw=2)
# Add a colorbar
cbar = plt.colorbar(im, pad=0.07)
cbar.set_label(
"".join(["Herschel 350" r"$\mu$m ", "(", herschel_header["BUNIT"], ")"]), size=16
)
# Overlay set of Galactic Coordinate Axes
overlay = ax.get_coords_overlay("galactic")
overlay.grid(color="black", ls="dotted", lw=1)
overlay[0].set_axislabel("Galactic Longitude", fontsize=14)
overlay[1].set_axislabel("Galactic Latitude", fontsize=14)
Ejercicio: Con los contenidos del catálogo de nubes oscuras de Barnard VII/220A/barnard haga un mosaico con los objetos del catálogo que tengan imágenes de Herschel disponibles.
Descargar y visualizar datos de HI de la SMC#
Utilizaremos datos de emisión HI de 21 cm del survey HI4Pi. Queremos observar la emisión de gas neutro de la Nubes Menor de Magallanes y aprender sobre la cinemática del sistema y las densidades de columna. El cubo de datos se encuentra completo en VizieR. Archivo relevante disponible vía ftp desde CDS Strasbourg.
Tenemos una versión reducida del mismo, que será un cubo de datos FITS en coordenadas galácticas utilizando la proyección tangencial del cielo.
!pip install spectral-cube
from astropy.utils.data import download_file
from astropy.io import fits # We use fits to open the actual data file
from spectral_cube import SpectralCube
import astropy.units as u
# Downloads the HI data in a fits file format
hi_datafile = download_file(
"http://data.astropy.org/tutorials/FITS-cubes/reduced_TAN_C14.fits",
cache=True,
show_progress=True,
)
El paquete spectral_cube realiza gran parte del trabajo detallado necesario para manipular estos datos e incluso explorarlos rápidamente. Esto requiere que el archivo FITS esté estandarizado, lo cual a veces puede no ocurrir.
En este caso el archivo fue producido de forma profesional, por lo que podemos abrir nuestro archivo de datos y leer la información como un SpectralCube.
La variable cube contiene los datos utilizando SpectralCube, y hi_data es el cubo de datos del archivo FITS sin el formato especial proporcionado por SpectralCube.
hi_data = fits.open(hi_datafile) # Open the FITS file for reading
cube = SpectralCube.read(hi_data) # Initiate a SpectralCube
hi_data.close() # Close the FITS file - we already read it in and don't need it anymore!
print(cube)
Algunas cosas a las que debemos prestar atención aquí:#
Un cubo de datos tiene tres ejes. En este caso, hay Longitud Galáctica (x), Latitud Galáctica (y) y un eje espectral en términos de una velocidad LSR (z — listado como s en spectral_cube).
Los datos contenidos en el cubo viven como un ndarray con forma (n_s, n_y, n_x), de manera que el eje 0 corresponde al eje espectral, el eje 1 corresponde al eje de Latitud Galáctica y el eje 2 corresponde al eje de Longitud Galáctica.
Cuando hacemos print(cube) podemos ver la forma, tamaño y unidades de todos los ejes, así como los datos almacenados en el cubo. En este cubo, las unidades de los datos son temperaturas (K). Los ejes espaciales están en grados y el eje espectral está en (metros / segundo).
El cubo también contiene información sobre las coordenadas correspondientes a los datos en forma de un objeto WCS (World Coordinate System).
SpectralCube es inteligente y mantiene todos los datos enmascarados hasta que realmente los necesites, de modo que puedas trabajar con grandes conjuntos de datos. ¡Así que veamos cómo lucen realmente nuestros datos!
SpectralCube tiene un método quicklook() que puede proporcionar una vista previa rápida y conveniente de los datos. Es útil cuando solo necesitas echar un vistazo a un corte o espectro sin conocer ninguna otra información (por ejemplo, para asegurarte de que los datos no estén corruptos o que correspondan a la región correcta).
Para hacer esto, indexamos nuestro cubo a lo largo de un eje (para un corte) o dos ejes (para un espectro):
cube[
300, :, :
].quicklook() # Slice the cube along the spectral axis, and display a quick image
cube[:, 75, 75].quicklook() # Extract a single spectrum through the data cube
Crear un cubo más pequeño, enfocándonos en la SMC#
El cubo de datos HI que descargamos es más grande de lo que realmente necesitamos. Intentemos acercarnos únicamente a la parte que nos interesa y crear un nuevo sub_cube.
La forma más sencilla de hacerlo es recortando una parte del cubo utilizando índices o coordenadas.
Podemos extraer las coordenadas del mundo desde el cubo utilizando el método .world().
Advertencia: usar .world() extraerá coordenadas de cada posición que solicites. Esto puede representar una CANTIDAD ENORME de datos si no haces un corte del cubo. Una solución alternativa es hacer un corte a lo largo de dos ejes y extraer coordenadas únicamente a lo largo de una sola dimensión.
La salida de .world() es un objeto Quantity de Astropy que representa las coordenadas de los píxeles e incluye unidades. Puedes extraer estos objetos Quantity de Astropy utilizando cortes (slicing) sobre los datos.
_, lat, _ = cube.world[0, :, 0] # extract latitude world coordinates from cube
_, _, lon = cube.world[0, 0, :] # extract longitude world coordinates from cube
Ahora podemos extraer un sub_cube en las coordenadas espaciales del cubo.
# Define desired latitude and longitude range
lat_range = [-46, -40] * u.deg
lon_range = [306, 295] * u.deg
# Create a sub_cube cut to these coordinates
sub_cube = cube.subcube(
xlo=lon_range[0], xhi=lon_range[1], ylo=lat_range[0], yhi=lat_range[1]
)
print(sub_cube)
Realmente no necesitamos datos de un rango de velocidades tan grande, así que extraigamos solamente una pequeña porción. Podemos hacer esto en cualquier unidad que queramos utilizando el método .spectral_slab().
sub_cube_slab = sub_cube.spectral_slab(-300.0 * u.km / u.s, 300.0 * u.km / u.s)
print(sub_cube_slab)
Mapas de Momentos#
Los mapas de momentos son una herramienta de análisis muy útil para estudiar cubos de datos. En resumen, un momento es una integral ponderada a lo largo de un eje (típicamente el eje espectral) que puede proporcionar información sobre la intensidad total (o densidad de columna), la velocidad media o la dispersión de velocidades a lo largo de las líneas de visión.
SpectralCube hace esto muy sencillo mediante el método .moment().
Podemos convertir a unidades espectrales más convenientes, como km/s, y estas nuevas proyecciones 2D también pueden guardarse como nuevos archivos FITS, incluyendo información WCS modificada.
moment_0 = sub_cube_slab.with_spectral_unit(u.km / u.s).moment(
order=0
) # Zero-th moment
moment_1 = sub_cube_slab.with_spectral_unit(u.km / u.s).moment(order=1) # First moment
# Write the moments as a FITS image
# moment_0.write('hi_moment_0.fits')
# moment_1.write('hi_moment_1.fits')
print("Moment_0 has units of: ", moment_0.unit)
print("Moment_1 has units of: ", moment_1.unit)
# Convert Moment_0 to a Column Density assuming optically thin media
hi_column_density = moment_0 * 1.82 * 10**18 / (u.cm * u.cm) * u.s / u.K / u.km
Recordemos que el framework WCSAxes de Astropy nos permite mostrar imágenes con diferentes ejes de coordenadas y proyecciones.
Siempre que tengamos un objeto WCS asociado a los datos, podemos transferir esa proyección a un eje de matplotlib. SpectralCube permite acceder directamente al objeto WCS asociado a un cubo de datos.
print(moment_1.wcs) # Examine the WCS object associated with the moment map
Como era de esperarse, la imagen del primer momento que creamos solo tiene dos ejes (Longitud Galáctica y Latitud Galáctica). Podemos pasar este objeto WCS directamente a una instancia de eje de matplotlib.
import matplotlib.pyplot as plt
# Initiate a figure and axis object with WCS projection information
fig = plt.figure(figsize=(18, 12))
ax = fig.add_subplot(111, projection=moment_1.wcs)
# Display the moment map image
im = ax.imshow(moment_1.hdu.data, cmap="RdBu_r", vmin=0, vmax=200)
ax.invert_yaxis() # Flips the Y axis
# Add axes labels
ax.set_xlabel("Galactic Longitude (degrees)", fontsize=16)
ax.set_ylabel("Galactic Latitude (degrees)", fontsize=16)
# Add a colorbar
cbar = plt.colorbar(im, pad=0.07)
cbar.set_label("Velocity (km/s)", size=16)
# Overlay set of RA/Dec Axes
overlay = ax.get_coords_overlay("fk5")
overlay.grid(color="white", ls="dotted", lw=2)
overlay[0].set_axislabel("Right Ascension (J2000)", fontsize=16)
overlay[1].set_axislabel("Declination (J2000)", fontsize=16)
# Overplot column density contours
levels = (1e20, 5e20, 1e21, 3e21, 5e21, 7e21, 1e22) # Define contour levels to use
ax.contour(hi_column_density.hdu.data, cmap="Greys_r", alpha=0.5, levels=levels)
Superponer contornos de HI de 21 cm sobre la Imagen Infrarroja de 350 Micras de Herschel#
Para comparar visualmente el gas neutro y el polvo trazados mediante emisión HI de 21 cm y emisión infrarroja de 350 micras, podemos utilizar contornos e imágenes en escala de color producidas con el framework WCSAxes y el método .get_transform().
El método WCSAxes.get_transform() devuelve una transformación desde un sistema de referencia especificado hacia las coordenadas de píxel/datos. Este método acepta una cadena de texto que especifique el sistema de referencia o un objeto WCS.
# Initiate a figure and axis object with WCS projection information
fig = plt.figure(figsize=(18, 12))
ax = fig.add_subplot(111, projection=herschel_wcs)
# Display the moment map image
im = ax.imshow(
herschel_data_nonans, cmap="viridis", norm=LogNorm(vmin=5, vmax=50), alpha=0.8
)
# ax.invert_yaxis() # Flips the Y axis
# Add axes labels
ax.set_xlabel("Right Ascension", fontsize=16)
ax.set_ylabel("Declination", fontsize=16)
ax.grid(color="white", ls="dotted", lw=2)
# Extract x and y coordinate limits
x_lim = ax.get_xlim()
y_lim = ax.get_ylim()
# Add a colorbar
cbar = plt.colorbar(im, fraction=0.046, pad=-0.1)
cbar.set_label(
"".join(["Herschel 350" r"$\mu$m ", "(", herschel_header["BUNIT"], ")"]), size=16
)
# Overlay set of RA/Dec Axes
overlay = ax.get_coords_overlay("galactic")
overlay.grid(color="black", ls="dotted", lw=1)
overlay[0].set_axislabel("Galactic Longitude", fontsize=14)
overlay[1].set_axislabel("Galactic Latitude", fontsize=14)
hi_transform = ax.get_transform(
hi_column_density.wcs
) # extract axes Transform information for the HI data
# Overplot column density contours
levels = (2e21, 3e21, 5e21, 7e21, 8e21, 1e22) # Define contour levels to use
ax.contour(
hi_column_density.hdu.data,
cmap="Greys_r",
alpha=0.8,
levels=levels,
transform=hi_transform,
) # include the transform information with the keyword "transform"
# Overplot velocity image so we can also see the Gas velocities
im_hi = ax.imshow(
moment_1.hdu.data,
cmap="RdBu_r",
vmin=0,
vmax=200,
alpha=0.5,
transform=hi_transform,
)
# Add a second colorbar for the HI Velocity information
cbar_hi = plt.colorbar(im_hi, orientation="horizontal", fraction=0.046, pad=0.07)
cbar_hi.set_label("HI " r"$21$cm Mean Velocity (km/s)", size=16)
# Apply original image x and y coordinate limits
ax.set_xlim(x_lim)
ax.set_ylim(y_lim)
Archivos HDF y observación terrestre#
HDF (Hierarchical Data Format) es un modelo de datos y formato de archivo versátil diseñado para gestionar y almacenar conjuntos de datos masivos y complejos. En los contextos de Astronomía y Observación de la Tierra, es ampliamente utilizado por su capacidad tipo “contenedor” para agrupar arreglos multidimensionales (como cubos de datos, cubos de simulación, y observaciones con múltiples capas), metadatos y tablas relacionadas dentro de un único archivo autodescriptivo y autocontenido.
A diferencia de formatos más simples, HDF permite operaciones de entrada/salida de alto rendimiento y el uso de chunking, que es una operación que permite el acceso a subconjuntos específicos de un conjunto de datos de escala TB —como una región geográfica particular o una línea espectral específica— sin necesidad de cargar todo el archivo en memoria.
Es el formato principal utilizado por grandes misiones como el sistema EOS (Earth Observing System) de la NASA y también se usa en simulaciones cosmológicas de gran escala, donde el volumen de datos y la complejidad estructural superan los límites de los archivos FITS tradicionales.
NASA Earthdata#
El portal NASA Earthdata sirve como la principal puerta de acceso a datos gratuitos sobre ciencias de la Tierra, administrado por el sistema EOSDIS (Earth Observing System Data and Information System). Los datos disponibles incluyen composición atmosférica, oceanografía, biología terrestre y estudios de la criosfera.
Herramientas como Earthdata Search y Worldview, permiten realizar consultas espaciales y temporales, extraer subconjuntos de enormes archivos HDF para regiones específicas de interés e interactuar con imágenes satelitales casi en tiempo real.
Al centralizar datos provenientes de distintos DAACs (Distributed Active Archive Centers), el portal facilita la investigación interdisciplinaria y garantiza que parámetros como el Vapor de Agua Precipitable (PWV) y la profundidad óptica de aerosoles (AOD), estén disponibles de manera abierta para el modelado climático, el monitoreo global y otras aplicaciones.
Para este ejercicio no es necesario, pero es recomendado crear su propia cuenta en el portal, desde la cual se puede crear un token de acceso para hacer descargas usando la API de Earthdata.
Medición satelital del vapor de agua precipitable#
El vapor de agua en la atmósfera terrestre es la principal fuente de opacidad e inestabilidad de fase en longitudes de onda milimétricas y submilimétricas. Las moléculas de agua poseen numerosas líneas de transición rotacionales y vibracionales que absorben la radiación cósmica entrante, “cerrando” efectivamente las ventanas atmosféricas necesarias para la observación.
El PWV (Precipitable Water Vapor, o Vapor de Agua Precipitable) se utiliza como una métrica en la caracterización de sitios para observaciones en estas longitudes de onda. Un valor más bajo de PWV se correlaciona directamente con una mayor transparencia atmosférica.
Por lo tanto, cuantificar el PWV mediano y su estabilidad estacional ayuda a determinar si un sitio de gran altitud puede satisfacer los exigentes requerimientos de la instrumentación submilimétrica.
Los instrumentos MODIS (Moderate Resolution Imaging Spectroradiometer) a bordo de los satélites Terra y Aqua (en LEO) hacen estimaciones diarias del Vapor de Agua Precipitable (PWV, Precipitable Water Vapor) utilizando mediciones tanto en el infrarrojo cercano (NIR) como en el infrarrojo térmico (IR).
Al comparar la emisión de bandas de absorción del vapor de agua con bandas de “ventana” no absorbentes, MODIS puede recuperar el contenido total de vapor de agua en la columna atmosférica sobre tierra y océano con alta resolución espacial.
!pip install netCDF4 pyhdf cartopy
A continuación usaremos la API de Earthdata para descargar los datos mensuales promediados tomados por los instrumentos MODIS para toda la Tierra en una malla de 1x1 grado, a partir del producto MYD08_M3. Interactuaremos con Earthdata a través de la librería requests, autenticándonos con un Token personalizado, generado desde el portal Earthdata.
En este enlace hay una descripción del producto.
# Se recomienda generar su propio token, este perderá validez
TOKEN='eyJ0eXAiOiJKV1QiLCJvcmlnaW4iOiJFYXJ0aGRhdGEgTG9naW4iLCJzaWciOiJlZGxqd3RwdWJrZXlfb3BzIiwiYWxnIjoiUlMyNTYifQ.eyJ0eXBlIjoiVXNlciIsInVpZCI6ImdjaGFwYXJybyIsImV4cCI6MTc4MzY4ODYyOCwiaWF0IjoxNzc4NTA0NjI4LCJpc3MiOiJodHRwczovL3Vycy5lYXJ0aGRhdGEubmFzYS5nb3YiLCJpZGVudGl0eV9wcm92aWRlciI6ImVkbF9vcHMiLCJhY3IiOiJlZGwiLCJhc3N1cmFuY2VfbGV2ZWwiOjN9.ClIy9LvUGX9Tev1vuF3_ppdWiG_EMmHjoYRBwPoyxgNbm9ie1OjyzaH6Qe7sKvu6TAitpVczLU0cgOMlQNnHgytqILDNYvNZFhaTchxbCz7qX-aXKTv3q_XYFHK6-x2I-OCg8xesYJsSsjJR_iU-hiO0_cnhTGGoT-E1HrjuSay8Jc67r_6DE-TivHIwo0vRPhq7Yh68ETIOFf7UFiymk4IEWiMfVUM0_uv7uFviet4oj7ay8if800MaspDdrhQMvCKma4ax4tM7DtMLKrRQy5VCscYiTqTUoCr3z6FdQfK6_92LShzQkWGFZZirnDwluTL-9aR6owp2vqSMh11THg'
import os
import requests
from datetime import datetime
from pathlib import Path
from dateutil.relativedelta import relativedelta
# Configuration
HEADERS = {"Authorization": f"Bearer {TOKEN}"}
PRODUCT = "MYD08_M3" # MODIS/Aqua Aerosol Cloud Water Vapor Ozone Monthly L3 Global 1Deg
COLLECTION = "61" # Versión más actual
# Define the range of months you want
START = datetime(2024, 1, 1)
END = datetime(2024, 12, 1)
foldername="modis_pwv_monthly"
OUTDIR = Path(foldername)
OUTDIR.mkdir(exist_ok=True)
def month_range(start, end):
current = start
while current <= end:
yield current
current += relativedelta(months=1)
for date in month_range(START, END):
year = date.strftime("%Y")
doy = date.strftime("%j") # For Monthly, this will be 001, 032, 060, etc.
# 1. Get the directory content list from the API
api_url = (
f"https://ladsweb.modaps.eosdis.nasa.gov/api/v2/"
f"content/details/allData/{COLLECTION}/"
f"{PRODUCT}/{year}/{doy}"
)
try:
r = requests.get(api_url, headers=HEADERS)
r.raise_for_status()
data = r.json()
# 2. Iterate through files in that month's folder
for item in data.get("content", []):
name = item["name"]
# Filter for the HDF file (ignoring .xml or .jpg)
if name.endswith(".hdf"):
file_url = item["downloadsLink"]
print(f"Downloading Monthly Data for {year}-Month {date.month}: {name}")
with requests.get(file_url, headers=HEADERS, stream=True) as rr:
rr.raise_for_status()
with open(OUTDIR / name, "wb") as f:
for chunk in rr.iter_content(chunk_size=16384):
f.write(chunk)
print(f"Successfully saved {name}")
except Exception as e:
print(f"Error processing {year}-{doy}: {e}")
fname=os.path.join(foldername, sorted(os.listdir(foldername))[0])
fname
# revisemos que el archivo sea del tipo correcto (para detectar errores de descarga)
# !file NOMBRE DE ARCHIVO
import os
print(os.path.getsize(fname) / 1024**2, "MB")
Revisemos los contenidos de un archivo.
from pyhdf.SD import SD, SDC
hdf = SD(fname, SDC.READ)
hdf.datasets()
sds_name = 'Atmospheric_Water_Vapor_Mean_Mean'
sds_obj = hdf.select(sds_name)
data = sds_obj.get().astype(float)
attrs = sds_obj.attributes()
attrs
scale_factor = attrs.get("scale_factor") # scale, here = 1e-3
add_offset = attrs.get("add_offset") # possible offset, here = 0
fill_value = attrs.get("_FillValue") # invalid data
data[data == fill_value] = np.nan # fill invalid data with nan
data = (data * scale_factor) + add_offset # reconstruct data
La siguiente es una función que extrae la información de PWV de un archivo y extrae la información latitud, longitud, pwv para hacer un mapa más adelante.
def process_hdf(file_path):
# Open HDF file
hdf = SD(file_path, SDC.READ)
sds_name = 'Atmospheric_Water_Vapor_Mean_Mean'
if sds_name not in hdf.datasets():
print(f"Error: SDS '{sds_name}' not found in {file_path}.")
print("Available datasets in this HDF file:")
for name in hdf.datasets():
print(f" - {name}")
raise ValueError(f"SDS '{sds_name}' not found in {file_path}. Please check the correct SDS name for MYD08_M3 products.")
sds_obj = hdf.select(sds_name)
data = sds_obj.get().astype(float)
attrs = sds_obj.attributes()
scale_factor = attrs.get("scale_factor")
add_offset = attrs.get("add_offset")
fill_value = attrs.get("_FillValue")
data[data == fill_value] = np.nan
data = (data * scale_factor) + add_offset # Standard MODIS scaling: (raw_value * scale_factor) + add_offset
# generate Lat/Lon grids (MODIS M3 is a global 1x1 degree grid)
# The grid is 180 rows (Lat) x 360 cols (Lon)
lats = np.linspace(90, -90, 180) # lat_start, lat_end, num_lat
lons = np.linspace(-180, 180, 360) # lon_start, lon_end, num_lon
lon_grid, lat_grid = np.meshgrid(lons, lats)
return lon_grid, lat_grid, data
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
# Colombia Geographic Extent [min_lon, max_lon, min_lat, max_lat]
COLOMBIA_EXTENT = [-82, -66, -5, 14]
lons, lats, pwv = process_hdf(fname)
fig = plt.figure(figsize=(10, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent(COLOMBIA_EXTENT, crs=ccrs.PlateCarree())
# Add Geographic Features
ax.add_feature(cfeature.COASTLINE, linewidth=1)
ax.add_feature(cfeature.BORDERS, linestyle=':', linewidth=1)
ax.add_feature(cfeature.STATES, linewidth=0.5, edgecolor='gray')
# Plot Data
mesh = ax.pcolormesh(lons, lats, pwv, cmap='Blues', transform=ccrs.PlateCarree(), vmin=0, vmax=6)
plt.colorbar(mesh, label='Precipitable Water Vapor (cm)', fraction=0.035, pad=0.04)
date_str = fname.split('.')[1][1:] # Extracts AYYYYDDD -> YYYYDDD
plt.title(f"Colombia PWV - {date_str}", fontsize=14)
plt.show()
Ejercicio: Amplíe el rango de cobertura para visualizar a México/Chile junto con Colombia.
Ejercicio: De forma programática, repita la gráfica para todos los meses de 2024 para Colombia.
Ejercicio: Encuentre el lugar sobre Colombia con el menor median-PWV del año y con el mínimo PWV durante el año 2024. Hacer una gráfica de PWV vs mes para el median-PWV y el mínimo para cada uno de esos lugares.
Avanzado: Repita esas gráficas para esos dos lugares, pero con promedio multianual usando todos los años disponibles de MODIS (~2002-2025), incluyendo barras de error.
Otras formas de acceder a datos satelitales#
La librería earthaccess permite hacer búsquedas con coincidencias espaciales a partir de algún rango geográfico. Requiere autenticarse interactivamente en la celda con las credenciales del portal Earthdata. En este corto ejemplo descargaremos datos de MODIS de nivel 2 (producto MOD21: MODIS/Terra Land Surface Temperature/3-Band Emissivity 5-Min L2 1km V061).
! pip install earthaccess
import earthaccess
earthaccess.login() # al ejecutar pide usuario/contraseña
results = earthaccess.search_data(
short_name="MOD21",
version="061",
temporal=(
"2016-01-01",
"2016-01-03"
),
bounding_box=(
-79.1,
-4.3,
-66.8,
13.5
)
)
print(f"Found {len(results)} granules.\n")
for r in results[:5]:
print(r)
earthaccess.download(
results,
local_path="mod21_colombia_week"
)
!pip install python-geotiepoints basemap basemap-data-hires
import geotiepoints as gp
from pyhdf import SD
import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import matplotlib as mpl
import os
%matplotlib inline
INPUT_DIR='mod21_colombia_week'
filename=sorted(os.listdir(INPUT_DIR))[2]
hdf=SD.SD(INPUT_DIR+'/'+filename)
# Read datasets
pwvds = hdf.select('PWV')
latitude = hdf.select('Latitude')[:]
longitude = hdf.select('Longitude')[:]
# Interpolate MODIS geolocation from 5 km to 1 km
lon1km, lat1km = gp.modis5kmto1km(longitude, latitude)
# Read PWV data and apply scaling
scale = pwvds.attributes()['_Scale']
z = pwvds[:] * scale
# Handle fill values
fill_value = pwvds.attributes().get('_FillValue')
if fill_value is not None and fill_value != 'n/a':
z = np.ma.masked_where(pwvds[:] == fill_value, z)
# Define discrete 1 cm levels
levels = np.arange(0, 11, 1)
# Discrete colormap + normalization
cmap = plt.cm.get_cmap('Blues', len(levels) - 1)
norm = mpl.colors.BoundaryNorm(
boundaries=levels,
ncolors=cmap.N,
clip=True
)
# Create figure
plt.figure(figsize=(10, 10))
m = Basemap(
projection='cyl',
resolution='l',
llcrnrlat=0,
urcrnrlat=13,
llcrnrlon=-78,
urcrnrlon=-70
)
# Coastlines and grid
m.drawcoastlines(linewidth=0.5)
m.drawparallels(
np.arange(-90., 120., 2.),
labels=[1, 0, 0, 0],
linewidth=0.5,
dashes=[1, 1]
)
m.drawmeridians(
np.arange(-180., 180., 2.),
labels=[0, 0, 0, 1],
linewidth=0.5,
dashes=[1, 1]
)
# Plot discrete PWV field
pcm = m.pcolormesh(
lon1km,
lat1km,
z,
latlon=True,
cmap=cmap,
norm=norm,
shading='nearest'
)
# Colorbar
from mpl_toolkits.axes_grid1 import make_axes_locatable
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", size="5%", pad=0.1)
cbar = plt.colorbar(
pcm,
cax=cax,
ticks=levels
)
cbar.set_label('Unit: ' + pwvds.attributes()['units'])
plt.title(filename)
plt.show()
Avanzado: Repita los ejercicios para datos de Optical Depth of Field, que se relaciona con presencia de material particulado (polución) en la atmósfera.