Como llegar desde un dataset netCDF a una representación en RAWgraphs.

Hace unas semanas desde la ESA se lanzó el concurso “little pictures” Se pide a los participantes transformar datos derivados de observaciones satelitales a pequeñas visualizaciones que llamen la atención sobre el cambio climático.

Ofrecen algunos datasets ya preparados pero quise aprovechar para para investigar un poco los datasets de “Climate Data Store”

Hay montones de datasets de muchos tipos Vienen en un formato llamado netCDF que nunca habíamos trabajado, ¡otro motivo más para dedicar un rato a “little pictures”! Sin ser muy experto ni en R ni en netCDF, voy a tratar compartir lo que he aprendido.

Lo primero que llama la atención de los archivos netCDF es que pesan bastante. No es un excel de unos pocos megas sino que es fácil que ocupen 100 mb por archivo. ¿Por qué?

Vamos a ver un ejemplo de una representación de uno de estos archivos:

Para abrirlo he usado Panoply que es fácil aunque muy básico. QGIS también lo abre.
Observamos una gráfica con latitud-longitud en el que con cierta resolución (décima de grado, centésima de grado…) se representa un valor.

Por tanto vamos a tener ya un montón de puntos ahí. Pero es que además, lo habitual es tener no solo una “foto fija” sino también el valor de la evolución a lo largo de un periodo de tiempo. Y no solo de una variable sino incluso de varias. Ya vamos viendo de donde salen los 100 megas del archivo.

Estos archivos además son “autodescriptivos” es decir llevan incorporada la metadata sobre su contenido.

Las captura y el resto del artículo los voy a hacer con este dataset Sobre la duración de la actividad del mosquito tigre en europa

Simplificando

Si hacemos visualizaciones web no nos gusta obligar a la gente a descargar 100 megas de datos. Sobre todo si lo que nos interesa es responder a una pregunta simple, por ejemplo: ¿Cual va a ser la evolución del mosquito tigre en una ciudad?

En este caso voy a querer pasar de este formato de datos donde tengo un grid de puntos y un valor (llamemosle un formato 3D) + el tiempo (4D)A un formato 2D donde solo quiero la evolución temporal del valor dado un lugar en el mapa.

Código

Desde el proyecto de Climate Data Store ofrecen muchos ejemplos de código en python y aplicaciones para copiar. Pero para estas cosas me manejo un poco mejor en R. Encontré estos recursos que os recomiendo leer para entender mejor lo de los archivos netCDF

• https://towardsdatascience.com/how-to-crack-open-netcdf-files-in-r-and-extract-data-as-time-series-24107b70dcd
• https://rpubs.com/boyerag/297592

Trabajando con el paquete NETCDF

Antes de todo vamos a definir las ciudades que nos interesan, por ejemplo estas.

POIS <- list(
  'Berlin'= c(13.404954, 52.520008), #long lat
  'Budapest'= c(19.040236, 47.497913),
  'Paris'= c(2.354478, 48.860506)
)

Como decía arriba, el propio archivo contiene su metadata. Así que en este caso si lo abrimos nos dice su contenido:

nc_file <- nc_open("mosquito_seas_rcp85_mean_v1.0.nc") # change 85 for 45 for other experiment
print(nc_file)

# saca entre otras cosas:
2 variables (excluding dimension variables):        
    float season_length[lon,lat,time]   (Contiguous storage)              
    _FillValue: NaN            
    units: 1            
    long_name: Ensemble members average of VBD season_length  for future climate under rcp85            
    coordinates: height
....

3 dimensions:        
lat  Size:425             
        _FillValue: NaN            
        units: degrees_north            
        standard_name: latitude        
lon  Size:599             
        _FillValue: NaN            
      units: degrees_east            
        standard_name: longitude        
time  Size:100             
        standard_name: time            
        units: days since 1986-01-01 00:00:00            
        calendar: proleptic_gregorian

Donde vemos que nos dice las dimensiones con sus unidades. Ojo que aquí no hay estándar. Por ejemplo, las unidades de tiempo no tienen por qué ser UNIX time, aquí dice que son días transcurridos desde 1 de Enero de 1986.

En los artículos de arriba vienen mejor explicadas las formas de acceder al contenido del archivo usando ncvar_get y ncatt_get

#ncvar_get nos devuelve la estructura de datos de una de las variables
time <- ncvar_get(nc_file,"time") # days after first datum
##
> head(time)
[1]    0  365  730 1096 1461 1826
> length(time)
[1] 100

