1 Zusammengefasst

In dieser Lektion geht es um das praktische Arbeiten mit dem Datenformat des Arrays, das sowohl in R als auch in Python dazu dient, multidimensionale Daten zu speichern.

2 Legende

Wie alle anderen Lektionen benutzen wir auch hier die folgende Notation:

R-Code (den ihr über copy & paste in Euer Script überführen und ausführen könnt) wird so dargestellt:

# Hashtag kennzeichnet einen Kommentar; wird von R ignoriert = nicht ausgeführt

Die Ausgabe von Befehlen in R wird hier so dargestellt:

## [1] "Die Ausgabe eines R-Befehls"

Aufgaben zu den einzelnen Kapiteln werden eingerückt dargestellt:

Aufgaben

Lösungen, Zusatzaufgaben sowie sowie wichtige Hinweise und Tipps haben folgende Buttons:

3 Multidimensionale vs. multivariate Daten

In den Erd- und Umweltwissenschaften verlassen wir uns oft auf Beobachtungen und Messungen. Wenn es sich dabei um eine Meßgröße wie z.B. Temperatur, Niederschlag, Bodentyp oder Bestandsdichte handelt, die sich räumlich oder zeitlich ändert, sprechen wir von multdimensionalen Daten. Hierbei ist die Definition von Dimension meist physikalisch motiviert (also Raum und Zeit). Sobald unsere Beobachtungen mehr als nur eine Eigenschaft systematisch aufnehmen, sprechen wir üblicherweise von multivariaten Daten. Letztlich können Datensätze auch sowohl multidimensional wie multivariat sein, wenn wir uns z.B. für das gleichzeitige Ausbreitungsmuster mehrerer invasiver Arten interessieren.

Digitale Videos wie z.B. Regenradarfilme sind multidimensionale Datenformate, und bieten ihre Information in einer rechteckigen Anordnung von Pixelwerten (üblicherweise rot, grün and blau, auch RGB) an. Während das rechteckige Format des Videos unverändert bleibt, ändern sich ortsfeste RGB-Pixelwerte mit der Zeit. Hingegen sind topographische Karten, die Punkte, Linien, Polygone und Textzeichen in einem kompakten Format kombinieren, multivariate Daten. Geschriebenen Text kann man auch als multivariate Daten auffassen. Wir können die Anzahl einzelner Wörter bestimmen, testen ob sie korrekte Rechtschreibung haben, die Sprache festlegen, usw. Social media Konten beherbergen sowohl multidimensionale als auch multivariate Daten, die oft lückenhaft sein können: Benutzer*innen möchten vielleicht nur begrenzt multivariate Informationen (Alter, Hobbies, etc.) in ihren Profilen bereitstellen, wohingegen z.B. ihre Likes etwa über Ortungsdienste in Zeit und Raum, also multidimensional, geloggt werden.

Wissenswertes über multidimensionale (und multivariate) Daten

  • Die Sichtung und Auswahl bestimmter Informationen (= Variablen) ist besonders wichtig, aber auch herausfordernd; manchmal sind alle Daten in einem Datensatz interessant, oft ist es aber eben nur ein Teil davon.

  • Die statistische Interpretation von multidimensionalen Daten ist etwas schwieriger und mit mehr Fallstricken versehen, verglichen mit Datensätzen, die nur ein oder zwei Variablen haben.

  • Multidimensionale Daten kombinieren oft unterschiedliche Typen oder Messeinheiten. Oft müssen wir die Daten transformieren, um eine sinnvolle Vergleichsbasis zu haben.

  • Multimensionale Daten können auf mehrere Arten gespeichert werden; einige dieser Arten sind für ein bestimmtes Problem vorteilhafter als andere. Es lohnt sich, eine passende Datenstruktur zu wählen.

  • Die Dokumentation vom multidimensionalen Daten (auch bekannt als Metadaten) ist essenziell, damit Nutzer*innen alle Details über die Datenbeschaffung, -qualität, -auflösung, usw. nachvollziehen können.

  • Datenlücken brauchen spezielle Aufmerksamkeit. Oft überwiegen in Datensätzen die fehlenden Einträge.

  • Erfreulicherweise gibt es sowohl in R als auch in Python viele Algorithmen und Bibliotheken, die auf das Arbeiten mit mutltidimensionalen Daten spezialisiert sind.

