Ekosystem(y)

Z prezentacji Rogera Bivanda 'A practical history of R'

Ekosystem(y)

Z prezentacji Rogera Bivanda 'A practical history of R-sig-geo'

Ekosystem(y)

Ekosystem(y)

Ekosystem(y)

Wymagania wstępne

pkgs = c(
  "sf",                
  "raster",            
  "dplyr",            
  "tmap",              
  "landscapemetrics",
  "landscapetools",
  "usethis"
)
to_install = !pkgs %in% installed.packages()
if(any(to_install)) {
  install.packages(pkgs[to_install])
}

usethis::use_course("http://bit.ly/biogis19")

Zadania

1. Zastanów się w czym R może tobie pomóc? Dlaczego chcesz się nauczyć go używać do analiz geoinformacyjnych?
2. Przejrzyj CRAN Task View: Analysis of Spatial Data i sprawdź czy istnieją tam pakiety, które mogą się Tobie przydać w codziennej pracy?

Reprezentacja danych przestrzennych

sf:

• Pakiet sf to następca pakietu sp oparty o OGC standard Simple Features.
• Łączy on możliwości zarówno pakietu sp jak i pakietów rgdal i rgeos.
• Większość jego funkcji zaczyna się od prefiksu st_.
• Ten pakiet pozwala na obsługę dodatkowych typów danych wektorowych (np. poligon i multipoligon to dwie oddzielne klasy), łatwiejsze przetwarzanie danych, oraz obsługę przestrzennych baz danych takich jak PostGIS.
• https://r-spatial.github.io/sf/ oraz https://github.com/rstudio/cheatsheets/blob/master/sf.pdf

raster:

• Pakiet raster zawiera klasy i metody reprezentujące obiekty i operacje rastrowe.
• Pozwala on na wczytywanie i zapisywanie danych rastrowych.
• Umożliwia algebrę rastrów i przetwarzanie rastrów.
• W jego skład wchodzi szereg dodatkowych funkcji, np. do analiz charakterystyk terenu.
• Pozwala na pracę na dużych zbiorach danych.
• ?raster-package oraz http://www.rpubs.com/etiennebr/visualraster

Wczytywanie danych przestrzennych - wektor

library(sf)
polska = st_read("data/polska.gpkg")

## Reading layer `polska' from data source `/home/jn/Science/conferences/BioGIS_19/workshop-pl/data/polska.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 3 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 171842.5 ymin: 132323.8 xmax: 861798.4 ymax: 775305.9
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

plot(polska)

Struktura obiektów przestrzennych - wektor

polska

## Simple feature collection with 1 feature and 3 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 171842.5 ymin: 132323.8 xmax: 861798.4 ymax: 775305.9
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
##     name  pop_est gdp_md_est                           geom
## 1 Poland 38476269    1052000 POLYGON ((487928.3 182549, ...

Wczytywanie danych przestrzennych - raster

library(raster)
dem = raster("data/polska_srtm.tif")

plot(dem)

Struktura obiektów przestrzennych - raster

dem

## class      : RasterLayer 
## dimensions : 702, 1202, 843804  (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333  (x, y)
## extent     : 14.125, 24.14167, 48.99167, 54.84167  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /home/jn/Science/conferences/BioGIS_19/workshop-pl/data/polska_srtm.tif 
## names      : polska_srtm 
## values     : -9, 2220  (min, max)

Zadania

1. Wczytaj do R plik wlkp_pn.gpkg zawierający granice Wielkopolskiego Parku Narodowego. Wyświetl ten obiekt i przejrzyj jego strukturę. Co możesz powiedzieć o zawartości tego pliku? Jaki typ danych on przechowuje? Jaki ma układ współrzędnych? Ile atrybutów zawiera? Jaka jest jego geometria?
2. Wczytaj do R plik wlkp_pn_dem.tif zawierający cyfrowy model wysokościowy dla obszaru Wielkopolskiego Parku Narodowego. Wyświetl ten obiekt i przejrzyj jego strukturę. Co możesz powiedzieć o zawartości tego pliku? Jaki typ danych on przechowuje? Jaki ma układ współrzędnych? Ile atrybutów zawiera? Jakie ma wymiary i rozdzielczość?

