Como llegar desde un dataset netCDF a una representación en RAWgraphs.
Tags: netcdf, r, datascience, datavisualization
Published: 2023-10-24 00:06:30


Hace unas semanas desde la ESA se lanzó el concurso “[little pictures](https://climate.esa.int/en/littlepicturescompetition/)” 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](https://github.com/littlepictures/datasets) pero quise aprovechar para para investigar un poco los datasets de “**[Climate Data Store](https://cds.climate.copernicus.eu/#!/home)”** 

Hay [montones de datasets de muchos tipos](https://cds.climate.copernicus.eu/cdsapp#!/search?type=dataset) 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:

![Representación NetCDF en panoply](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/alewwo23sl7o5mttyfb0.png)

Para abrirlo he usado [Panoply](https://www.giss.nasa.gov/tools/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](https://cds.climate.copernicus.eu/cdsapp#!/software/app-health-aedes-albopictus-projections?tab=app) 


## 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](https://cds.climate.copernicus.eu/cdsapp#!/home) 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://towardsdatascience.com/how-to-crack-open-netcdf-files-in-r-and-extract-data-as-time-series-24107b70dcd)
- [https://rpubs.com/boyerag/297592](https://rpubs.com/boyerag/297592)

Trabajando con el paquete [NETCDF](https://cran.r-project.org/web/packages/ncdf4/index.html)

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

```r
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:

```r
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

```r
#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 :

```r
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

```r
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)


![Diagrama de los datos temporales de un archivo netCDF](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/nhvf2fd4l6wjqtgkae1m.png)



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.

```r
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)
}
```


![Tabla de datos](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/w8kxcsn823fp1jkm1w8m.png)


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.


![Gráfica de evolución del dato en una ciudad determinada](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/gqjww5cc22m5v81ee3v6.png)




> 💡 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](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](https://climate.esa.int/en/little-pictures-gallery/) 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](https://es.wikipedia.org/wiki/Trayectorias_de_concentraci%C3%B3n_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](https://vdt-labs.com/easy-csv-editor/))  y añadir una columna nueva para diferenciar el modelo de emisiones utilizado ( rcp85 o rcp45)

```r
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](https://www.rawgraphs.io) Importar el archivo CSV y echar un vistazo


![Mapeo de datos en Raw Graphs](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/xz0kupgqsde0quv7245u.png)

![Gráficas de líneas en Raw Graphs](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/ukswmge702b31ca1eoe9.png)


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. 


![Animación mosquito tigre](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/bww7ah71auvsr6flu748.gif)

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)

![Exportar a SVG en la interface de raw graphs](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/78n9ucgc0s5fyasjowf2.png)


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.