# install.packages(c('data.table','raster','rgeos','rgdal','maptools','ggplot2','gtools','gridExtra'))

setwd('~/ownCloud/Active Docs/vyuka/gis_r')

library(data.table) # jednodussi operace s daty 
library(raster) # prace s rastry
library(rgdal) # prace s prostorovymi daty
library(rgeos) # prace s projekcemi
library(ggplot2)  # efektivni nastroj pro tvorbu grafu
library(RColorBrewer) # barvicky... :)
library(gtools) # uzitecne predprogramovane fce
library(gridExtra) # moznost layoutu pro ggplot

krovak <- '+init=epsg:5513'
wgs84 <- '+init=epsg:4326'
etrs <- '+init=epsg:3035'

povodi <- readOGR('dib_A11_Povodi_vodomer_stanic/A11_Povodi_vodomer_stanic.shp', 'A11_Povodi_vodomer_stanic', verbose = F) # načtení prostorovych dat
povodi@proj4string # zobrazení projekce SHPfilu - pripadne NA je nutno nastavit na správnou projekci

# povodi@proj4string <- CRS(krovak) # nastavení projekce SHPfilu

toky  <- readOGR('dib_A03_Vodni_tok_HU/A03_Vodni_tok_HU.shp', 'A03_Vodni_tok_HU', verbose = F)

stanice <- readOGR('dib_E04_Vodomerne_stanice/E04_Vodomerne_stanice.shp', 'E04_Vodomerne_stanice', verbose = F)

pe_etrs <- spTransform(povodi, CRS(etrs)) # příprava polygonů pro tvorbu processing extent - změna projekce SHPfilu
pe_etrs <- extent(pe_etrs) # processing extent - min, max souradnice
pe_etrs <- 1.15*pe_etrs # zvětšení processing extentu

corine <- raster('g250_06/g250_06.tif') # načtení rastru
projection(corine) # zobrazení projekce rastru, rastry se většinou načitají s projekcí pokud ne ->

# crs(corine) <- etrs # nastavení projekce rastru
corine_s <- crop(corine,pe_etrs) # oriznuti rastru na velikost 'okna'
system.time({ # čas výpočtu
  corine_k <- projectRaster(corine_s, crs = krovak, method = 'ngb') # změna projekce rastru - venujte pozornost pouzite metode
})

pe_wgs <- spTransform(povodi, CRS(wgs84))
pe_wgs <- extent(pe_wgs)
pe_wgs <- 1.15*pe_wgs

dem <- raster(paste(getwd(), '/DMR3_wgs84_cr/DMR3_wgs84_cr.tif', sep = '')) 
projection(dem) 

dem_s <- crop(dem,pe_wgs)
system.time({
  dem_k <- projectRaster(dem_s, crs = krovak)
})

prec <- mixedsort(list.files(paste(getwd(), '/prec_16_tif/', sep = ''), full.names = T, pattern = '.tif')) # načtení a seřazení všech souborů z adresáře s koncovkou *.tif
prec <- stack(prec) # uloženi načtených souborů do proměnné 'tprec'
srazky <- sum(prec) # výpočet ročního úhrnu
projection(srazky)

srazky_s <- crop(srazky,pe_wgs)
system.time({
  srazky_k <- projectRaster(srazky_s, crs = krovak)
})

att <- stanice@data # načtení atributové tabulky do proměnné 'att'
head(att)

vyber <- c('016000') # výběr uzávěrového profilu

stanice_s <- stanice[stanice$DBCN %in% vyber,] # select uzávěrového profilu
stanice_b <- gBuffer(stanice_s, width = 10000, quadsegs = 20) # příklad tvorby bufferu 10km

povodi_s <- gIntersects(povodi, stanice_b, byid=TRUE)[1,] # oříznutí toků podle polygonu selektovaného povod - výsledek jsou hodnoty TRUE/FALSE
povodi_s <- povodi[povodi_s,] # 'export' hodnot TRUE/FALSE do polygonů

toky_s <- toky[povodi_s, ] # elegantnejsi verze intersectu (pouhy subset matice dat)