Wizualizacja przestrzenna

library(tmap)

• Pakiet służący do tworzenia map tematycznych
• Pozwala on na budowanie map statycznych, map animowanych oraz map interaktywnych
• Jego obsługa oparta jest o dodawanie kolejnych warstw do zwizualizowania a następnie ich modyfikacji
• https://github.com/mtennekes/tmap - lista dodatkowych materiałów do pakietu tmap

Wizualizacja przestrzenna

tm_shape(polska) + 
  tm_polygons()

Wizualizacja przestrzenna

tm_shape(polska) + 
  tm_polygons() + 
  tm_scale_bar(position = c("left", 
                            "bottom"))

Wizualizacja przestrzenna

tm_shape(polska) + 
  tm_polygons() + 
  tm_scale_bar(position = c("left", 
                            "bottom")) + 
  tm_compass()

Wizualizacja przestrzenna

tm_shape(polska) + 
  tm_polygons() + 
  tm_scale_bar(position = c("left", 
                            "bottom")) + 
  tm_compass() +
  tm_layout(title = "Polska")

Wizualizacja przestrzenna

tmap_mode("view")

tm_shape(polska) + 
  tm_polygons() + 
  tm_layout(title = "Polska")

tmap_mode("plot")

Wizualizacja przestrzenna

tm_shape(dem) + 
  tm_raster()

Wizualizacja przestrzenna

tm_shape(dem) + 
  tm_raster(style = "cont",
            midpoint = NA)

Wizualizacja przestrzenna

tm_shape(dem) + 
  tm_raster(style = "cont", 
            midpoint = NA,
            palette = "-RdYlGn")

tmaptools::palette_explorer()

Wizualizacja przestrzenna

tm_shape(dem) + 
  tm_raster(style = "fixed",
            breaks = c(-99, 0, 300,
                       600, 9999), 
            midpoint = NA,
            palette = "-RdYlGn",
            title = "") + 
  tm_layout(legend.position = c("LEFT", 
                                "BOTTOM"))

Wizualizacja przestrzenna

tm_shape(dem) + 
  tm_raster(style = "fixed",
            breaks = c(-99, 0, 300,
                       600, 9999), 
            labels = c("Depresje", 
                       "Niziny", 
                       "Wyżyny",
                       "Góry"),
            midpoint = NA,
            palette = "-RdYlGn",
            title = "") + 
  tm_layout(legend.position = c("LEFT", 
                                "BOTTOM"))

Wizualizacja przestrzenna

tm_shape(dem) + 
  tm_raster(style = "fixed",
            breaks = c(-99, 0, 300,
                       600, 9999), 
            labels = c("Depresje",
                       "Niziny",
                       "Wyżyny", 
                       "Góry"),
            midpoint = NA,
            palette = c("#5E8B73",
                        "#DAE97A", 
                        "#EADC70", 
                        "#AF8D5C"),
            title = "") + 
  tm_layout(legend.position = c("LEFT",
                                "BOTTOM"))

Wizualizacja przestrzenna

mapa1 = tm_shape(dem) + 
  tm_raster(style = "fixed",
            breaks = c(-99, 0, 300, 600, 9999),
            labels = c("Depresje", "Niziny", "Wyżyny", "Góry"),
            midpoint = NA, 
            palette = c("#5E8B73", "#DAE97A", "#EADC70", "#AF8D5C"),
            title = "") +
  tm_layout(legend.position = c("LEFT", "BOTTOM"))

tmap_save(mapa1, "moja_mapa1.png")

Zadania

1. Ulepsz obiekt mapa1 poprzez np. dodanie skali, strzałki północy, czy tytułu. Spróbuj też dodać granicę Polski do tej mapy.
2. Wczytaj do R pliki wlkp_pn.gpkg oraz wlkp_pn_dem.tif. Stwórz mapę przedstawiającą model terenu dla dla obszaru Wielkopolskiego Parku Narodowego. Zapisz uzyskaną mapę do pliku "WLKP_NAZWISKO.png".
3. Dodatkowo: spróbuj powtórzyć kroki z artykułu Geocomputation with R: maps extended w celu wykonania mapy Polski z rzeźbą terenu wykonaną metodą cieniowania.

