Z prezentacji Rogera Bivanda 'A practical history of R'
Z prezentacji Rogera Bivanda 'A practical history of R-sig-geo'
Pakiety (moduły) R:
Klasy obiektów przestrzennych:
Wizualizacje:
Klasy obiektów przestrzennych:
Wizualizacje:
Wszystkie pakiety używane w poniższych przykładach można zainstalować używając kodu poniżej:
pkgs = c( "sf", "raster", "dplyr", "tmap", "landscapemetrics", "landscapetools", "usethis")to_install = !pkgs %in% installed.packages()if(any(to_install)) { install.packages(pkgs[to_install])}
Do pobrania materiałów na te warsztaty służy poniższy kod:
usethis::use_course("http://bit.ly/biogis19")
sf:
st_
.raster:
raster-package
oraz http://www.rpubs.com/etiennebr/visualrasterlibrary(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
st_write()
.plot(polska)
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, ...
library(raster)dem = raster("data/polska_srtm.tif")
writeRaster()
.plot(dem)
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)
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?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ść?library(tmap)
tm_shape(polska) + tm_polygons()
tm_shape(polska) + tm_polygons() + tm_scale_bar(position = c("left", "bottom"))
tm_shape(polska) + tm_polygons() + tm_scale_bar(position = c("left", "bottom")) + tm_compass()
tm_shape(polska) + tm_polygons() + tm_scale_bar(position = c("left", "bottom")) + tm_compass() + tm_layout(title = "Polska")
tmap_mode("view")
tm_shape(polska) + tm_polygons() + tm_layout(title = "Polska")
tmap_mode("plot")
tm_shape(dem) + tm_raster()
tm_shape(dem) + tm_raster(style = "cont", midpoint = NA)
tm_shape(dem) + tm_raster(style = "cont", midpoint = NA, palette = "-RdYlGn")
tmaptools::palette_explorer()
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"))
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"))
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"))
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")
mapa1
poprzez np. dodanie skali, strzałki północy, czy tytułu.
Spróbuj też dodać granicę Polski do tej mapy.
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".library(dplyr)
dane_meteo = read.csv("data/polska_meteo_2017.csv")head(dane_meteo)
Kod.stacji | Nazwa.stacji | Rok | Miesiac | Dzien | tavg | precip | vapor_pres | humidity | pressure |
---|---|---|---|---|---|---|---|---|---|
354150100 | KOŁOBRZEG | 2017 | 1 | 1 | 3.6 | 4.6 | 7.4 | 92.5 | 1013.5 |
354150100 | KOŁOBRZEG | 2017 | 1 | 2 | 3.5 | 4.7 | 6.0 | 76.9 | 1011.1 |
354150100 | KOŁOBRZEG | 2017 | 1 | 3 | 1.9 | 10.0 | 6.7 | 92.3 | 1007.1 |
354150100 | KOŁOBRZEG | 2017 | 1 | 4 | 2.7 | 2.3 | 5.8 | 76.0 | 993.9 |
354150100 | KOŁOBRZEG | 2017 | 1 | 5 | -1.7 | 6.2 | 3.9 | 74.5 | 1024.3 |
354150100 | KOŁOBRZEG | 2017 | 1 | 6 | -4.5 | 4.0 | 3.5 | 77.3 | 1038.7 |
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))
dane_meteo_wyb = dane_meteo %>% filter(Nazwa.stacji %in% unique(stacje_meteo$NAZWA_ST))
Kod.stacji | Nazwa.stacji | Rok | Miesiac | Dzien | tavg | precip | vapor_pres | humidity | pressure |
---|---|---|---|---|---|---|---|---|---|
354160105 | KOSZALIN | 2017 | 1 | 1 | 3.4 | 7.0 | 7.5 | 95.4 | 1013.5 |
354160105 | KOSZALIN | 2017 | 1 | 2 | 2.3 | 0.9 | 5.8 | 80.6 | 1010.7 |
354160105 | KOSZALIN | 2017 | 1 | 3 | 1.4 | 8.0 | 6.4 | 93.0 | 1007.0 |
354160105 | KOSZALIN | 2017 | 1 | 4 | 2.3 | 6.9 | 5.6 | 77.0 | 992.7 |
354160105 | KOSZALIN | 2017 | 1 | 5 | -2.9 | 1.5 | 2.8 | 55.6 | 1023.6 |
354160105 | KOSZALIN | 2017 | 1 | 6 | -5.7 | 0.9 | 3.2 | 79.5 | 1038.8 |
dane_meteo_wyb = dane_meteo %>% filter(Nazwa.stacji %in% unique(stacje_meteo$NAZWA_ST)) %>% mutate(data = as.Date(paste0(Rok, "-", Miesiac, "-", Dzien)))
Kod.stacji | Nazwa.stacji | Rok | Miesiac | Dzien | tavg | precip | vapor_pres | humidity | pressure | data |
---|---|---|---|---|---|---|---|---|---|---|
354160105 | KOSZALIN | 2017 | 1 | 1 | 3.4 | 7.0 | 7.5 | 95.4 | 1013.5 | 2017-01-01 |
354160105 | KOSZALIN | 2017 | 1 | 2 | 2.3 | 0.9 | 5.8 | 80.6 | 1010.7 | 2017-01-02 |
354160105 | KOSZALIN | 2017 | 1 | 3 | 1.4 | 8.0 | 6.4 | 93.0 | 1007.0 | 2017-01-03 |
354160105 | KOSZALIN | 2017 | 1 | 4 | 2.3 | 6.9 | 5.6 | 77.0 | 992.7 | 2017-01-04 |
354160105 | KOSZALIN | 2017 | 1 | 5 | -2.9 | 1.5 | 2.8 | 55.6 | 1023.6 | 2017-01-05 |
354160105 | KOSZALIN | 2017 | 1 | 6 | -5.7 | 0.9 | 3.2 | 79.5 | 1038.8 | 2017-01-06 |
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)
Kod.stacji | Nazwa.stacji | tavg | precip | vapor_pres | humidity | pressure | data |
---|---|---|---|---|---|---|---|
354160105 | KOSZALIN | 3.4 | 7.0 | 7.5 | 95.4 | 1013.5 | 2017-01-01 |
354160105 | KOSZALIN | 2.3 | 0.9 | 5.8 | 80.6 | 1010.7 | 2017-01-02 |
354160105 | KOSZALIN | 1.4 | 8.0 | 6.4 | 93.0 | 1007.0 | 2017-01-03 |
354160105 | KOSZALIN | 2.3 | 6.9 | 5.6 | 77.0 | 992.7 | 2017-01-04 |
354160105 | KOSZALIN | -2.9 | 1.5 | 2.8 | 55.6 | 1023.6 | 2017-01-05 |
354160105 | KOSZALIN | -5.7 | 0.9 | 3.2 | 79.5 | 1038.8 | 2017-01-06 |
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))
tavg | pressure |
---|---|
8.685298 | 960.5246 |
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))
Kod.stacji | tavg | pressure |
---|---|---|
349190600 | 9.2780822 | 1017.513 |
349190625 | 6.5290411 | 0.000 |
349190650 | 0.1978082 | 0.000 |
349200660 | 9.2180822 | 1014.784 |
349210670 | 8.7728767 | 1017.553 |
349220690 | 8.3227397 | 1017.737 |
meteo = stacje_meteo %>% left_join(dane_meteo_wyb, by = c("KOD_SZS" = "Kod.stacji"))
KOD_SZS | NAZWA_ST | tavg | pressure | geom |
---|---|---|---|---|
249190020 | BRZEŹNICA | NA | NA | c(545976.260510018, 233149.457458902) |
249199999 | OSIEK | NA | NA | c(519109.405047427, 231411.908795745) |
249190040 | GIERAŁTOWICE | NA | NA | c(527197.290270776, 230151.641154815) |
249190050 | RADZISZÓW | NA | NA | c(558306.218144104, 230885.909302064) |
249190060 | SIEPRAW | NA | NA | c(569259.6149073, 226821.435334128) |
249190070 | KĘTY | NA | NA | c(516182.656254972, 223533.556193128) |
meteo = stacje_meteo %>% inner_join(dane_meteo_wyb, by = c("KOD_SZS" = "Kod.stacji"))
KOD_SZS | NAZWA_ST | tavg | pressure | geom |
---|---|---|---|---|
349190625 | ZAKOPANE | 6.5290411 | 0.000 | c(569804.515915351, 158924.036586795) |
349190650 | KASPROWY WIERCH | 0.1978082 | 0.000 | c(571455.865783155, 152132.170449405) |
349200660 | NOWY SĄCZ | 9.2180822 | 1014.784 | c(621910.953224735, 196894.713701549) |
349210670 | KROSNO | 8.7728767 | 1017.553 | c(699590.639267207, 208053.055658539) |
349220690 | LESKO | 8.3227397 | 1017.737 | c(742035.646054704, 183036.071495296) |
350200570 | KIELCE-SUKÓW | 8.4849315 | 1016.895 | c(619185.547042676, 328424.723335625) |
plot(dem)
plot(st_geometry(meteo), axes = TRUE)
meteo = stacje_meteo %>% inner_join(dane_meteo_wyb, by = c("KOD_SZS" = "Kod.stacji")) %>% st_transform(4326)
tm_shape(dem) + tm_raster() + tm_shape(meteo) + tm_symbols(col = "tavg", palette = "RdBu")
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
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)
wlkp_pn_lc = crop(lc, wlkp_pn)plot(wlkp_pn_lc)
wlkp_pn_lc2 = mask(wlkp_pn_lc, wlkp_pn)plot(wlkp_pn_lc2)
BIAŁYSTOK
, KRAKÓW-BALICE
, oraz POZNAŃ
. (Porada: użyj do tego, między innymi funkcji right_join()
).library(landscapemetrics)library(landscapetools)
show_landscape(wlkp_pn_lc, discrete = TRUE)
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
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
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
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
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
mapa_p_para = get_lsm(wlkp_pn_lc, what = "lsm_p_para")
## Warning: only using 'what' argument
show_landscape(mapa_p_para[[1]][[1]])
show_lsm(wlkp_pn_lc, what = "lsm_p_para")
punkt = st_sf(st_sfc(geom = st_point(c(347500, 493000))))plot(wlkp_pn_lc)plot(st_geometry(punkt), add = TRUE, cex = 2)
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.
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.polska_lc.tif
do zasięgu (obwiedni) Białowieskiego Parku Narodowego.unique()
)?TA
i CA
).
Jakie są trzy kategorie o największej powierzchni?NP
).Dodatkowe materiały na https://nowosad.github.io/elp/ergosum.html#resources:
#rspatial
, #geocompr
Twitter: jakub_nowosad
Email: nowosad.jakub@gmail.com
Keyboard shortcuts
↑, ←, Pg Up, k | Go to previous slide |
↓, →, Pg Dn, Space, j | Go to next slide |
Home | Go to first slide |
End | Go to last slide |
Number + Return | Go to specific slide |
b / m / f | Toggle blackout / mirrored / fullscreen mode |
c | Clone slideshow |
p | Toggle presenter mode |
t | Restart the presentation timer |
?, h | Toggle this help |
Esc | Back to slideshow |