R jako GIS

Podpůrný text pro předmět: GIS pro biologické aplikace
Autor: Matěj Man
Aktualizace: 01. 12. 2024
Odkaz na R skript: Kopírovat a vložit do R
Odkaz na data: Stáhnout data zip Odkaz na github repozirář: Repozitář

Další zdroje:

Knihovny

# list.of.packages <- c("sf", "terra", "mapview", "randomcoloR", "leaflet", "RColorBrewer")
# install.packages(list.of.packages)

library(RColorBrewer)
library(sf)
library(terra)
library(mapview)
library(randomcoloR)
library(leaflet)

Načítání vektorových GIS dat do R

## Načítání vektorových dat shp
# Nastavte, kde máte u sebe na PC data
# Pozor, musíte zdvojit nebo otočit lomítka
cesta <- "d:/Git/gisproba/data/"
# Zkonstruuje cestu, kde leží vrstva Brdy
data.path <- paste0(cesta, "CHKO_Brdy.shp")
# Načte .shp do R
brdy <- st_read(data.path, stringsAsFactors = FALSE, quiet = TRUE)

# Obdobně načteme třeba hranici ČR
data.path <- paste0(cesta, "hrcr1_wgs.shp")
hrcr <- st_read(data.path, stringsAsFactors = FALSE, quiet = TRUE)

# Prostý obrázek, bez interaktivity
plot(st_geometry(hrcr)) 
# Parametr add přidá další vrstvu do existujícího obrázku
plot(st_geometry(brdy), col = "red", add = TRUE) 

# Nastaví pracovní adresář (Working Directory)
# Dále nemusíme data z WD volat celou cestou, stačí názvy 
setwd(cesta) 
## Načítání vektorových dat z tabulky
# Načíst prostou tabulku
chmu <- read.table("staniceCHMUtablecoma.csv", header = TRUE, sep = ",")
# Převést tabulky na prostorová data
chmu_sf <- st_as_sf(chmu, coords = c("Xcoo", "Ycoo"), crs = 4326) 

## Vizualizovat prostorová data interaktivně
# Pokročilejší balíček leaflet

# Definujeme paletu o 17 barvách podle typu stanice
distCol <- colorFactor(distinctColorPalette(17), chmu_sf$Station.ty)

# Definujeme okno s mapou 
leaflet(chmu_sf) %>% 
  addTiles() %>%  
  addCircleMarkers(popup = chmu_sf$Name, color = ~distCol(Station.ty),
                   stroke = FALSE, fillOpacity = 0.8, radius = 4) %>%
  addLegend(pal = distCol, values = ~Station.ty)

Projekce a CRS

# st_crs(brdy) # Jaký CRS má vrstva brdy?
# st_crs(hrcr) # Jaký CRS má vrstva hranice ČR?

## Vektory
# Transformace podle EPSG
brdy_32633 <- st_transform(brdy, 32633)
# st_crs(brdy_32633) # Jaký CRS má vrstva brdy_32633?

# Transformace podle proj4 string
brdy_5514 <- st_transform(brdy, "+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56")
# st_crs(brdy_5514) # Jaký CRS má vrstva brdy_5514?

# Transformace podle jiné existující vrstvy
brdy_krovak <- st_transform(brdy, st_crs(brdy_5514))
# st_crs(brdy_krovak) # Jaký CRS má vrstva brdy_krovak?

# Kde jsou ty Brdy? Nevidíme je, protože to má jiný CRS
plot(st_geometry(hrcr), main = "Kde jsou ty Brdy?") 
plot(st_geometry(brdy_32633), add = TRUE) 

# Musíme tedy transformovat jedno nebo druhé
hrcr_32633 <- st_transform(hrcr, 32633)
plot(st_geometry(hrcr_32633), main = "No jo, už máme správný CRS")
plot(st_geometry(brdy_32633), add = TRUE, col = "red")

Načítání rastrových GIS dat do R