Przetwarzanie danych

library(dplyr)

• Pakiet dplyr pozwala na łatwe przetwarzanie danych występujących w postaci tabelarycznej
• Jego możliwości obejmują wybieranie konkretnych kolumn czy wierszy, dodawane nowych kolumn, czy wykonywanie podsumowań
• Więcej informacji można znaleźć na https://dplyr.tidyverse.org/ oraz https://r4ds.had.co.nz/transform.html

Przetwarzanie atrybutów

dane_meteo = read.csv("data/polska_meteo_2017.csv")
head(dane_meteo)

Przetwarzanie atrybutów

stacje_meteo = st_read("data/polska_stacje.gpkg")

## Reading layer `polska_stacje' from data source `/home/jn/Science/conferences/BioGIS_19/workshop-pl/data/polska_stacje.gpkg' using driver `GPKG'
## Simple feature collection with 1998 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 172691.6 ymin: 139488.3 xmax: 854361.7 ymax: 846760.3
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +units=m +no_defs

plot(st_geometry(stacje_meteo))

Przetwarzanie atrybutów

dane_meteo_wyb = dane_meteo %>% 
  filter(Nazwa.stacji %in% unique(stacje_meteo$NAZWA_ST))

Przetwarzanie atrybutów

dane_meteo_wyb = dane_meteo %>% 
  filter(Nazwa.stacji %in% unique(stacje_meteo$NAZWA_ST)) %>% 
  mutate(data = as.Date(paste0(Rok, "-", Miesiac, "-", Dzien)))

Przetwarzanie atrybutów

dane_meteo_wyb = dane_meteo %>% 
  filter(Nazwa.stacji %in% unique(stacje_meteo$NAZWA_ST)) %>% 
  mutate(data = as.Date(paste0(Rok, "-", Miesiac, "-", Dzien))) %>% 
  dplyr::select(-Rok, -Miesiac, -Dzien)

Przetwarzanie atrybutów

dane_meteo_wyb = dane_meteo %>% 
  filter(Nazwa.stacji %in% unique(stacje_meteo$NAZWA_ST)) %>% 
  mutate(data = as.Date(paste0(Rok, "-", Miesiac, "-", Dzien))) %>% 
  dplyr::select(-Rok, -Miesiac, -Dzien) %>% 
  summarize(tavg = mean(tavg), pressure = mean(pressure))

Przetwarzanie atrybutów

dane_meteo_wyb = dane_meteo %>% 
  filter(Nazwa.stacji %in% unique(stacje_meteo$NAZWA_ST)) %>% 
  mutate(data = as.Date(paste0(Rok, "-", Miesiac, "-", Dzien))) %>% 
  dplyr::select(-Rok, -Miesiac, -Dzien) %>% 
  group_by(Kod.stacji) %>% 
  summarize(tavg = mean(tavg), pressure = mean(pressure))

Łączenie danych

meteo = stacje_meteo %>% 
  left_join(dane_meteo_wyb, by = c("KOD_SZS" = "Kod.stacji"))

Łączenie danych

meteo = stacje_meteo %>% 
  inner_join(dane_meteo_wyb, by = c("KOD_SZS" = "Kod.stacji"))

Reprojekcje

plot(dem)

plot(st_geometry(meteo), axes = TRUE)

Reprojekcje

meteo = stacje_meteo %>% 
  inner_join(dane_meteo_wyb, by = c("KOD_SZS" = "Kod.stacji")) %>% 
  st_transform(4326)

Reprojekcje

tm_shape(dem) +
  tm_raster() +
  tm_shape(meteo) +
  tm_symbols(col = "tavg", palette = "RdBu")

Interakcja wektor-raster

meteo$elev = extract(dem, meteo)
head(meteo)

