Acceso a datos CALIPSO de lidar de nubes y aerosoles desde NASA Earthdata

Acceso a datos CALIPSO de lidar de nubes y aerosoles desde NASA Earthdata#

CAL_LID_L2_01kmCLay-Standard-V5-00 es el producto de datos Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) Lidar Level 2 1 km Cloud Layer. Este producto de datos fue recolectado usando el instrumento Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP). Dentro de este producto de capa de nubes, generado con una resolución horizontal de 1 km, hay dos clases generales de datos: Column Properties (incluyendo datos de posición y geometría de observación) y Layer Properties. Los productos de capas de nubes consisten en una secuencia de descriptores de columna, cada uno asociado con un número variable de descriptores de capas de nubes. Fuente: NASA Earthdata.

Requisitos#

Autenticación EDL (usuario/contraseña)
Obtener el archivo environment.yml e instalar el ambiente aislado de conda.
python>=3.12

Objetivos#

Descarga un archivo remoto con (sub)selection#

a) Por variables.
b) Por coordenadas geograficas.
c) Solo datos diurnos.

Descarga de múltiples archivos remotos#

A cada archivo se le aplica la subseleccion de una subregion de datos de interes antes de descargar.

Referencias#

Getzewich, B. (2025). CALIPSO Lidar Level 2 1 km Cloud Layer, V5-00 [Data set]. NASA Langley Atmospheric Science Data Center Distributed Active Archive Center. https://doi.org/10.5067/CALIOP/CALIPSO/CAL_LID_L2_01KMCLAY-STANDARD-V5-00

import xarray as xr
import datetime as dt
import earthaccess
import matplotlib.pyplot as plt
import numpy as np

# import pydap-specific tools
from pydap.client import get_cmr_urls, open_url, create_session
from pydap.client import to_netcdf as dap_to_netcdf

Como encontrar URLs de OPeNDAP#

Usando la API CMR de NASA#

Consultamos el Common Metadata Repository (CMR por sus siglas en ingles) de la NASA, para identificar archivos remotos que intersectan la siguiente área geográfica (caja delimitadora) y cubren el siguiente rango temporal:

-121 < longitud < -115, y 26.5 < latitud < 31
2 años de datos solo de primavera: del 1 de marzo al 31 de mayo (2022-2023).

Por último, nos interesan SOLO datos donde Daytime es True.

Calipso_L2_ccid = "C3463063995-LARC_CLOUD" # 
bbox = [-121,26.5,-115,31] # [west, south, east, north]

# 2 years of March data
time_ranges = [[dt.datetime(year, 3, 1), dt.datetime(year, 5, 31)] for year in range(2022, 2024)]

CMR_URLs = []
args = {
    "ccid": Calipso_L2_ccid,
    "bounding_box": bbox,
    "limit": 1000,
}
cmr_urls = [url for time_range in time_ranges for url in get_cmr_urls(**args, time_range=time_range)] # you can increase the limit of results

print("################################################ \n We found a total of ", len(cmr_urls), "OPeNDAP URLS!!!\n################################################")

################################################ 
 We found a total of  90 OPeNDAP URLS!!!
################################################

Autenticación EDL mediante earthaccess y OPeNDAP#

Puedes autenticarte mediante earthaccess como se muestra a continuacion. Debes tener una cuenta EDL válida. Hay dos estrategias para autenticarte con earthaccess:

strategy="interactive". Esto pedirá tu usuario y contraseña de EDL.
strategy="netrc". Usa esta opción si el notebook se ejecuta en un entorno donde se puede recuperar un .netrc con tus credenciales.

A continuación, el valor predeterminado será netrc, asumiendo que el usuario ejecutó el notebook Authenticate.ipynb. Si no es así, puedes cambiar la estrategia a "interactive".

from earthaccess.exceptions import LoginStrategyUnavailable
try:
    auth = earthaccess.login(strategy="netrc", persist=True) # you will be promted to add your EDL credentials
except LoginStrategyUnavailable:
    auth = earthaccess.login(strategy="interactive", persist=True)

# pass Token Authorization to a new Session.
my_session = create_session(session=auth.get_session())

Acceso SOLO a metadatos con PyDAP#

Podemos acceder a metadatos producidos por OPeNDAP para identificar las variables de interés. En particular, aquellas asociadas con valores de latitud y longitud.

%%time
pyds = open_url(cmr_urls[0], protocol="dap4", session=my_session)
pyds.tree()

