R jako GIS

Podpůrný text pro předmět: GIS pro biologické aplikace
Autor: Matěj Man
Aktualizace: 21. 11. 2022
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:

• Online kniha Geocomputation with R
• Online kniha Spatial Data Science with R
• Velice pěkný návod pro R knihovnu sf
• The Visual Raster Cheat Sheet Rpubs dokument

Knihovny

# list.of.packages <- c("sf","raster","mapview","whitebox","randomcoloR","leaflet","Rcpp")
# install.packages(list.of.packages)
# install.packages("devtools")
# install.packages('raster', repos='https://rspatial.r-universe.dev')
# library(devtools)
# install_github("r-spatial/sf")

library(sf)
library(raster)
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/"
# zkonstuuje cestu kde leží vrstva Brdy
data.path<-paste0(cesta,"CHKO_Brdy.shp")
# načte .shp do R
brdy<-st_read(data.path,stringsAsFactors = F,quiet = T)

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

# 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=T) 

# Nastaví pracovní adresář working dorectory (WD)
# 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 = T, 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)

# definoujeme 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 hrnice č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_5514 ? ------------------------------

# kde jsou ty brdy? no nevidím je protože to má jiný CRS
plot(st_geometry(hrcr),main="Kde jsou ty Brdy?") 
plot(st_geometry(brdy_32633),add=T) 

# 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=T,col="red")

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

# načítání jdnoduchou funkcí "raster"
DEM<-raster(paste0(cesta,"DEM_Jested_100m.tif"))
# kontrolí 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 to wgs
# http://freegis.fsv.cvut.cz/gwiki/S-JTSK

DEM@crs<-CRS("+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(DEM) %>% 
   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]

# délka linií (obvod polygonů)
st_length(brdy)

## 0 [m]

# 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=T,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=T,axes=T)
plot(st_geometry(brdy_cen),add=T,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=T,axes=T)
plot(st_geometry(chmu_brdy),add=T,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 = T)
parky.files<-parky.files[-1]

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

# připojí všechny další shp file do jendoho 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=T,col="green")

# připojit metadat GIS join
meta.path<-paste0(cesta,"metadata_parky.csv")
metadata<-read.table(meta.path,sep=";",header=T,stringsAsFactors = T)
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=T,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 whitebox, což je implementace jinak v Python/Rust napsaných geo algoritmů.

• Homepage whitebox tools
• Manuál pro whitebox
• Homepage r whitebox

## instalace balíčku whitebox pro rastrové analýzy
# install.packages("whitebox", repos="http://R-Forge.R-project.org") # instalace
# whitebox::wbt_init() # inicializace
library(whitebox)

# pozor na diakritiku a mezery v cestě. Pro WBT nesmí být
# cesta k rastrovým datům
dem <-("d:/Git/gisproba/data/DEM_Jested_100m.tif")

## stínovaný reliéf
# cesta k budoucímu výsledku
output<- ("c:/Users/Public/Documents/Hillshade_100m.tif")

# spustí algoitmus pro výpočet stinovaného reliéfu 
wbt_hillshade(dem, output, azimuth = 315, altitude = 30, 
                    zfactor = 1,verbose_mode = FALSE)

## [1] "hillshade - Elapsed Time (excluding I/O): 0.1s"

# načtu výsledek jako rastr
hill<-raster(output)
# načtu původní dem jako rastr
elev<-raster(dem)
# vytisknu oba obrázky pro kontrolu
plot(hill, col=grey(0:100/100), legend=FALSE,main="stínovaný reliéf")
plot(elev, col=rainbow(25, alpha=0.35), add=TRUE)

## Sklon svahu
output<- ("c:/Users/Public/Documents/slope_100m.tif")
wbt_slope(dem,output,zfactor = 1, verbose_mode = FALSE)

## [1] "slope - Elapsed Time (excluding I/O): 0.1s"

slp<-raster(output)
plot(hill, col=grey(0:100/100), legend=FALSE,main="Sklon svahu")
plot(slp, col=heat.colors(25, alpha=0.2), add=TRUE)

## The terrain ruggedness index (TRI)
output<- ("c:/Users/Public/Documents/TRI_100m.tif")
wbt_ruggedness_index(dem,output,zfactor = 1, verbose_mode = FALSE)

## [1] "ruggedness_index - Elapsed Time (excluding I/O): 0.1s"

tri<-raster(output)
plot(hill, col=grey(0:100/100), legend=FALSE,main="terrain ruggedness index")
plot(tri, col=cm.colors(25, alpha=0.3), add=TRUE)