## Simple feature collection with 6 features and 5 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 19.96032 ymin: 49.23253 xmax: 22.3417 ymax: 50.81048
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##     KOD_SZS        NAZWA_ST      tavg pressure                      geom
## 1 349190625        ZAKOPANE 6.5290411    0.000 POINT (19.96032 49.29382)
## 2 349190650 KASPROWY WIERCH 0.1978082    0.000 POINT (19.98182 49.23253)
## 3 349200660       NOWY SĄCZ 9.2180822 1014.784  POINT (20.6886 49.62714)
## 4 349210670          KROSNO 8.7728767 1017.553 POINT (21.76917 49.70674)
## 5 349220690           LESKO 8.3227397 1017.737  POINT (22.3417 49.46647)
## 6 350200570    KIELCE-SUKÓW 8.4849315 1016.895 POINT (20.69221 50.81048)
##   elev
## 1  844
## 2 1796
## 3  275
## 4  323
## 5  384
## 6  247

Interakcja raster-wektor

lc = raster("data/polska_lc.tif")
plot(lc)

wlkp_pn = st_read("data/wlkp_pn.gpkg", quiet = TRUE)
plot(st_geometry(wlkp_pn), axes = TRUE)

Interakcja raster-wektor

wlkp_pn_lc = crop(lc, wlkp_pn)
plot(wlkp_pn_lc)

Interakcja raster-wektor

wlkp_pn_lc2 = mask(wlkp_pn_lc, wlkp_pn)
plot(wlkp_pn_lc2)

Zadania

1. Wczytaj dane meteorologiczne dla roku 2017 oraz dane o położeniu stacji meteorologicznych. Stwórz obiekt przestrzenny zawierający wszystkie pomiary dla trzech stacji: BIAŁYSTOK, KRAKÓW-BALICE, oraz POZNAŃ. (Porada: użyj do tego, między innymi funkcji right_join()).
2. Wylicz średnią temperaturę w styczniu oraz w lipcu dla tych trzech stacji.
3. Określ na jakim pokryciu terenu znajdują się te trzy stacje.
4. Dodatkowe: ulepsz mapę temperatur średnich zaprezentowaną w części "Reprojekcje".

Przykład - metryki krajobrazowe

library(landscapemetrics)
library(landscapetools)

show_landscape(wlkp_pn_lc, discrete = TRUE)

Przykład - metryki krajobrazowe

list_lsm()

## # A tibble: 132 x 5
##    metric name                         type             level function_name
##    <chr>  <chr>                        <chr>            <chr> <chr>        
##  1 area   patch area                   area and edge m… patch lsm_p_area   
##  2 cai    core area index              core area metric patch lsm_p_cai    
##  3 circle related circumscribing circ… shape metric     patch lsm_p_circle 
##  4 contig contiguity index             shape metric     patch lsm_p_contig 
##  5 core   core area                    core area metric patch lsm_p_core   
##  6 enn    euclidean nearest neighbor … aggregation met… patch lsm_p_enn    
##  7 frac   fractal dimension index      shape metric     patch lsm_p_frac   
##  8 gyrate radius of gyration           area and edge m… patch lsm_p_gyrate 
##  9 ncore  number of core areas         core area metric patch lsm_p_ncore  
## 10 para   perimeter-area ratio         shape metric     patch lsm_p_para   
## # … with 122 more rows

Metryki krajobrazowe - krajobraz

• AI - Aggregation index

lsm_l_ai(wlkp_pn_lc)

## # A tibble: 1 x 6
##   layer level     class    id metric value
##   <int> <chr>     <int> <int> <chr>  <dbl>
## 1     1 landscape    NA    NA ai      70.3

Metryki krajobrazowe - klasa

• ED - Edge density

lsm_c_ed(wlkp_pn_lc)

## # A tibble: 12 x 6
##    layer level class    id metric  value
##    <int> <chr> <int> <int> <chr>   <dbl>
##  1     1 class    10    NA ed      8.27 
##  2     1 class    11    NA ed     10.6  
##  3     1 class    30    NA ed      3.12 
##  4     1 class    40    NA ed      0.767
##  5     1 class    60    NA ed      2.35 
##  6     1 class    70    NA ed      9.61 
##  7     1 class    90    NA ed      6.59 
##  8     1 class   100    NA ed      3.65 
##  9     1 class   130    NA ed      1.38 
## 10     1 class   180    NA ed      0.295
## 11     1 class   190    NA ed      5.32 
## 12     1 class   210    NA ed      1.69

