Dane czasoprzestrzenne w R. Część 2

Zbiór danych air: Dane o jakości powietrza w Niemczech

Author

Anna Dmowska

library(tidyverse)
library(sf)
library(spacetime)
library(sp)

1 Zbiór danych `air` z pakietu `spacetime`

Zbiór danych air z pakietu spacetime dostarcza średnie dzienne wartości PM10 dla stacji tła na obszarach wiejskich w Niemczech w latach 1998–2009.

# Ładujemy zbiór 'air', który zawiera pomiary pyłu zawieszonego (PM10) 
# z wiejskich stacji w Niemczech
data(list = "air", package = "spacetime")

#strona pomocy z opisem danych
?air

# Sprawdzamy załadowane obiekty:
ls()

[1] "air"      "dates"    "DE"       "DE_NUTS1" "stations"

Zbiór danych air składa się z 5 obiektów:

stations: współrzędne geograficzne stacji przechowywane jako SpatialPoints.
dates: wektor z datami pomiarów w formacie Date. Dane są dostarczone z krokiem dziennym (1 pomiar na dzień)
air: macierz z surowymi danymi pomiarowymi PM10, składa się z 70 wierszy odpowiadającym stacjom oraz 4383 momenty w czasie (kolumny).
DE: granica Niemiec
DE_NUTS1: granice regionów NUTS1 pomocne przy wizualizacji

Komponent przestrzenny: obiekt stations

head(stations)

SpatialPoints:
        coords.x1 coords.x2
DESH001  9.585911  53.67057
DENI063  9.685030  53.52418
DEUB038  9.791584  54.07312
DEBE056 13.647013  52.44775
DEBE062 13.296353  52.65315
DEBE032 13.225856  52.47309
Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84
+no_defs

Komponent czasowy: obiekt dates

head(dates)

[1] "1998-01-01" "1998-01-02" "1998-01-03" "1998-01-04" "1998-01-05"
[6] "1998-01-06"

Komponent dane: obiekt air

head(air[, 1:5])

        [,1] [,2] [,3] [,4] [,5]
DESH001   NA   NA   NA   NA   NA
DENI063   NA   NA   NA   NA   NA
DEUB038   NA   NA   NA   NA   NA
DEBE056   NA   NA   NA   NA   NA
DEBE062   NA   NA   NA   NA   NA
DEBE032   NA   NA   NA   NA   NA

1.0.1 Dane dodatkowe: DE_NUTS1 oraz DE

plot(DE_NUTS1)

2 Wizualizacja danych

Dane DE_NUTS1 oraz stations zostały przekształcone z klasy sp na obiekty klasy sf w celu ich wizualizacji za pomocą ggplot2.

ggplot() +
  geom_sf(data = st_as_sf(DE_NUTS1)) + 
  geom_sf(data = st_as_sf(stations), color = "red", size = 1) + 
  coord_sf() + 
  theme_bw() +
  labs(x = "Długość geogr.", y = "Szerokość geogr.", title = "Lokalizacja stacji na tle landów")

Druga mapa dodatkowo uwzględnia także etykiety z nazwami landów.

ggplot() +
  # 1. Rysowanie granic landów
  geom_sf(data = st_as_sf(DE_NUTS1)) + 
  # 2. Dodawanie etykiet tekstowych z kolumny NAME_1
  # aes(label = NAME_1) wskazuje, co ma być napisem
  geom_sf_label(data = st_as_sf(DE_NUTS1), 
               aes(label = NAME_1), 
               size = 3,           # wielkość czcionki
               color = "darkblue", # kolor tekstu
               check_overlap = TRUE) + # zapobiega nakładaniu się napisów
  # 3. Rysowanie stacji pomiarowych
  geom_sf(data = st_as_sf(stations), color = "red", size = 1) + 
  coord_sf() + 
  theme_bw() +
  labs(x = "Długość geogr.", y = "Szerokość geogr.", title = "Lokalizacja stacji na tle landów")

3 Przygotowanie danych do analizy geostatystycznej

3.1 Konstrukcja obiektu typu STFDF (Spatio-Temporal Full Data Frame)

