library(spdep)
library(readxl)
library(rgdal)

setwd('.')

# nacitanie dat
data <- read_excel('data.xlsx')
map <- readOGR(dsn = '.', 'SPH_OBEC')
names(map) <- c('uzkod', 'KOD_LAU1', 'KOD_ORP', 'KOD_OPU', 'KOD_LAU2', 'NAZEV_LAU2')
spatial_tab <- merge(map, data, by = 'uzkod')


# LISA klastrova mapa

# vytvorenie matice vah
coords <- coordinates(map)
IDs <- row.names(as(map, "data.frame"))
nb <- poly2nb(map, queen = TRUE)
listw <- nb2listw(nb, style = 'W')

local <- localmoran(x = spatial_tab$podielNezam_T, listw = listw) # vypocet lokalnej Moranovej statistiky
kvadrant <- vector(mode="numeric",length=nrow(local))
m.nezam <- spatial_tab$podielNezam_T - mean(spatial_tab$podielNezam_T) # centruje sledovanu premennu  
m.local <- local[,1] - mean(local[,1]) # centruje lokalne Moranovo I 
signif <- 0.1 # hladina vyznamnosti

# vytvorenie kvadrantov
kvadrant[m.nezam >0 & m.local>0] <- 4  
kvadrant[m.nezam <0 & m.local<0] <- 1      
kvadrant[m.nezam <0 & m.local>0] <- 2
kvadrant[m.nezam >0 & m.local<0] <- 3
kvadrant[local[,5]>signif] <- 0   
brks <- c(0,1,2,3,4)
colors <- c("lightgrey","blue",rgb(0,0,1,alpha=0.4),rgb(1,0,0,alpha=0.4),"red")
plot(map,border=NA,col=colors[findInterval(kvadrant,brks,all.inside=FALSE)])
legend(x = 'topright',legend=c("nevýznamné","nízka-nízka","nízka-vysoká","vysoká-nízka","vysoká-vysoká"),
       fill=colors,bty="n", cex = 0.7)