4 Vorbereitung

Wir beginnen mit den üblichen Schritten, d.h. wir leeren die Arbeitsumgebung und definieren das Arbeitsverzeichnis:

# Säubern
graphics.off()
rm(list = ls(all = TRUE))

# Arbeitsverzeichnis setzen
setwd("~/Desktop/FUND/MultiDim")

5 Datentypen und Datenformate

Wir erinnern uns daran, dass wir Datentypen von Datenformaten unterscheiden können. Dazu dienen z.B. die Befehle typeof() und class(). Die gängigsten Datenformate in R sind Vektoren, Matrizen und Dataframes. Letztere eignen sich besonders für die Speicherung multidimensionaler Daten, v.a. weil wir für jede Spalte einen unterschiedlichen Datentyp verwenden können. Dennoch bleibt die Struktur von Dataframes nur zweidimensional, d.h. es gibt nur Reihen und Spalten.

Spannend wird es also, neue Dimensionen hinzuzufügen. Hierzu gibt es in R das Format des array, das wir mit dem gleichnamigen Befehl erstellen können:

(bsp <- array(1:10, c(2, 3, 4)))
## , , 1
## 
##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6
## 
## , , 2
## 
##      [,1] [,2] [,3]
## [1,]    7    9    1
## [2,]    8   10    2
## 
## , , 3
## 
##      [,1] [,2] [,3]
## [1,]    3    5    7
## [2,]    4    6    8
## 
## , , 4
## 
##      [,1] [,2] [,3]
## [1,]    9    1    3
## [2,]   10    2    4

Aufgabe

Kannst Du anhand des Ergebnisses erkennen, welche Argumente im obigen Befehl was genau beschreiben?*

Teste, ob Du mit einem Array auch mehr als drei Dimensionen generieren kannst.

Wir sehen, dass wir mit obigem Befehl gleich große Matrizen “gestapelt” haben. Und im Nachgang fragen wir gleich noch nach Datentyp bzw. Datenformat, um nochmals den Unterschied zu erkennen. Zusätzlich benutzen wir dim(), um uns sowohl die Anzahl der Dimensionen als auch die Anzahl Dateneinträge pro Dimension anzeigen zu lassen:

typeof(bsp)
## [1] "integer"
class(bsp)
## [1] "array"
dim(bsp)
## [1] 2 3 4

5.1 Rechnen mit Arrays

Oft möchten wir Daten aus einem Array zusammenfassen. Hierzu ist der Befehl apply() sehr mächtig, denn er erlaubt es, bestimmte Operationen bzw. Funktionen entlang einer oder mehrerer Dimensionen anzuwenden:

# Mittelwerte der ersten Dimension des gesamten Arrays
apply(bsp, 1, mean)
## [1] 4.5 5.5
# Mittelwerte der zweiten Dimension des gesamten Arrays
apply(bsp, 2, mean)
## [1] 5.5 5.0 4.5
# Mittelwerte der dritten Dimension des gesamten Arrays
apply(bsp, 3, mean)
## [1] 3.500000 6.166667 5.500000 4.833333
# Punktweise Mittelwerte ENTLANG der dritten Dimension
apply(bsp, c(1, 2), mean)
##      [,1] [,2] [,3]
## [1,]    5  4.5    4
## [2,]    6  5.5    5

Um diese Operationen zu verstehen, ist es nützlich, von den Dimensionen der jeweiligen Ergebnisse auszugehen.

5.2 Elemente eines Arrays