# Načítání jednoduchou funkcí "rast" z balíčku terra
DEM <- rast(paste0(cesta, "DEM_Jested_100m.tif"))
# Kontrolní obrázek
# Data tu jsou
plot(DEM)

# Ale leží na správném místě? 
# Víme, že to je Křovák, chceme přesnější klíč pro transformaci do WGS
# http://freegis.fsv.cvut.cz/gwiki/S-JTSK

crs(DEM) <- "+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56"

# I rastr je možné vizualizovat interaktivně
# Definujeme barvy
pal <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), values(DEM),
                    na.color = "transparent")

leaflet() %>% 
   addTiles() %>%  
   addRasterImage(DEM, colors = pal, opacity = 0.8) %>%
   addLegend(pal = pal, values = values(DEM),
             title = "Elevation [m]")

Příklad vektorové analýzy

# Plocha polygonů
st_area(brdy)
## 343898692 [m^2]
# Buffer 5 km
brdy_5000 <- st_buffer(brdy_32633, 5000)
plot(st_geometry(brdy_5000), main = "Buffer 5 km")
plot(st_geometry(brdy_32633), add = TRUE, col = "red")

# Centroidy
brdy_cen <- st_centroid(brdy_32633) 
## Warning: st_centroid assumes attributes are constant over geometries
plot(st_geometry(brdy_32633), main = "Centroid", graticule = TRUE, axes = TRUE)
plot(st_geometry(brdy_cen), add = TRUE, col = "blue", pch = 19)

# Průnik
chmu_brdy <- st_intersection(brdy, chmu_sf)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
plot(st_geometry(brdy), main = "Intersection - CHMU stanice v Brdech", graticule = TRUE, axes = TRUE)
plot(st_geometry(chmu_brdy), add = TRUE, col = "red", pch = 19)

Hromadné zpracování dat

Aneb co by nám v QGIS trvalo klikat hodiny, můžeme v R naprogramovat během minut

# Načteme všechny cesty k SHP vrstvám v daném adresáři
parky.files <- list.files(path = cesta, pattern = "*.shp$", full.names = TRUE)
parky.files <- parky.files[-1]

## Když chceme všechny soubory do jednoho objektu
# Načte první SHP soubor
parky.init <- st_read(parky.files[1], quiet = TRUE)

# Připojí všechny další SHP soubory do jednoho objektu  
for (i in 2:length(parky.files)) {
  parky.next <- st_read(parky.files[i], quiet = TRUE)
  parky.init <- rbind(parky.init, parky.next)
}  

head(parky.init)
## Simple feature collection with 6 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 13.6537 ymin: 48.81219 xmax: 18.74347 ymax: 50.67169
## Geodetic CRS:  WGS 84
##   OBJECTID                       geometry
## 1     3900 POLYGON ((18.53572 49.66022...
## 2     3894 POLYGON ((17.88455 49.15734...
## 3     3873 POLYGON ((14.82518 49.66756...
## 4     3877 POLYGON ((14.19356 49.00166...
## 5     3904 POLYGON ((13.93738 49.81539...
## 6     3888 POLYGON ((16.10584 50.662, ...
plot(st_geometry(hrcr))
plot(st_geometry(parky.init), add = TRUE, col = "green")

# Připojit metadata (GIS join)
meta.path <- paste0(cesta, "metadata_parky.csv")
metadata <- read.table(meta.path, sep = ";", header = TRUE, stringsAsFactors = TRUE)
metadata$NAZEV <- iconv(metadata$NAZEV, from = "windows-1250", to = "UTF-8")
parky <- merge(parky.init, metadata)

