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.
Top comments (0)