Um bestimmte Element eines Arrays anzusprechen, benutzen wir in R wie für viele andere Datenformate auch die eckigen Klammern [], die auf eine oder mehrere Positionen im Array hinweisen. Das obige Ergebnis enthält hierzu entsprechende die Klammern [], während das Komma-Symbol , die verschiedenen Dimensionen voneinander trennt.

Zur Verdeutlichung folgen nochmals ein paar Beispiele, in denen wir bestimmte Teile eines Arrays isolieren:

(coolarray <- array(c("how", "cool", "is", "this", "!"), c(2, 5, 3)))
## , , 1
## 
##      [,1]   [,2]   [,3]  [,4]   [,5]  
## [1,] "how"  "is"   "!"   "cool" "this"
## [2,] "cool" "this" "how" "is"   "!"   
## 
## , , 2
## 
##      [,1]   [,2]   [,3]  [,4]   [,5]  
## [1,] "how"  "is"   "!"   "cool" "this"
## [2,] "cool" "this" "how" "is"   "!"   
## 
## , , 3
## 
##      [,1]   [,2]   [,3]  [,4]   [,5]  
## [1,] "how"  "is"   "!"   "cool" "this"
## [2,] "cool" "this" "how" "is"   "!"
coolarray[1, , ]
##      [,1]   [,2]   [,3]  
## [1,] "how"  "how"  "how" 
## [2,] "is"   "is"   "is"  
## [3,] "!"    "!"    "!"   
## [4,] "cool" "cool" "cool"
## [5,] "this" "this" "this"
coolarray[ , 2, 2]
## [1] "is"   "this"
coolarray[ , 2:3, ]
## , , 1
## 
##      [,1]   [,2] 
## [1,] "is"   "!"  
## [2,] "this" "how"
## 
## , , 2
## 
##      [,1]   [,2] 
## [1,] "is"   "!"  
## [2,] "this" "how"
## 
## , , 3
## 
##      [,1]   [,2] 
## [1,] "is"   "!"  
## [2,] "this" "how"

5.3 Recycling und Vektorisierung

Die Konzepte Recycling und Vektorisierung sind allgegenwärtig in R und werden oft automatisch und nutzbringend auf multidimensionale Datenformate angewandt. Das Konzept der Vektorisierung erlaubt es uns, mit Datenformaten unterschiedlicher Länge zu rechnen. Wir können z.B. mühelos eine Zahl mit einem Vektor multiplizieren:

2 * c(1, 2, 6)
## [1]  2  4 12

Ebenso gilt das für Vektoren gleicher Länge.

c(1, 2, 3) / c(4, 5, 6)
## [1] 0.25 0.40 0.50

Was passiert aber, wenn die Länge der beiden Vektoren unterschiedlich ist?

c(3, 4) * c(1, 1.2, 5)
## Warning in c(3, 4) * c(1, 1.2, 5): Länge des längeren Objektes
##       ist kein Vielfaches der Länge des kürzeren Objektes
## [1]  3.0  4.8 15.0

R warnt uns, bietet uns aber dennoch ein Ergebnis an.

  • Wie hat R dieses Ergebnis erstellt? Was sind die Vor- und Nachteile dieser Operation?

Das Auffüllen von fehlenden Einträgen im kürzeren Vektor mit bestehenden Werten in der ursprünglichen Reihenfolge wird als Recycling bezeichnet. Recycling und Vektorisierung sind sehr effektiv, können aber auch ungeplante Ergebnisse hervorrufen und dabei schwierig zu findende Bugs in Euer Skript bringen. Versucht doch einmal, die Ergebnisse der folgenden Befehle vorherzusagen, bevor Ihr Eure Schätzung überprüft.

a <- c(1, 2, 3)
b <- c(4, 5)
array(a, b)
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    2    3    1    2
## [2,]    2    3    1    2    3
## [3,]    3    1    2    3    1
## [4,]    1    2    3    1    2

Mit nur einer kleinen Änderung in diesem Befehl erhalten wir ein völlig anderes Ergebnis:

a <- c(1, 2, 3)
b <- c(3, 4, 5)
array(a, b)
## , , 1
## 
##      [,1] [,2] [,3] [,4]
## [1,]    1    1    1    1
## [2,]    2    2    2    2
## [3,]    3    3    3    3
## 
## , , 2
## 
##      [,1] [,2] [,3] [,4]
## [1,]    1    1    1    1
## [2,]    2    2    2    2
## [3,]    3    3    3    3
## 
## , , 3
## 
##      [,1] [,2] [,3] [,4]
## [1,]    1    1    1    1
## [2,]    2    2    2    2
## [3,]    3    3    3    3
## 
## , , 4
## 
##      [,1] [,2] [,3] [,4]
## [1,]    1    1    1    1
## [2,]    2    2    2    2
## [3,]    3    3    3    3
## 
## , , 5
## 
##      [,1] [,2] [,3] [,4]
## [1,]    1    1    1    1
## [2,]    2    2    2    2
## [3,]    3    3    3    3

Trotz der vielen Anwendungsmöglichkeiten von Arrays sind viele moderne geo- und umweltwissenschaftliche Daten in ihrem eigenen Format abgespeichert, v.a. um sowohl die Daten als Metadaten strukturiert und nachvollziehbar aufzubewahren. Geographische Informationssysteme bieten überdies eine Reihe von Datenformaten, deren Austausch über verschiedene Computerbetriebssysteme mittlerweile dank bestimmter Standards möglich ist.

6 Anwendungsbeispiel: Schneeklimatologie in Europa mit ERA5

ERA5 ist das jüngste globale Reanalyseprodukt des ECMWF (European Centre for Medium-Range Weather Forecasts). Aber was ist eigentlich eine Reanalyse? Schaut Euch bitte zunächst das folgende Video an: https://www.youtube.com/embed/FAGobvUGl24

Wir verwenden in diesem Kurs “ERA5-Land”, eine räumlich höher aufgelöstes Spezialprodukt, welches Klimavariablen für die Landoberfläche bereitstellt. Die Daten sind zwar stündlich und ab 1981 gemessen, aber wir nutzen das Produkt mit monatlichen Mittelwerten, um die Datenmenge zu reduzieren. Auch nutzen wir nur ein Fenster über Europa.

Die Daten können hier heruntergeladen werden - Ihr müsst Euch für den Download allerdings bei Copernicus registrieren. Aber das Formular für die Auswahl der Daten könnt Ihr Euch auch so anschauen. Eine Beschreibung der verfügbaren Parameter aus ERA5-Land gibt es hier.

Wir haben für unser Beispiel die Parameter Schneebedeckung (snowc) und Temperatur in 2m Höhe (t2m) heruntergeladen. Das sind für Europa allein immerhin schon über 300 MB.

6.1 Lesen der Daten aus einer NetCDF-Datei

Ohne hier in die Tiefe zu gehen, lernen wir ein wichtiges Datenformat zur Speicherung mehrdimensionaler Daten kennen, das netcdf-Format. Wir nutzen dazu die Bibliothek ncdf4.

library(ncdf4)
era5 <- nc_open("../data/era5_2mtemp_and_snowcover.nc")

Der Befehl nc_open() sollte hier selbsterklärend sein. Das Besondere am NetCDF-Format ist, dass die Daten bereits entlang der Dimensionen strukturiert sind. Diese Dimesionen und weitere Metadaten kann man sich wie folgt anzeigen lassen:

era5
## File ../data/era5_2mtemp_and_snowcover.nc (NC_FORMAT_64BIT):
## 
##      2 variables (excluding dimension variables):
##         short t2m[longitude,latitude,time]   
##             scale_factor: 0.000966278598202075
##             add_offset: 280.116757766463
##             _FillValue: -32767
##             missing_value: -32767
##             units: K
##             long_name: 2 metre temperature
##         short snowc[longitude,latitude,time]   
##             scale_factor: 0.00152594875864068
##             add_offset: 49.9992370256207
##             _FillValue: -32767
##             missing_value: -32767
##             units: %
##             long_name: Snow cover
## 
##      3 dimensions:
##         longitude  Size:401 
##             units: degrees_east
##             long_name: longitude
##         latitude  Size:401 
##             units: degrees_north
##             long_name: latitude
##         time  Size:485 
##             units: hours since 1900-01-01 00:00:00.0
##             long_name: time
##             calendar: gregorian
## 
##     2 global attributes:
##         Conventions: CF-1.6
##         history: 2021-08-18 09:39:00 GMT by grib_to_netcdf-2.20.0: /opt/ecmwf/mars-client/bin/grib_to_netcdf -S param -o /cache/data9/adaptor.mars.internal-1629279532.0576615-3357-15-5ee5f2ed-8159-4098-bf27-eb40d97538ce.nc /cache/tmp/5ee5f2ed-8159-4098-bf27-eb40d97538ce-adaptor.mars.internal-1629279143.1965292-3357-6-tmp.grib

Wir haben also die drei Dimensionen latitude, longitude und time. Entlang dieser Dimensionen kann eine NetCDF-Datei beliebig viele Variablen (also Dateneinträge eines bestimmten Typs) bereithalten. Wie sehen, dass unsere Datei zwei Variablen mit den Namen t2m und snowc besitzt.

Ein paar kleine Testfragen zum Verständnis:

  • Welche räumliche Auflösung haben die Daten?

  • Ab wann beginnt diese Datenreihe?

  • Gibt es Datenlücken? Wenn ja, wie erkennen wir diese?

Eine weitere Möglichkeit, durch die Details einer NetCDF-Datei zu navigieren bietet das $-Zeichen, das man an den Dateinamen (bzw. die darin eingefügten Ebenen) anfügt, z.B.:

era5$var$snowc$varsize
## [1] 401 401 485

Wir können ncvar_get() benutzen, um die Werte gewünschter Variablen zu erhalten. Hierbei entdecken wir eine weitere Eigenheit der NetCDF-Datei, nämlich dass die Koordinaten von Nord nach Süd, bzw. von West nach Ost angeordnet sind:

lat <- ncvar_get(era5, "latitude")
head(lat)  # zeigt nur die ersten paar Werte
## [1] 70.0 69.9 69.8 69.7 69.6 69.5
lon <- ncvar_get(era5, "longitude")
head(lon)  # zeigt nur die ersten paar Werte
## [1] -10.0  -9.9  -9.8  -9.7  -9.6  -9.5

6.2 Analyse der Schneebedeckung

Um einen groben Überblick über die räumlichen Verhältnisse zu erhalten, berechnen wir zunächst den Mittelwert der Schneebedeckung (snowc) über den gesamten Zeitraum für jede Gitterzelle. Da die Zeit die erste Dimension ist, berechnen wir den Mittelwert (mean) entlang der ersten Dimension. Mit ncvar_get() sprechen wir die Dateneinträge direkt an und erhalten diese dem Array-Format, das wir oben kennengelernt haben:

class(ncvar_get(era5, "snowc"))
## [1] "array"
dim(ncvar_get(era5, "snowc"))
## [1] 401 401 485

Da uns nun die Struktur der Daten bewusst ist, können wir nun die lokalen Mittelwerte der Schneebedeckung [%] berechnen:

snowc_time_avg <- apply(ncvar_get(era5, "snowc"), c(1, 2), mean, na.rm = TRUE)

Eine Reihe interessanter Möglichkeiten zum Plotten regelmäßiger Gitter (= Raster) bietet die Bibliothek ggplot2. Zuvor müssen wir allerdings unsere durchschnittlichen Schneehöhen umformen. Dazu nehmen wir passenderweise den Befehl melt() aus der Bibliothek reshape2; dieser Befehl macht aus einer Matrix einen Dataframe, den ggplot() erwartet. Auch müssen wir vorher noch die Reihen und Spalten in unserer Matrix tauschen, da NetCDF die Daten ursprünglich in einer anderen Richtung anlegt (s.o.). Mit expand.grid() können wir die Gitterkoordinaten durch die entsprechenden geographischen Kooordinaten ersetzen:

library(ggplot2)
library(reshape2)
## Warning: Paket 'reshape2' wurde unter R Version 4.2.3 erstellt
# Umformen von Matrix in Dataframe
mean_snow <- melt(t(snowc_time_avg), 
            c("Latitude", "Longitude"), value.name = "Percent")

# Ersetzen der Zellennummern durch echte Koordinatem
mean_snow$Latitude <- expand.grid(lat, lon)[ , 1]
mean_snow$Longitude <- expand.grid(lat, lon)[ , 2]

# Plot
ggplot(mean_snow, aes(Longitude, Latitude, fill = Percent)) +
  geom_tile() +
  labs(title = "Average snow cover (%) from 1981 to 2021",
       x = "Longitude E (°)", y = "Latitude N (°)") +
  theme(aspect.ratio = 1)

Aufgabe

Stelle nun den Zusammenhang zwischen Breitengrad und der langjährig mittleren Schneebedeckung dar, mit anderen Worten: Berechne den Mittelwert der Schneebedeckung entlang der Dimensionen time und longitude und plotte dann den resultierenden eindimensionalen Array gegen die Dimension latitude.

mean_snow_by_lat <- apply(ncvar_get(era5, "snowc"), 2, mean, na.rm = TRUE)

plot(lat, mean_snow_by_lat, type = "l", col = "blue",
    xlab = "Latitude N (°)", 
    ylab = "Average snow cover (%) 1981-2021", cex.lab = 1.2)
grid()

7 Geht das auch anders?

In R gibt es viele Bibliotheken, die auf die Darstellung und Analyse räumlicher (und daher auch multidimensionaler) Daten spezialisiert sind. Alternativ könnten wir aus der Bibliothek raster auch den Befehl brick() benutzen, um aus der NetCDF-Datei einen Stapel aus Rastern zu generieren. Jedes dieser Raster würde dann eine Zeitscheibe aus unseren Daten darstellen:

library(raster)
## Lade nötiges Paket: sp
alt_snow <- brick("../data/era5_2mtemp_and_snowcover.nc", varname = "snowc")
## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/francke/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/francke/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

Dieses Datenformat sieht wie folgt strukturiert aus und wird in R auch oft als Multiband-Raster bezeichnet:

alt_snow
## class      : RasterBrick 
## dimensions : 401, 401, 160801, 485  (nrow, ncol, ncell, nlayers)
## resolution : 0.1, 0.1  (x, y)
## extent     : -10.05, 30.05, 29.95, 70.05  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : era5_2mtemp_and_snowcover.nc 
## names      : X1981.01.01, X1981.02.01, X1981.03.01, X1981.04.01, X1981.05.01, X1981.06.01, X1981.07.01, X1981.08.01, X1981.09.01, X1981.10.01, X1981.11.01, X1981.12.01, X1982.01.01, X1982.02.01, X1982.03.01, ... 
## Date/time  : 1981-01-01, 2021-05-01 (min, max)
## varname    : snowc

Für Satellitenbilder, die oft viele Bänder, also Bildkanäle für bestimmte Wellenlängenbereiche elektromagnetischer Strahlung wie z.B. Licht, Radar, oder Mikrowellen haben, ist dieses Format sehr nützlich. Wir können dieses Objekt vom Typ RasterBrick einfach visualisieren. Um allerdings zu vermeiden, dass wir unbeabsichtigt 485 Zeitschritte plotten, benutzen wir subset(), um uns z.B. die ersten neun Zeitschritte (also die Monate Januar bis September 1981) herauszufiltern:

plot(subset(alt_snow, 1:9))

Der Befehl subset() macht letztlich das Gleiche wie die eckigen Klammern [], denn unsere Daten können wir immer noch als Array auffassen.