head(parky)
## Simple feature collection with 6 features and 10 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: 13.61426 ymin: 48.80132 xmax: 15.00087 ymax: 50.67842
## Geodetic CRS:  WGS 84
##   OBJECTID KOD  KAT                    NAZEV      ROZL  ZMENA_G  ZMENA_T
## 1     3873  21 CHKO                   Blaník  4029.195 20030101 20170901
## 2     3874  22 CHKO               Český kras 13225.701 20120824 20170901
## 3     3875  23 CHKO Kokořínsko - Máchův kraj 41037.143 20150223 20170901
## 4     3876  24 CHKO             Křivoklátsko 62497.146 20131107 20170901
## 5     3877  31 CHKO              Blanský les 21968.998 20091120 20170901
## 6     3878  32 CHKO                Třeboňsko 68744.620 20091218 20170901
##   PREKRYV SHAPEAREA  SHAPELEN                       geometry
## 1       0  40291954  31958.58 POLYGON ((14.82518 49.66756...
## 2       0 132257006  83895.23 POLYGON ((14.3244 50.00705,...
## 3       0 410371433 198507.51 MULTIPOLYGON (((14.4139 50....
## 4       0 624971460 138706.87 POLYGON ((13.81561 50.14976...
## 5       0 219689976  79790.35 POLYGON ((14.19356 49.00166...
## 6       0 687446196 143605.35 POLYGON ((14.71406 49.18412...
plot(st_geometry(hrcr))
plot(st_geometry(parky), add = TRUE, col = "red")

# Interaktivně
distCol <- colorFactor(distinctColorPalette(32), parky.init$OBJECTID)
leaflet(parky) %>% 
  addTiles() %>%  
  addPolygons(color = ~distCol(OBJECTID),
              stroke = FALSE, fillOpacity = 0.8,
              label = paste0(parky$KAT, "_", parky$NAZEV))

Rastrové analýzy

Velice dobrá knihovna pro rastrové analýzy v R je terra, která umožňuje efektivní práci s rastrovými daty a poskytuje mnoho funkcí pro terénní analýzy.

# Rastrové analýzy pomocí balíčku terra

# Načteme DEM
dem <- rast(paste0(cesta, "DEM_Jested_100m.tif"))
crs(dem) <- "+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56"

# Rastrové analýzy pomocí balíčku terra

# Načtení DEM
dem <- rast(paste0(cesta, "DEM_Jested_100m.tif"))
crs(dem) <- "+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +units=m +no_defs"

## Stínovaný reliéf
# Vypočítáme sklon a orientaci svahu v radiánech
slope_rad <- terrain(dem, 'slope', unit = 'radians')
aspect_rad <- terrain(dem, 'aspect', unit = 'radians')
# Vypočítáme stínovaný reliéf
hill <- shade(slope_rad, aspect_rad, angle = 30, direction = 315)

# Definujeme kartografickou paletu pro výšku
elevation_palette <- terrain.colors(25, alpha = 0.35)

# Vykreslíme výsledky
plot(hill, col = grey(0:100/100), legend = FALSE, main = "Stínovaný reliéf")
plot(dem, col = elevation_palette, add = TRUE)

## Sklon svahu
# Vypočítáme sklon v stupních
slope_deg <- terrain(dem, 'slope', unit = 'degrees')

# Definujeme barevnou paletu pro sklon

slope_palette <- brewer.pal(9, "YlOrRd")  # Žlutá až červená
slope_colors <- colorRampPalette(slope_palette)(100)

# Vykreslíme sklon svahu
plot(hill, col = grey(0:100/100), legend = FALSE, main = "Sklon svahu")
plot(slope_deg, col = slope_colors, add = TRUE)

## Terrain Ruggedness Index (TRI)
# Vypočítáme TRI
tri_raster <- terrain(dem, 'TRI')

# Definujeme barevnou paletu pro TRI
tri_palette <- brewer.pal(9, "PuBuGn")  # Fialová až modrozelená
tri_colors <- colorRampPalette(tri_palette)(100)

# Vykreslíme TRI
plot(hill, col = grey(0:100/100), legend = FALSE, main = "Terrain Ruggedness Index")
plot(tri_raster, col = tri_colors, add = TRUE)