Wprowadzenie

Celem tej analizy jest wykonanie interpolacji przestrzennej danych klimatycznych Państwowego Instytutu Badawczego Instytutu Meteorologii i Gospodarki Wodnej (IMGW-PIB) z wykorzystaniem pakietu climate w R. Interpolacja przestrzenna ma na celu określenie wartości zmiennej, w punkcie w którym nie była ona zmierzona. Często stosuje się w analizach i prognozach przestrzennych, opracowaniach map izoliniowych, badaniach klimatologicznych (np. dla opracowania map rozkładu temperatury) itd.

Metody

Najpierw należy wczytać pakiety potrzebne do analizy, w tym pakiet climate, który ma na celu automatyczne pobieranie danych meteorologicznych i hydrologicznych z publicznie dostępnych baz danych (w tym z IMGW).

library(sf)
## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(maps)
library(climate)

Za pomocą funkcji meteo_imgw_monthly pobierzemy dane miesięczne ze stacji synoptycznych dla 2018 roku w następujący sposób:

stacje = meteo_imgw_monthly(rank = "synop", year = 2018, coords = T)
## [1] "https://dane.imgw.pl/data/dane_pomiarowo_obserwacyjne/dane_meteorologiczne/miesieczne/synop/s_m_d_format.txt"
## C:\Users\daria\AppData\Local\Temp\RtmpyqpSMi\file17f4657366a6
## [1] "https://dane.imgw.pl/data/dane_pomiarowo_obserwacyjne/dane_meteorologiczne/miesieczne/synop/s_m_t_format.txt"
## C:\Users\daria\AppData\Local\Temp\RtmpyqpSMi\file17f42ff330
## [1] "https://dane.imgw.pl/data/dane_pomiarowo_obserwacyjne/dane_meteorologiczne/miesieczne/synop/"
## C:\Users\daria\AppData\Local\Temp\RtmpyqpSMi\file17f42bba422f
## [1] "https://dane.imgw.pl/data/dane_pomiarowo_obserwacyjne/dane_meteorologiczne/miesieczne/synop/2018/2018_m_s.zip"
## C:\Users\daria\AppData\Local\Temp\RtmpyqpSMi\file17f47cd45c05

Dodajemy coords = T, ponieważ trzeba będzie dostać wartości współrzędnych stacji, żeby następnie przekształcić je w obiekt przestrzenny. Z uzyskanych danych wybierzemy dane dotyczące liczby dni deszczowych oraz wskaźnika nasłonecznienia (insolacji). Statystyki miesięczne są konwertowane na dane średnie roczne za pomocą funkcji summarize, wyniki są zaokrąglone za pomocą round i zapisywane jako imgw_roczne.

imgw_roczne = stacje %>%
  group_by(station, X, Y) %>%
  summarize(srednia_roczna_liczba_dni_deszczowych = round(mean(rain_days)),
            naslonecznienie = round(mean(insolation_monthly)))
## `summarise()` has grouped output by 'station', 'X'. You can override using the `.groups` argument.

Otrzymane kolumny w tabeli imgw_roczne noszą nazwy srednia_roczna_liczba_dni_deszczowych i naslonecznienie odpowiednio. Aby przekonwertować obiekt imgw_roczne do obiektu przestrzennego można skorzystać z funkcji st_as_sf (coords definiuje kolumny które zawierają współrzędne x i y):

imgw_sf = st_as_sf(imgw_roczne, coords = c(2,3))

Dalej trzeba nadać obiektowi imgw_sf układ współrzędnych WGS-84 za pomocą kodu EPSG 4326 i następnie przetransformować geograficzne współrzędne do kartograficznych, do układu współrzędnych 1992:

st_crs(imgw_sf) = 4326
imgw_92 = st_transform(imgw_sf, 2180)

Teraz można dodać do analizy granice polskich województw, wykorzystając funkcję st_read:

granice_woj = st_read("C:\\Users\\daria\\Desktop\\Studia EGP\\AiWDP\\Cwiczenia\\Interpolacja przestrzenna\\woj\\Wojewodztwa.shp")
## Reading layer `Wojewodztwa' from data source `C:\Users\daria\Desktop\Studia EGP\AiWDP\Cwiczenia\Interpolacja przestrzenna\woj\Wojewodztwa.shp' using driver `ESRI Shapefile'
## Simple feature collection with 16 features and 29 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 171677.5 ymin: 133223.7 xmax: 861895.8 ymax: 775019.1
## Projected CRS: ETRS89 / Poland CS92

Za pomocą funkcji st_union zagregujemy granice województw do granic Polski:

pl = st_union(granice_woj)

