Tutorial de 5 minutos

Tutorial de 5 minutos#

OPeNDAP: la visión#

La visión original de OPeNDAP (Cornillion, et al 1993) es la de democratizar el acceso remoto a datos, estableciendo las equivalencias

\( \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \boxed{\text{URL} \approx \text{Conjunto de datos remoto} }\)

\( \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \boxed{\text{URL + Restricciones} \approx \text{Subconjunto del conjunto de datos remoto}} \)

Eso llevó al desarrollo del protocolo DAP2 (antes conocido como DODS). Actualmente, los servidores OPeNDAP y Unidata implementan el protocolo moderno y más amplio DAP4 (consulta la especificación DAP4), para seguir habilitando la visión original de OPeNDAP.

Pydap#

La lógica interna de PyDAP permite construir expresiones de restricción para cada URL de forma interactiva, de tal manera que el usuario ya no debe de hacerlo manualmente. Esto permite el acceso a cientos de archivos remotos de manera interactiva y eficiente. Aun así, una comprensión básica del uso de expresiones de restricción (CEs por su siglas en ingles) puede resultar muy útil al trabajar con múltiples archivos remotos, y conducir a flujos de trabajo más eficientes.

Objetivos:#

  • Demostrar cómo especificar el protocolo DAP4 al servidor remoto mediante PyDAP.

  • Demostrar el uso de Xarray con PyDAP para la exploracion de datos remotos de la NASA.

  • Demostrar distintas formas de usar expresiones de restricción (CEs), y cómo asegurar que el servidor remoto es el que implementa el segmento de subconjunto (slice en ingles), de forma próxima a los datos, sin pérdida de rendimiento del lado del usuario.

Requisitos#

  • Archivos de datos cientificos detrás de un servidor de OPeNDAP que implemente el protocolo DAP4. En este caso se usara el servidor de pruebas: http://test.opendap.org/opendap/.

  • pydap

  • xarray

  • python>=3.12

Note

La gran mayoría de los servidores OPeNDAP de NASA implementan el protocolo DAP4.

from pydap.client import open_url, consolidate_metadata, create_session
import xarray as xr
import numpy as np

# create a session to inspect downloads. cache_name must have `debug`
session = create_session(use_cache=True, cache_kwargs={"cache_name":'data/debug_case1'})
session.cache.clear()

Segmento de subconjunto de dos archivos distintos.#

Considera el caso en el que el servidor remoto permite el acceso a dos archivos distintos que forman parte de una misma simulacion, or proyecto. En este case usaremos coads_climatology y coads_climatology2. Estos dos archivos de datos comparten dimensiones espaciales idénticas, se pueden agregar en el tiempo, y comparten casi todas las mismas variables.

Note

Es importante siempre comprobar que los datos the dos distintos archivos se puedan agregar. PyDAP y Xarray tienen lógica interna para verificar si dos o más conjuntos de datos se pueden aggregar.

Un paso importante será usar Expresiones de Restricción (CEs) para asegurar que solo se concatenen las variables de interés.

Warning

Uno de estos archivos tiene variables extra que no están presentes en el otro archivo y que serán descartadas mediante el uso de CEs.

urls = ["http://test.opendap.org/opendap/data/nc/coads_climatology.nc", "http://test.opendap.org/opendap/data/nc/coads_climatology2.nc"]
dap4_urls = [url.replace("http","dap4") for url in urls]

# constraint expression
dap4_CE = "?dap4.ce=" + ";".join(["/SST", "/COADSX", "/COADSY", "/TIME"])

# Final list of OPeNDAP URLs
dap4ce_urls =[url+dap4_CE for url in dap4_urls]
print("====================================================\nThe following are the DAP4 OPeNDAP URLs \n", dap4ce_urls)

====================================================
The following are the DAP4 OPeNDAP URLs 
 ['dap4://test.opendap.org/opendap/data/nc/coads_climatology.nc?dap4.ce=/SST;/COADSX;/COADSY;/TIME', 'dap4://test.opendap.org/opendap/data/nc/coads_climatology2.nc?dap4.ce=/SST;/COADSX;/COADSY;/TIME']

Note

¿Por qué usar CEs cuando Xarray tiene un método .drop_variables? Porque Xarray primero necesita analizar TODOS los metadatos asociados con cada archivo remoto para después descartar las variables. En algunos archivos remotos puede haber miles de variables. En ese caso Xarray analizaría todas las variables, para luego descartar las que no son de interes. Con la CE, el servidor solo envía metadatos restringidos asociados solamente con las variables deseadas.