ggplot() +
  geom_polygon(data = povodi_s, aes(x = long , y = lat, group = id), fill = NA, color = 'black', alpha = .75) +
  geom_path(data = toky_s, aes(x = long, y = lat, group = id), color = 'steelblue4') +
  geom_polygon(data = stanice_b, aes(x = long, y = lat, group = id), color = 'red4', alpha = .25) +
  geom_point(aes(x = coordinates(stanice_s)[1] ,y = coordinates(stanice_s)[2]), color = 'red4') +
  coord_fixed() +
  theme_classic()

centroid <- as.data.frame(coordinates(povodi_s)) # těžiště jednotlivych polygonů
names(centroid) <- c('long','lat')

centroid <- as.data.frame(gCentroid(povodi_s)@coords) # celkove těžiště polygonu
names(centroid) <- c('long','lat')

pe_pov <- extent(povodi_s)
pe_pov <- 1.15*pe_pov

dem_s <- crop(dem_k,pe_pov)
dem_s <- mask(dem_s,povodi_s) # oriznuti rastru podle polygonu

dem_s_p <- data.frame(rasterToPoints(dem_s)) # potřeba převést raster na body pro plot v gg
colnames(dem_s_p) <- c('X','Y','Z')

cols <- terrain.colors(100, alpha=0.1)

ggplot(dem_s_p) + 
  geom_raster(aes(x = X, y = Y, fill = Z)) +
  scale_fill_gradientn(name = 'dem', colours = cols) +
  geom_contour(aes(x = X, y = Y, z = Z), color = 'grey13', lwd = .1) +
  geom_polygon(data = povodi_s, aes(x = long , y = lat, group = id), fill = NA, color = 'black', alpha = .75) +
  geom_path(data = toky_s, aes(x = long, y = lat, group = id), color = 'darkblue') +
  coord_fixed() +
  theme_classic() +
  ylab(NULL) +
  xlab(NULL)  +
  ggtitle('')

teren_sl <- terrain(dem_s, opt = 'slope', unit = 'degrees') # sklonitost, osvit atd. terénu - ??terrain   
# teren_as <- terrain(dem_s, opt = 'aspect', unit = 'degrees') 

teren_sl_p <- data.frame(rasterToPoints(teren_sl))
colnames(teren_sl_p) <- c('X','Y','Z')

brk_sl <- seq(min(teren_sl_p$Z), max(teren_sl_p$Z), .5) 

ggplot(teren_sl_p) + 
  geom_raster(aes(x = X, y = Y, fill = Z)) +
  scale_fill_gradientn(name = 'Slope [°]',
                       colours = cols,
                       breaks = brk_sl,
                       guide = guide_legend(title.position = 'top',
                                            label.position = 'bottom',
                                            label.hjust = .5,
                                            label.vjust = .5,
                                            label.theme = element_text(angle = 45, size = 7.5),
                                            nrow = 1,
                                            keywidth = .75,
                                            keyheight = .5,
                                            byrow = TRUE)) +
  geom_polygon(data = povodi_s, aes(x = long , y = lat, group = id), fill = NA, color = 'black', alpha = .75) +
  geom_path(data = toky_s, aes(x = long, y = lat, group = id), color = 'darkblue') +
  coord_fixed() +
  theme_classic() +
  ylab(NULL) +
  xlab(NULL) +
  ylim(extent(teren_sl)@ymin*1.0025, extent(teren_sl)@ymax) +
  theme(legend.direction = 'horizontal',
        legend.box = 'horizontal',
        legend.position = c(0, 0),
        legend.justification = c(0, 0),
        legend.background = element_rect(fill = 'transparent'),
        plot.title = element_text(hjust = 0.5)) +
  ggtitle('Slope')

srazky_s <- crop(srazky_k,pe_pov)
srazky_s <- mask(srazky_s,povodi_s)

srazky_s_p <- data.frame(rasterToPoints(srazky_s))
colnames(srazky_s_p) <- c('X','Y','Z')

cols <-  colorRampPalette(brewer.pal(9, 'Blues'))(round(max(srazky_s_p$Z) - min(srazky_s_p$Z)))