Do wizualizacji wyników wykorzystujemy pakiety tmap i RColorBrewer, otrzymane wykresy są przedstawione w następnym rozdziale.

library(tmap)
library(RColorBrewer)

Do interpolacji przestrzennej nasłonecznienia w Polsce wczytamy pakiet stars zapewniający infrastrukturę dla analizy danych rastrowych i pakiet gstat, który zawiera różne modele geostatystyczne, w tym metody interpolacji przestrzennej:

library(stars)
## Loading required package: abind
library(gstat)

Używając funkcji st_as_stars utworzymy pusty raster, zasięg którego odpowiada zasięgu granic Polski, argumenty dx i dy definiują rozmiar pikseli w rastrze. Dalej utworzonemu rastrowi został nadany układ 1992:

raster = st_as_stars(st_bbox(pl),
                   dx = 500,
                   dy = 500)
st_crs(raster) = 2180

Do interpolacji skorzystamy z metody odwrotnych odległości IDW, która wylicza wartość dla każdej komórki na podstawie wartości punktów obokległych ważonych odwrotnością ich odległości. To znaczy, że czym punkt jest bardziej oddalony, tym mniejszy jest jego wpływ na interpolowaną wartość. Wykonujemy interpolację za pomocą funkcji idw (wartość wykładnika potęgowego idp = 2):

slonce_idw = idw(naslonecznienie ~ 1, locations = imgw_92,
                 newdata = raster, idp = 2)
## [inverse distance weighted interpolation]

Obetniemy otrzymany obiekt do granic Polski za pomocą funkcji st_crop:

slonce_idw_crop = st_crop(slonce_idw, pl)

Wyniki interpolacji są przedstawione na wykresie w następnym rozdziale.

Wyniki, wizualizacje i wnioski

Przedstawimy wizualizację obiektu punktowego imgw_92 i granic województw Polski na poniższym wykresie, żeby sprawdzić czy stacje lokalizują się w odpowiednich miejscach.

plot(st_geometry(granice_woj))
plot(st_geometry(imgw_92), add = T)

Dalej zwizualizujemy granice Polski utworzone za pomocą funkcji st_union:

plot(st_geometry(pl))

Wykorzystując pakiet tmap, przedstawimy wykres średniej rocznej liczby dni deszczowych w Polsce w 2018 roku zgodnie z pobranymi z IMGW danymi:

tm_shape(granice_woj)+
  tm_borders()+
  tm_shape(imgw_92)+
  tm_dots(col = "srednia_roczna_liczba_dni_deszczowych", palette = "RdYlGn", size = 0.3,
          title = "Srednia roczna liczba dni deszczowych w Polsce (2018)")+
  tm_text("srednia_roczna_liczba_dni_deszczowych", size = 0.4, xmod = 1.0)+
  tm_layout(legend.title.size = 3, 
            inner.margins = c(0,.0,.05,.0), asp = 1.5,
            legend.outside = T, legend.outside.position = 'bottom',
            legend.frame = T)

W ten sam sposób zwizualizujemy nasłonecznienie w Polsce w 2018 roku:

tm_shape(granice_woj)+
  tm_borders()+
  tm_shape(imgw_92)+
  tm_dots(col = "naslonecznienie", palette = "Set1", size = 0.3, 
          title = "Średnie naslonecznienie w Polsce (2018)")+
  tm_text("naslonecznienie", size = 0.4, xmod = 1.0)+ 
  tm_layout(legend.title.size = 3, 
            inner.margins = c(0,.0,.05,.0), asp = 1.5,
            legend.outside = T, legend.outside.position = 'right',
            legend.frame = T)

Następny wykres pokazuje wyniki interpolacji danych nasłonecznienia metodą IDW.

plot(slonce_idw["var1.pred"], main = "IDW (insolacja)")

Dalej obiekt został obciągnięty do granic Polski i zwizualizowany za pomocą palety kolorów pakietu RColorBrewer:

plot(slonce_idw_crop, col = brewer.pal(10, "RdYlGn"), main = "IDW (insolacja)")

W ramach tej analizy pobrano i zwizualizowano dane IMGW dotyczące średniej liczby dni deszczowych w 2018 roku, a także średniego nasłonecznienia w 2018 roku. Oprócz tego, dane o nasłonecznieniu średnim interpolowano metodą odwrotnych odległości IDW (model deterministyczny), która pozwala na obliczenie brakujących wartości wskaźnika na podstawie średniej ważonej z wartości z pobliskich lokalizacji (gdzie pomiar jest wykonany). Na powstałym wykresie widać interpolowane nasłonecznienie (np. w południowo-wschodniej części Polski wartości insolacji są niższe, niż w południowo-zachodniej).