Pierwszym zadaniem jest połączenie informacji przestrzennych i czasowych w jedną strukturę danych: obiekt klasy STFDF (ang. space-time full-grid data frame). STFDF (Spatio-Temporal Full Data Frame) to struktura, w której każda stacja posiada pomiar w każdym kroku czasowym.

W tym celu użyjemy funkcji STFDF() z pakietu spacetime, która wymaga podania:

Punktów przestrzennych: obiekt stations,
Czasów obserwacji: obiekt dates,
Zmiennej zmierzonej w każdym punkcie w przestrzeni i w każdym momencie czasu: obiekt air

Uwaga! Liczba punktów przestrzennych oraz dat musi odpowiadać odpowiednio liczbie wierszy i kolumn w macierzy. .

rural = STFDF(sp = stations, time = dates, data = data.frame(PM10 = as.vector(air)))

str(rural)

Formal class 'STFDF' [package "spacetime"] with 4 slots
  ..@ data   :'data.frame': 306810 obs. of  1 variable:
  .. ..$ PM10: num [1:306810] NA NA NA NA NA NA NA NA NA NA ...
  ..@ sp     :Formal class 'SpatialPoints' [package "sp"] with 3 slots
  .. .. ..@ coords     : num [1:70, 1:2] 9.59 9.69 9.79 13.65 13.3 ...
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:70] "DESH001" "DENI063" "DEUB038" "DEBE056" ...
  .. .. .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
  .. .. ..@ bbox       : num [1:2, 1:2] 6.28 47.81 14.79 54.92
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
  .. .. .. .. ..$ : chr [1:2] "min" "max"
  .. .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs"
  .. .. .. .. ..$ comment: chr "GEOGCRS[\"unknown\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.25722"| __truncated__
  ..@ time   :An xts object on 1998-01-01 / 2009-12-31 containing: 
  Data:    integer [4383, 1]
  Columns: timeIndex
  Index:   Date [4383] (TZ: "UTC")
  ..@ endTime: POSIXct[1:4383], format: "1998-01-02" "1998-01-03" ...

3.2 Selekcja danych do analizy

Do analizy wybieramy podzbiór danych z lat 2005-2010, i usuwamy stacje dla których dla tych lat nie wykonano żadnych pomiarów.

# Składnia "od::do" pozwala na szybką selekcję zakresów czasowych w pakietach 'xts' i 'spacetime'
rr = rural[,"2005::2010"]

# Identyfikacja oraz usunięcie stacji bez żadnych pomiarów
unsel = which(apply(as(rr, "xts"), 2, function(x) all(is.na(x))))
r5to10 = rr[-unsel,]

dim(r5to10)

    space      time variables 
       53      1826         1

Po selekcji zbiór danych składa się z 53 punktów pomiarowych, 1826 dziennych pomiarów oraz jednej zmiennej (PM10).

# dodatkowe informacje o obiekcie r5to10
str(r5to10)
head(r5to10@sp)
r5to10@sp@bbox
head(r5to10@time)
head(r5to10@data)

#zapis przygotowanych danych do pliku 
save(r5to10, file = "r5to10.RData")

3.3 Wizualizacja danych

Opcja mode = “xt wyświetla tzw. diagram Hovmollera, gdzie na osi X są lokalizacje, a na osi Y czas, kolor odpowiada wartością wyświetlanej zmiennej. Diagram ten pozwala m.in na identyfikację braków danych oraz trendów w przestrzeni/czasie.

library(RColorBrewer)
cuts = seq(0, 300, 50)
ck = list(at = cuts, labels = as.character(cuts))
stplot(r5to10, mode = "xt", 
       col.regions = rev(brewer.pal(6, "Spectral")), 
       cuts = cuts, colorkey = ck, asp = 0.5,
       scales = list(x = list(rot = 90)),
       xlab = "Lokalizacje", ylab = "Czas")

3.4 Rozkład wartości zmiennej PM10 w stacjach pomiarowych w wybranych dniach

stplot(r5to10[, "2005-03-20::2005-03-25"])

stplot(r5to10[, 400:402])

3.5 Rozkład wartości zmiennej PM10

ggplot(r5to10@data, aes(x = PM10)) + geom_histogram() + theme_bw()

1 Zbiór danych air z pakietu spacetime