ggplot(srazky_s_p) + 
  geom_raster(aes(x = X, y = Y, fill = Z)) +
  scale_fill_gradientn(name = '[mm]', colours = cols) +
  geom_contour(aes(x = X, y = Y, z = Z), color = 'grey13', lwd = .1) +
  geom_polygon(data = povodi_s, aes(x = long , y = lat, group = id), fill = NA, color = 'black', alpha = .75) +
  coord_fixed() +
  theme_classic() +
  ylab(NULL) +
  xlab(NULL)  +
  ggtitle('Precipitation')

denst <- data.frame(values(srazky_s)[!is.na(values(srazky_s))])
names(denst) <- 'value'

ggplot() +
  geom_density(data = denst, aes(x = value), fill = cols[median(seq_along(cols))]) +
  theme_classic()

legenda_corine <- fread('g250_06/clc_legend.csv')

corine_s <- crop(corine_k,povodi_s)
corine_s <- mask(corine_s,povodi_s)

clc <- unlist(extract(corine_s, povodi_s))
clc <- table(clc)
clc <- data.frame(clc/sum(clc)*100)
names(clc) = c('GRID_CODE','%')
clc <- merge(clc, legenda_corine[,c(1,5)], by = 'GRID_CODE')

head(clc)

barvicka <- strsplit(as.character(legenda_corine$RGB),'-') # určení barev pro corine podle RGB v legendě
barva_clc <- c()
for (i in 1:length(barvicka)) {
  if(is.finite(as.numeric(barvicka[[i]][1]))) {
    barva_clc[[i]] = rgb(as.numeric(barvicka[[i]][1])/255,as.numeric(barvicka[[i]][2])/255,as.numeric(barvicka[[i]][3])/255)
  }
  else {
    barva_clc[[i]] = rgb(1,1,1)
  }
}

corine_p <- data.frame(rasterToPoints(corine_s))
colnames(corine_p) = c('x','y','value')

legenda_corine$LABEL3 <- sub(', ', '\n', legenda_corine$LABEL3)  

cor_hist <- ggplot(clc) + 
  geom_bar(aes(x = GRID_CODE, y = `%`, fill = GRID_CODE), stat = 'identity') +
  geom_bar(aes(x = GRID_CODE, y = `%`, fill = GRID_CODE), stat = 'identity', colour = 'black', show.legend = FALSE, lwd = .25) +
  scale_fill_manual(name = '',
                    values = barva_clc[as.numeric(legenda_corine$GRID_CODE) %in% clc$GRID_CODE],
                    labels = legenda_corine[as.numeric(legenda_corine$GRID_CODE) %in% clc$GRID_CODE, LABEL3]) +
  coord_fixed(4/length(clc$GRID_CODE)) +
  theme_classic() +
  ylab('%') +
  xlab('Clc code') +
  ggtitle('Corine bar plot')

cor_map <- ggplot() +
  geom_raster(data = corine_p, aes(x, y, fill = factor(value))) +
  geom_polygon(data = povodi_s, aes(x = long , y = lat, group = id), fill = NA, color = 'black', alpha = .75) +
  scale_fill_manual(name = 'Land Cover',
                    values = barva_clc[as.numeric(legenda_corine$GRID_CODE) %in% unique(corine_p$value)],
                    labels = paste0(legenda_corine[as.numeric(legenda_corine$GRID_CODE) %in% unique(corine_p$value), LABEL3],
                                    ' (', legenda_corine[as.numeric(legenda_corine$GRID_CODE) %in% unique(corine_p$value), GRID_CODE], ')')) +
  coord_fixed() +
  theme_classic() +
  ylab(NULL) +
  xlab(NULL) +
  ggtitle('Corine map')

get_legend <- function(myggplot) {
  
  tmp <- ggplot_gtable(ggplot_build(myggplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == 'guide-box')
  legend <- tmp$grobs[[leg]]
  
  return(legend)
}

cor_leg <- get_legend(cor_map)

grid.arrange(cor_map + theme(legend.position = 'none'), 
             cor_leg,
             cor_hist + theme(legend.position = 'none'), 
             layout_matrix = matrix(c(1,2,3,2), nrow = 2, byrow = T))