.CAL_LID_L2_01kmCLay-Standard-V5-00.2022-03-02T10-34-17ZN.hdf
├──Lidar_Surface_Detection
│  ├──Surface_Top_Altitude_532
│  ├──Surface_Base_Altitude_532
│  ├──Surface_Integrated_Attenuated_Backscatter_532
│  ├──Surface_532_Integrated_Depolarization_Ratio
│  ├──Surface_532_Integrated_Attenuated_Color_Ratio
│  ├──Surface_Detection_Flags_532
│  ├──Surface_Overlying_Integrated_Attenuated_Backscatter_532
│  ├──Surface_Scaled_RMS_Background_532
│  ├──Surface_Peak_Signal_532
│  ├──Surface_Detections_333m_532
│  ├──Surface_Top_Altitude_1064
│  ├──Surface_Base_Altitude_1064
│  ├──Surface_Integrated_Attenuated_Backscatter_1064
│  ├──Surface_1064_Integrated_Depolarization_Ratio
│  ├──Surface_1064_Integrated_Attenuated_Color_Ratio
│  ├──Surface_Detection_Flags_1064
│  ├──Surface_Overlying_Integrated_Attenuated_Backscatter_1064
│  ├──Surface_Scaled_RMS_Background_1064
│  ├──Surface_Peak_Signal_1064
│  └──Surface_Detections_333m_1064
├──Ocean_Derived_Column_Optical_Depth
│  ├──ODCOD_Effective_Optical_Depth_532
│  ├──ODCOD_Effective_Optical_Depth_532_Uncertainty
│  ├──ODCOD_QC_Flag_532
│  ├──ODCOD_Surface_Wind_Speeds_10m
│  └──ODCOD_Surface_Wind_Speed_Correction
├──Lidar_Data_Altitudes
├──Profile_ID
├──Latitude
├──Longitude
├──Profile_Time
├──Profile_UTC_Time
├──Day_Night_Flag
├──Off_Nadir_Angle
├──Solar_Zenith_Angle
├──Solar_Azimuth_Angle
├──Scattering_Angle
├──Spacecraft_Position
├──Parallel_Column_Reflectance_532
├──Parallel_Column_Reflectance_Uncertainty_532
├──Perpendicular_Column_Reflectance_532
├──Perpendicular_Column_Reflectance_Uncertainty_532
├──Column_Integrated_Attenuated_Backscatter_532
├──Column_IAB_Cumulative_Probability
├──Column_Particulate_Optical_Depth_Above_Opaque_Water_Cloud_532
├──Column_Particulate_Optical_Depth_Above_Opaque_Water_Cloud_Uncertainty_532
├──Tropopause_Height
├──Tropopause_Temperature
├──IGBP_Surface_Type
├──Snow_Ice_Surface_Type
├──DEM_Surface_Elevation
├──Minimum_Laser_Energy_532
├──Low_Energy_Mitigation_Column_QC_Flag
├──Number_Layers_Found
├──Scene_Flag
├──Low_Energy_Mitigation_Feature_QC_Flag
├──Layer_Top_Altitude
├──Layer_Base_Altitude
├──Layer_Top_Pressure
├──Midlayer_Pressure
├──Layer_Base_Pressure
├──Layer_Top_Temperature
├──Layer_Centroid_Temperature
├──Midlayer_Temperature
├──Layer_Base_Temperature
├──Opacity_Flag
├──Attenuated_Scattering_Ratio_Statistics_532
├──Attenuated_Backscatter_Statistics_532
├──Integrated_Attenuated_Backscatter_532
├──Integrated_Attenuated_Backscatter_Uncertainty_532
├──Attenuated_Backscatter_Statistics_1064
├──Integrated_Attenuated_Backscatter_1064
├──Integrated_Attenuated_Backscatter_Uncertainty_1064
├──Volume_Depolarization_Ratio_Statistics
├──Integrated_Volume_Depolarization_Ratio
├──Integrated_Volume_Depolarization_Ratio_Uncertainty
├──Attenuated_Total_Color_Ratio_Statistics
├──Integrated_Attenuated_Total_Color_Ratio
├──Integrated_Attenuated_Total_Color_Ratio_Uncertainty
├──Overlying_Integrated_Attenuated_Backscatter_532
├──Layer_IAB_QA_Factor
├──Feature_Classification_Flags
├──CAD_Score
├──Initial_CAD_Score
├──metadata.Product_ID
├──metadata.Date_Time_at_Granule_Start
├──metadata.Date_Time_at_Granule_End
├──metadata.Date_Time_of_Production
├──metadata.Number_of_Good_Profiles
├──metadata.Number_of_Bad_Profiles
├──metadata.Initial_Subsatellite_Latitude
├──metadata.Initial_Subsatellite_Longitude
├──metadata.Final_Subsatellite_Latitude
├──metadata.Final_Subsatellite_Longitude
├──metadata.Orbit_Number_at_Granule_Start
├──metadata.Orbit_Number_at_Granule_End
├──metadata.Orbit_Number_Change_Time
├──metadata.Path_Number_at_Granule_Start
├──metadata.Path_Number_at_Granule_End
├──metadata.Path_Number_Change_Time
├──metadata.Lidar_L1_Production_Date_Time
├──metadata.Number_of_Single_Shot_Records_in_File
├──metadata.Number_of_Average_Records_in_File
├──metadata.Number_of_Features_Found
├──metadata.Number_of_Cloud_Features_Found
├──metadata.Number_of_Aerosol_Features_Found
├──metadata.Number_of_Indeterminate_Features_Found
├──metadata.Ocean_Fresnel_Reflection_Coefficient_532
├──metadata.MERRA2_Wind_Uncertainty
├──metadata.AMSR_Wind_Correction_Uncertainty
├──metadata.Lidar_Data_Altitudes
├──metadata.GEOS_Version
├──metadata.GMAO_Files_Used
├──metadata.Classifier_Coefficients_Version_Number
├──metadata.Classifier_Coefficients_Version_Date
└──metadata.Production_Script
CPU times: user 142 ms, sys: 4.89 ms, total: 147 ms
Wall time: 1.95 s