Warning

Xarray espera la presencia de dimensiones en los metadatos. Al construir CEs, el usuario debe asegurarse de incluir en la CE todas las dimensiones asociadas con las variables de interés. En el ejemplo anterior, COASX, COADSY y TIME son las dimensiones de SST.

Consolidate Metadata acelera la generación del Dataset.#

consolidate_metadata(dap4ce_urls, session=session, concat_dim="TIME")

datacube has dimensions ['COADSX[0:1:179]', 'COADSY[0:1:89]'] , and concat dim: `['TIME']`

Note

consolidate_metadata(dap4_urls, concat_dim='...', session=session) descarga las dimensiones del archivo remoto y las almacena como SQLite para reutilizarlas. Ademas de actuar como gestor de base de datos, la sesión permite autenticar. Esta práctica puede producir una mejora de rendimiento de flujos de trabajo aproximadamente 10 a 100 veces más rápidos.

Usar la lógica de Xarray para acceso a metadatos.#

ds = xr.open_mfdataset(
    dap4ce_urls, 
    engine='pydap',
    concat_dim='TIME',
    session=session,
    combine="nested",
    parallel=True,
    decode_times=False,
    chunks={},
)
ds

<xarray.Dataset> Size: 2MB
Dimensions:  (TIME: 24, COADSY: 90, COADSX: 180)
Coordinates:
  * TIME     (TIME) float64 192B 366.0 1.096e+03 ... 7.671e+03 8.401e+03
  * COADSY   (COADSY) float64 720B -89.0 -87.0 -85.0 -83.0 ... 85.0 87.0 89.0
  * COADSX   (COADSX) float64 1kB 21.0 23.0 25.0 27.0 ... 375.0 377.0 379.0
Data variables:
    SST      (TIME, COADSY, COADSX) float32 2MB dask.array<chunksize=(12, 90, 180), meta=np.ndarray>
Attributes:
    history:  FERRET V4.30 (debug/no GUI) 15-Aug-96
ds['SST']

<xarray.DataArray 'SST' (TIME: 24, COADSY: 90, COADSX: 180)> Size: 2MB
dask.array<concatenate, shape=(24, 90, 180), dtype=float32, chunksize=(12, 90, 180), chunktype=numpy.ndarray>
Coordinates:
  * TIME     (TIME) float64 192B 366.0 1.096e+03 ... 7.671e+03 8.401e+03
  * COADSY   (COADSY) float64 720B -89.0 -87.0 -85.0 -83.0 ... 85.0 87.0 89.0
  * COADSX   (COADSX) float64 1kB 21.0 23.0 25.0 27.0 ... 375.0 377.0 379.0
Attributes:
    long_name:  SEA SURFACE TEMPERATURE
    history:    From coads_climatology
    units:      Deg C
    Maps:       ('/TIME', '/COADSY', '/COADSX')

Note

El chunking de SST implica que toda la variable SST de cada arvhido es un solo chunk. Esta es una interpretación estándar que Xarray hace de las URLs de OPeNDAP. ¿Qué ocurre si descargamos un solo punto espacial desde un único archivo remoto?

session.cache.clear()

%%time
ds['SST'].isel(TIME=0, COADSX=0, COADSY=0).load() # this should download a single point one of the files

CPU times: user 22.6 ms, sys: 11 ms, total: 33.6 ms
Wall time: 478 ms

<xarray.DataArray 'SST' ()> Size: 4B
array(-1.e+34, dtype=float32)
Coordinates:
    COADSX   float64 8B 21.0
    COADSY   float64 8B -89.0
    TIME     float64 8B 366.0
Attributes:
    long_name:  SEA SURFACE TEMPERATURE
    history:    From coads_climatology
    units:      Deg C
    Maps:       ('/TIME', '/COADSY', '/COADSX')
print("====================================== \n Request sent to the Remote Server:\n ", session.cache.urls()[0].split("?")[-1].split("&dap4.checksum")[0].replace("%5B","[").replace("%5D","]").replace("%3A",":").replace("%2F","/"), "\n====================================== ")

====================================== 
 Request sent to the Remote Server:
  dap4.ce=/SST[0:1:11][0:1:89][0:1:179] 