Vemos que tenemos el array de tiempo con 100 valores (100 años) desde 0 y sumando 365 (o 366) días

si hacemos :

seasonLength_array <- ncvar_get(nc_file,"season_length")

Ahí vamos a tener un montón de datos. ¡¡TODOS los valores de season_length!! Así que vamos a ir simplificando. Podemos extraer un fotograma fijo de esos datos eligiendo un momento en el tiempo. Por ejemplo el año 1 o el año 20

seasonLength_slice <- seasonLength_array[,,20]
image(lon,lat, seasonLength_slice, col=brewer.pal(8,"YlOrRd"))

Al fijar un año lo que tenemos es el conjunto de todos los valores para ese año en función de (Lat, Lon)

Pero lo que queremos hacer no es fijar el año sino la longitud y latitud. Para ello necesitamos encontrar el indice de las coordenadas más cercanas a las coordenadas de cada ciudad.

seasonLength_df <- data.frame(city=character(), year=integer(), value=numeric())
seasonLength_byTime <- array()

for (city in names(POIS)) {
  # get the lat and lon for the current city
  sample_lat <- POIS[[city]][2]
  sample_lon <- POIS[[city]][1]

  # find the closest lat/lon index in the grid
  lon_index <- which.min(abs(lon - sample_lon))
  lat_index <- which.min(abs(lat - sample_lat))

  # get the seasonLength_byTime for the current city
  seasonLength_byTime <- seasonLength_array[lon_index, lat_index, ]

  #Adjust the data to the dataframe format
  seasonLength_df_city <- data.frame(city=rep(city, length(time_obs)), year=format(time_obs, "%Y"), value=seasonLength_byTime)
  seasonLength_df <- rbind(seasonLength_df, seasonLength_df_city)
}

El dataframe de resultados está ya listo para exportarse a CSV, y tenemos las tres ciudades con su nombre en la primera columna.

Si representamos el valor de seasonLength_byTime en cada iteración, tenemos la evolución en una de las ciudades.

💡 Por último indicar que esto es una forma un poco de andar por casa, seguramente sea más potente usar el package “Raster” https://pjbartlein.github.io/REarthSysSci/raster_intro.html Uno de los problemas que se me ocurre que pueden aparecer es que el dataset sea muy variable en la dimensión espacial. Es decir que una centésima de grado mas allá o mas acá oscile mucho el valor. En este caso habría que calcular una media de una región, algo que costaría más hacer con este método

Siguientes pasos

Para hacer un cartel estilo los de little pictures quedaría representar los datos con un objetivo comunicativo. Por ejemplo podríamos mirar a ver que pasa al comparar el escenario RCP 4.5 y 8.5 y el año 0 ( 1985)

💡 Cuando los científicos pasan de “observaciones” a modelos de predicción utilizan distintos escenarios futuros llamados “trayectoria de concentración representativa” https://es.wikipedia.org/wiki/Trayectorias_de_concentración_representativas que se diferencian en cuantas emisiones de CO2 vamos a ser capaces de recortar en este siglo. Este dataset se ofrece para trayectorias 4.5 y 8.5

Podemos generar dos CSVs cambiando el archivo de entrada y luego pegarlos a mano usando una hoja de cálculo (por aquí somos muy fans del easy CSV editor) y añadir una columna nueva para diferenciar el modelo de emisiones utilizado ( rcp85 o rcp45)

experiment  city    year    value
rcp85   Berlin  1986    79.9151763916016
rcp85   Berlin  1987    79.9485092163086
rcp85   Berlin  1988    81.0975112915039
rcp85   Berlin  1989    80.5817642211914

Luego podemos ir RAWgraphs Importar el archivo CSV y echar un vistazo

Si lo comparamos con el mapa se aprecia como la duración de la estación del mosquito crece, sobre todo en países hacia el este, donde ahora es muy corta.

Tal y como lo veo, RAWgraphs no está diseñado para obtener una gráfica ya perfecta, sino más bien como una herramienta intermedia para generar visualizaciones básicas que puedan exportarse a SVG (y desde ahí trabajarlas en otro software, tipo Figma o Illustrator)

A partir de aquí, ahora llegaría el momento de hacer algo bonito con estas gráficas: darle una vuelta, pensar de qué manera los datos crean mas impacto, comunican mejor... o lo que sea que queramos conseguir. Ahí ya depende de cada uno.