Metryki krajobrazowe - płat

• PARA - Perimeter-Area ratio

lsm_p_para(wlkp_pn_lc)

## # A tibble: 234 x 6
##    layer level class    id metric   value
##    <int> <chr> <int> <int> <chr>    <dbl>
##  1     1 patch    10     1 para   0.0169 
##  2     1 patch    10     2 para   0.0169 
##  3     1 patch    10     3 para   0.0117 
##  4     1 patch    10     4 para   0.00620
##  5     1 patch    10     5 para   0.0113 
##  6     1 patch    10     6 para   0.00893
##  7     1 patch    10     7 para   0.0169 
##  8     1 patch    10     8 para   0.0169 
##  9     1 patch    10     9 para   0.0127 
## 10     1 patch    10    10 para   0.0126 
## # … with 224 more rows

Metryki krajobrazowe - zbiorczo

calculate_lsm(wlkp_pn_lc, type = "aggregation metric")

## # A tibble: 428 x 6
##    layer level class    id metric value
##    <int> <chr> <int> <int> <chr>  <dbl>
##  1     1 class    10    NA ai      43.4
##  2     1 class    11    NA ai      85.0
##  3     1 class    30    NA ai      25.7
##  4     1 class    40    NA ai      40.6
##  5     1 class    60    NA ai      63.6
##  6     1 class    70    NA ai      68.7
##  7     1 class    90    NA ai      70.9
##  8     1 class   100    NA ai      33.5
##  9     1 class   130    NA ai      66.7
## 10     1 class   180    NA ai      37.5
## # … with 418 more rows

Metryki krajobrazowe - wizualizacja

mapa_p_para = get_lsm(wlkp_pn_lc, what = "lsm_p_para")

## Warning: only using 'what' argument

show_landscape(mapa_p_para[[1]][[1]])

Metryki krajobrazowe - wizualizacja

show_lsm(wlkp_pn_lc, what = "lsm_p_para")

Metryki krajobrazowe - próbki

punkt = st_sf(st_sfc(geom = st_point(c(347500, 493000))))
plot(wlkp_pn_lc)
plot(st_geometry(punkt), add = TRUE, cex = 2)

Metryki krajobrazowe - próbki

sample_lsm(wlkp_pn_lc, punkt, shape = "circle", size = 2000, what = "lsm_c_ed")

## # A tibble: 10 x 8
##    layer level class    id metric  value plot_id percentage_inside
##    <int> <chr> <int> <int> <chr>   <dbl>   <int>             <dbl>
##  1     1 class    10    NA ed     11.5         1              100.
##  2     1 class    11    NA ed     19.9         1              100.
##  3     1 class    30    NA ed      5.26        1              100.
##  4     1 class    60    NA ed      2.14        1              100.
##  5     1 class    70    NA ed      9.34        1              100.
##  6     1 class    90    NA ed      7.50        1              100.
##  7     1 class   100    NA ed      7.80        1              100.
##  8     1 class   180    NA ed      0.795       1              100.
##  9     1 class   190    NA ed      2.14        1              100.
## 10     1 class   210    NA ed      0.795       1              100.

Zadania

1. W pliku polska_pn.gpkg znajdują się granice parków narodowych w Polsce. Wczytaj ten plik do R i wybierz z niego tylko Białowieski Park Narodowy.
2. Przytnij mapę pokrycia terenu z pliku polska_lc.tif do zasięgu (obwiedni) Białowieskiego Parku Narodowego.
3. Jakie kategorie pokrycia terenu można znaleźć na tym obszarze (możesz do tego celu użyć funkcji unique())?
4. Wylicz powierzchnię całego obszaru oraz powierzchnie kolejnych klas (TA i CA). Jakie są trzy kategorie o największej powierzchni?
5. Sprawdź ile płatów znajduje się na tym obszarze (NP).

Źródła informacji