======================================

Toda la variable completa se descarga innecesariamente !!#

Esto es muy ineficiente. Lo ideal seria que Xarray solo solicite (request en ingles) mediante una expresión de restricción el subconjunto de datos de interest. En este caso, lo ideal hubiera sido lo siguiente:

dap4.ce=/SST[0:1:0][0:1:0][0:1:0]

Apesar de que Pydap puede generar las expressiones de restriccion, cuando se utiliza xr.open_mfdataset este no instruye a pydap el generar la condicion de restriccion relevante. En su lugar, el metodo arriba descarga toda la variable (todo el chunk) y ya estando la variable en memoria, la subdivide y descarta lo que no es de interes. Este proceso no es muy eficiente y conlleva a descargar innecesariamente variables completas.

Cómo pasar el la segmentacion del archivo al servidor remoto cuando se usa Xarray?#

La respuesta es definir chunks al momento de crear el Xarray.Dataset. chunks debera declarar, por cada variable, el tamaño esperado de tu subconjunto.

Warning

Si fragmentas el dataset con un tamaño menor que tu descarga esperada, activarás muchas descargas por archivo remoto, obligando a Xarray a hacer trabajo extra para ensamblar los datos.

A continuacion se demuestra como lograr ese objectivo

# consolidate metadata again, since the cached metadata was cleared before
session = create_session(use_cache=True, cache_kwargs={"cache_name":'data/debug_case2'})
session.cache.clear()
consolidate_metadata(dap4ce_urls, session=session, concat_dim="TIME")

datacube has dimensions ['COADSX[0:1:179]', 'COADSY[0:1:89]'] , and concat dim: `['TIME']`

# For a single element in all dimensions, the expected size of the download is:
expected_sizes = {"TIME":1, "COADSX":1, "COADSY":1}

%%time
ds = xr.open_mfdataset(
    dap4ce_urls, 
    engine='pydap',
    concat_dim='TIME',
    session=session,
    combine="nested",
    parallel=True,
    decode_times=False,
    chunks=expected_sizes,
)
session.cache.clear()

CPU times: user 248 ms, sys: 17.6 ms, total: 265 ms
Wall time: 311 ms

ds['SST'] # inspect chunks before download

<xarray.DataArray 'SST' (TIME: 24, COADSY: 90, COADSX: 180)> Size: 2MB
dask.array<concatenate, shape=(24, 90, 180), dtype=float32, chunksize=(1, 1, 1), chunktype=numpy.ndarray>
Coordinates:
  * TIME     (TIME) float64 192B 366.0 1.096e+03 ... 7.671e+03 8.401e+03
  * COADSY   (COADSY) float64 720B -89.0 -87.0 -85.0 -83.0 ... 85.0 87.0 89.0
  * COADSX   (COADSX) float64 1kB 21.0 23.0 25.0 27.0 ... 375.0 377.0 379.0
Attributes:
    long_name:  SEA SURFACE TEMPERATURE
    history:    From coads_climatology
    units:      Deg C
    Maps:       ('/TIME', '/COADSY', '/COADSX')
%%time
ds['SST'].isel(TIME=0, COADSX=0, COADSY=0).load() # triggers download of an individual chunk

CPU times: user 144 ms, sys: 14.1 ms, total: 158 ms
Wall time: 398 ms

<xarray.DataArray 'SST' ()> Size: 4B
array(-1.e+34, dtype=float32)
Coordinates:
    COADSX   float64 8B 21.0
    COADSY   float64 8B -89.0
    TIME     float64 8B 366.0
Attributes:
    long_name:  SEA SURFACE TEMPERATURE
    history:    From coads_climatology
    units:      Deg C
    Maps:       ('/TIME', '/COADSY', '/COADSX')
print("====================================== \n Request sent to the Remote Server:\n ", session.cache.urls()[0].split("?")[-1].split("&dap4.checksum")[0].replace("%5B","[").replace("%5D","]").replace("%3A",":").replace("%2F","/"), "\n====================================== ")

====================================== 
 Request sent to the Remote Server:
  dap4.ce=/SST[0:1:0][0:1:0][0:1:0] 
======================================

Exito? Una advertencia cuando se trabaja con chunks#