Descargar variables mínimas para identificar el subconjunto espacial y datos diurnos#

Las coordenadas son (nombres completamente calificados):

Latitude

Longitude

Day_Night_Flag

Antes de descargar, necesitamos identificar cualquier dimensión que también sea un arreglo del dataset.

(También puede haber dimensiones con nombre, es decir, dimensiones que solo tienen nombre y que NO están

asociadas con datos. No necesitamos declarar esas variables.)

DIMS = list(set(pyds['Latitude'].dims + pyds['Longitude'].dims + pyds['Day_Night_Flag'].dims))
dims = [dim for dim in DIMS if dim.split("/")[-1] in pyds[("/").join(DIMS[1].split('/')[:-1])].variables()] 
print("Dimensions that are also arrays: ", dims)

Dimensions that are also arrays:  []

output_path = "./data/"

Descarga de datos#

%%time
dap_to_netcdf(cmr_urls, session=my_session, 
              keep_variables = ["/Longitude", "/Day_Night_Flag"], 
              output_path=output_path)

CPU times: user 139 ms, sys: 224 ms, total: 363 ms
Wall time: 13.3 s

Inspeccionar todos los archivos descargados#

Aquí identificamos con más detalle el subconjunto de datos necesario en el archivo remoto, que descargará SOLO datos dentro de nuestra caja delimitadora, para cualquier variable de interés posible.

# Get data from Bounding Box
minLon, maxLon = bbox[0], bbox[2]

slices=[]
final_urls = []
for url in cmr_urls:
    filename = output_path+f"{url.split('/')[-1][:-4]}.nc4"
    dt1 = xr.open_datatree(filename).load()
    daytime_flag = dt1['Day_Night_Flag']
    # find index /data_01/longitude
    longitude = dt1['/Longitude']
    mask = (longitude >= minLon) & (longitude <= maxLon)
    idx = np.nonzero(mask.values)[0]
    daytime_flag = dt1['Day_Night_Flag'].isel(Record_Number=slice(idx[0], idx[-1]))==1
    if all(daytime_flag==0):
        final_urls.append(url)
        slices.append({"/Record_Number":(idx[0], idx[-1])})

print(f"\nOnly {len(final_urls)} out of the {len(cmr_urls)} remote files satisfy our Daylight Criteria\n")
print("Sample subsetting slices:")
slices[:4]

Only 42 out of the 90 remote files satisfy our Daylight Criteria

Sample subsetting slices:

[{'/Record_Number': (np.int64(11471), np.int64(14314))},
 {'/Record_Number': (np.int64(14165), np.int64(16225))},
 {'/Record_Number': (np.int64(12422), np.int64(14979))},
 {'/Record_Number': (np.int64(10329), np.int64(13187))}]

Inspeccionar visualmente el subconjunto#

# Subset data
Lon = dt1['Longitude'].isel(Record_Number=slice(idx[0], idx[-1]))

# Generate masked data to visualize only

Lon_masked = xr.full_like(dt1['Longitude'], np.nan)
Lon_masked.loc[dict(
    Record_Number = Lon['Record_Number'] + idx[0]
)] = Lon

# Visualize: Plot subset of data over original data
fig, axes = plt.subplots(figsize=(10,4))
dt1['Longitude'].plot(lw=5, color='k', alpha=0.75);
Lon_masked.plot(lw=10, color="#7f00ff")
axes.set_title(r"Longitude Subset $[^\circ$E]")
plt.show()

../_images/c93df6d12cf2e09378eb358eab5139b23254e35fe19147904af5e75b37f322e5.png

Descargar todos los datos de interés#

PRIMERO: necesitamos borrar todos los archivos descargados previamente para evitar colisiones de nombres.

import os
import glob

fnames = [output_path+f"{fname.split('/')[-1][:-4]}.nc4" for fname in cmr_urls]
for filename in fnames:
    try:
        os.remove(filename)
    except FileNotFoundError:
        print(f"The file '{filename}' is not in there anymore")    

# Will Download 34 Variables!
keep_variables = [
    '/Lidar_Surface_Detection', # <----- ALL Variables inside Group
    "/Ocean_Derived_Column_Optical_Depth", # < -- All varibles inside Group
    "/Lidar_Data_Altitudes", "/Profile_ID", "/Latitude", "/Longitude", 
    "/Profile_Time", "/Profile_UTC_Time", "/Day_Night_Flag", "/Tropopause_Height", 
    "/Tropopause_Temperature",
]

%%time
dap_to_netcdf(final_urls, session=my_session,
              keep_variables = keep_variables,
              dim_slices = slices,
              output_path=output_path)

CPU times: user 94.5 ms, sys: 291 ms, total: 386 ms
Wall time: 51.6 s