En el caso de arriba, si se logro descargar exactamente lo que solicitamos. Sin embargo, en algunos escenarios, el tiempo de descarga puede ser igual o mas más lento en comparación con el caso en que se solicita toda la variable. La razón puede atribuirse al número de chunks que genera Dask al orquestrar el proceso.

  • Sin chunks. Descarga todo el arreglo del archivo. 2 chunks en 5 graficas dask (uno por archivo).`

  • Con chunks. Descargar solo el elemento deseado de un archivo. 388800 chunks en 5 graficas de dask`.

Idealmente, Dask, el gestor de chunks, solo debería activar la descarga de un único chunk. Sin embargo, se crearon 388800 chunks virtuales para asegurar que la seleccion del subconjunto se pasara al servidor. Esto a veces puede causar que el proceso tarde mas para el usuario.

En el escenario anterior fuimos al extremo de definir chunks con el elemento mas pequeno posible. A continuacion demostraremos algo mas eficiente para este caso, subdividiendo en la dimension a lo largo the la dimension COADSY (en ambos archivos).

session = create_session(use_cache=True, cache_kwargs={"cache_name":'data/debug_case3'})
session.cache.clear()

consolidate_metadata(dap4ce_urls, session=session, concat_dim="TIME")

datacube has dimensions ['COADSX[0:1:179]', 'COADSY[0:1:89]'] , and concat dim: `['TIME']`

download_sizes = {"COADSY":1} # note that we will subset across all time

%%time
ds = xr.open_mfdataset(
    dap4ce_urls, 
    engine='pydap',
    concat_dim='TIME',
    session=session,
    combine="nested",
    parallel=True,
    decode_times=False,
    chunks=download_sizes,
)
session.cache.clear()

CPU times: user 46.3 ms, sys: 9.3 ms, total: 55.6 ms
Wall time: 94.7 ms

ds['SST']

<xarray.DataArray 'SST' (TIME: 24, COADSY: 90, COADSX: 180)> Size: 2MB
dask.array<concatenate, shape=(24, 90, 180), dtype=float32, chunksize=(12, 1, 180), chunktype=numpy.ndarray>
Coordinates:
  * TIME     (TIME) float64 192B 366.0 1.096e+03 ... 7.671e+03 8.401e+03
  * COADSY   (COADSY) float64 720B -89.0 -87.0 -85.0 -83.0 ... 85.0 87.0 89.0
  * COADSX   (COADSX) float64 1kB 21.0 23.0 25.0 27.0 ... 375.0 377.0 379.0
Attributes:
    long_name:  SEA SURFACE TEMPERATURE
    history:    From coads_climatology
    units:      Deg C
    Maps:       ('/TIME', '/COADSY', '/COADSX')
%%time
ds['SST'].isel(COADSX=0, COADSY=0).load()

CPU times: user 24.2 ms, sys: 7.62 ms, total: 31.8 ms
Wall time: 352 ms

<xarray.DataArray 'SST' (TIME: 24)> Size: 96B
array([-1.e+34, -1.e+34, -1.e+34, -1.e+34, -1.e+34, -1.e+34, -1.e+34,
       -1.e+34, -1.e+34, -1.e+34, -1.e+34, -1.e+34, -1.e+34, -1.e+34,
       -1.e+34, -1.e+34, -1.e+34, -1.e+34, -1.e+34, -1.e+34, -1.e+34,
       -1.e+34, -1.e+34, -1.e+34], dtype=float32)
Coordinates:
  * TIME     (TIME) float64 192B 366.0 1.096e+03 ... 7.671e+03 8.401e+03
    COADSX   float64 8B 21.0
    COADSY   float64 8B -89.0
Attributes:
    long_name:  SEA SURFACE TEMPERATURE
    history:    From coads_climatology
    units:      Deg C
    Maps:       ('/TIME', '/COADSY', '/COADSX')
print("====================================== \n Parallel Requests sent to the Remote Server:\n ", [url.split("?")[-1].split("&dap4.checksum")[0].replace("%5B","[").replace("%5D","]").replace("%3A",":").replace("%2F","/") for url in session.cache.urls()], "\n====================================== ")

====================================== 
 Parallel Requests sent to the Remote Server:
  ['dap4.ce=/SST[0:1:11][0:1:0][0:1:179]', 'dap4.ce=/SST[0:1:11][0:1:0][0:1:179]'] 
======================================

¡Éxito! Una descarga mas rapida y mucho más pequeña que la original!#