#####################################################################
#### CARGAR MAPAS DESDE DISTINTOS GEOSERVICIOS ######################
##################################################################

rm(list=ls())

library(sf)
library(mapview)

setwd(choose.dir(getwd(), "Seleccione Dirección de Trabajo"))

################### Geoservicio IDECOR#

## Visualizar las capas disponibles en WFS
idecor <- "WFS:http://idecor-ws.mapascordoba.gob.ar/geoserver/idecor/wms?request=GetCapabilities"

capas_idecor <- st_layers(idecor)

## Cargar capa de parcelas
parcelas <- st_read(idecor,"idecor:parcelas_cba")

names(parcelas)
parcelas = parcelas[,c("Superficie_Mejoras", "Superficie_Tierra_Urbana",
                       "vut_vigente", "Cantidad_Cuentas")]

summary(parcelas$geom)
parcelas = st_cast(parcelas, "GEOMETRYCOLLECTION") %>%    # Cambiar geometria a poligono
  st_collection_extract("POLYGON")

summary(parcelas)

## Cargar capa de escuelas mediante archvo .shp
escuelaspriv <- st_read("establecimientos_privado.shp")
escuelasest  <- st_read("establecimientos_estatales.shp")

names(escuelaspriv)
table(escuelasest$departamen)
table(escuelaspriv$departamen)

escuelaspriv <- subset(escuelaspriv, departamen=="CAPITAL")
escuelasest<- subset(escuelasest, departamen=="CAPITAL")

escuelas<- rbind(escuelasest[,c("sector", "geometry")],escuelaspriv[,c("sector","geometry")])
escuelas <- st_transform(escuelas, 22174)
table(escuelas$sector)

################## GEOSERVICIOS INDEC#

## Visualizar las capas disponibles en WFS
indec <- "WFS:http://geoservicios.indec.gov.ar/geoserver/ows?service=wfs&version=1.3.0&request=GetCapabilities"

capas_indec <- st_layers(indec)

## Cargar radios censales
radios = st_read(indec, "geocenso2010:radios_codigo")

names(radios)
radios = subset(radios, radios$codpcia=="14" & radios$coddpto=="014") # Unicamente Ciudad de Cordoba

summary(radios$geom)
radios = st_cast(radios, "GEOMETRYCOLLECTION") %>%   # Cambiar geometria a poligono
  st_collection_extract("POLYGON")

radios = st_transform(radios, 22174)   # Cambiar sistema de coordenadas

mapview(radios)

## Cargar NBI y otras variables para caracterizar los radios censales
nbi <- st_read(indec,"geocenso2010:nbi_radio")
nbi <- st_drop_geometry(nbi)
names(nbi)
nbi<- nbi[,c("link", "personas_con_nbi_porc", "total_pob", "hogares_con_nbi_porc", "total_hog")]

calidad_construccion <- st_read(indec, "geocenso2010:incalcons_radios")
calidad_construccion <- st_drop_geometry(calidad_construccion)
names(calidad_construccion)
calidad_construccion<- calidad_construccion[,c("link", "satisfactoria_porcentaje", "basico_porcentaje",
                                               "insuficiente_porcentaje")]

names(calidad_construccion)[names(calidad_construccion)=="satisfactoria_porcentaje"] <- "cons_sat"
names(calidad_construccion)[names(calidad_construccion)=="basico_porcentaje"] <- "cons_bas"
names(calidad_construccion)[names(calidad_construccion)=="insuficiente_porcentaje"] <- "cons_insf"


calidad_servicios <- st_read(indec, "geocenso2010:incalserv_radio")
calidad_servicios <- st_drop_geometry(calidad_servicios)
names(calidad_servicios)
calidad_servicios<- calidad_servicios[,c("link", "satisfactoria_porcentaje", "basica_porcentaje",
                                         "insuficiente_porcentaje")]

names(calidad_servicios)[names(calidad_servicios)=="satisfactoria_porcentaje"] <- "serv_sat"
names(calidad_servicios)[names(calidad_servicios)=="basica_porcentaje"] <- "serv_bas"
names(calidad_servicios)[names(calidad_servicios)=="insuficiente_porcentaje"] <- "serv_insf"



calidad_materiales <- st_read(indec, "geocenso2010:inmat_radio")
calidad_materiales <- st_drop_geometry(calidad_materiales)
names(calidad_materiales)
calidad_materiales<- calidad_materiales[,c("link", "calidad_1_porcentaje", "calidad_2_porcentaje",
                                           "calidad_3_porcentaje", "calidad_4_porcentaje")]

names(calidad_materiales)[names(calidad_materiales)=="calidad_1_porcentaje"] <- "mat_1"
names(calidad_materiales)[names(calidad_materiales)=="calidad_2_porcentaje"] <- "mat_2"
names(calidad_materiales)[names(calidad_materiales)=="calidad_3_porcentaje"] <- "mat_3"
names(calidad_materiales)[names(calidad_materiales)=="calidad_4_porcentaje"] <- "mat_4"


actividad <- st_read(indec, "geocenso2010:actividad_radio")
actividad <- st_drop_geometry(actividad)
actividad<- actividad[,c("link", "ocupada", "desocupada", "inactiva", "pea", "población_14_años_y_más",
                         "tasa_actividad", "tasa_empleo", "tasa_desocupacion")]
names(actividad)

##### UNION DE BASES - UNIFICACION DE LA INFORMACION #

## Escuelas con radios censales
aux <- st_join(radios["link"], escuelas, join=st_intersects)
names(aux)
summary(aux$sector)

aux$escuelaspriv <- ifelse(aux$sector=="Privado", 1, 0)
aux$escuelasest <- ifelse(aux$sector=="Estatal", 1, 0)

table(aux$escuelaspriv)

# agrupar las escuelas por radios (variable link)
library(tidyverse)
radios_escuela = aux %>%
  group_by(link) %>%
  summarise(priv = sum(escuelaspriv, na.rm=TRUE),
            est = sum(escuelasest, na.rm=TRUE))

table(radios_escuela$priv) 
table(radios_escuela$est) 


class(radios_escuela)
radios_escuela <- st_drop_geometry(radios_escuela)


## Parcelas con radios censales
parcelas_link = st_join(parcelas, radios["link"], join = st_intersects)
names(parcelas_link)
summary(parcelas_link)
parcelas_link = subset(parcelas_link, is.na(link)==FALSE)
parcelas_link = st_drop_geometry(parcelas_link)


# agrupar las parcelas por radios (variable link)
library(tidyverse)
radios_parcelas = parcelas_link %>%
  group_by(link) %>%
  summarise(vut = mean(vut_vigente, na.rm=TRUE),
            edif = mean(Superficie_Mejoras, na.rm=TRUE),
            sup = mean(Superficie_Tierra_Urbana, na.rm=TRUE),
            ctas = sum(Cantidad_Cuentas, na.rm=TRUE))
summary(radios_parcelas)

radios_parcelas <- na.omit(radios_parcelas)

## Union de parcelas y escuelas por radio censal (variable = link)
aux2 <-left_join(radios_parcelas, radios_escuela, by="link")
names(aux2)

## Union de NBI (variables = link)
aux3 <- left_join(aux2, nbi, by="link")
names(aux3)

##
aux4 <- left_join(aux3, calidad_construccion, by="link")
names(aux4)

aux5 <- left_join(aux4, calidad_servicios, by="link")
names(aux5)

aux6 <- left_join(aux5, calidad_materiales, by="link")
names(aux6)

aux7 <- left_join(aux6, actividad, by="link")
names(aux7)
summary(aux7)

names(radios)
base_final <-left_join(radios[,c("link", "geom")], aux7, by="link")
names(base_final)
summary(base_final)
base_final <- na.omit(base_final)


# Se agregan las coordenadas como variables

aux8 <-st_centroid(base_final)

library(tidyverse)
coords <- do.call(rbind, st_geometry(aux8)) %>%
  as_tibble() %>% setNames(c("x","y"))

base_final$x <- coords$x
base_final$y <- coords$y

# Eliminación de outlier

base_final <- subset(base_final, link != "140140515")

# Mapa de la base final

mapview::mapview(base_final)
names(base_final)

# Guardar base final

st_write(base_final, "base_final.gpkg", delete_dsn = T, delete_layer = T)


#####################################################################
#### CLUSTERIZACION - Tecnica FUZZY C-MEANs #########################
##################################################################

rm(list=ls())

setwd(choose.dir(getwd(), "Seleccione Dirección de Trabajo"))

library(sf)
library(clValid)
library(cluster)
library(factoextra)
library(devtools)
library(corrplot)
library(e1071)
library(caret)
library(dplyr)

base_final <- st_read("base_final.gpkg")

## Para clusterizar es necesario que todas las variables sean numéricas. Entonces:

# Se elimina la geometría

datos <-st_drop_geometry(base_final)

# Se eligen las variables (numericas) para realizar la clusterizacion

names(datos)

datos <- datos[c("vut" , "edif" , "sup" , "ctas" , "priv" , "est" , 
                 "personas_con_nbi_porc" , "hogares_con_nbi_porc" , 
                 "cons_sat" , "cons_bas" , "cons_insf" , "serv_sat" , 
                 "serv_bas" , "serv_insf" , "mat_1" , "mat_2" , "mat_3" ,
                 "mat_4" , "x" , "y" , "ocupada" , "desocupada" , "inactiva" ,
                 "pea" , "población_14_años_y_más" , "tasa_actividad" , 
                 "tasa_empleo" , "tasa_desocupacion")]


# Se estandariza las variables, por alguna tecnica de centrado

pre_proceso <- preProcess(datos, method = c("center", "scale"))

# Se calculan los datos estandarizado - quedan todas las variables en la misma escala

datos_est <- predict(pre_proceso, datos) 

## Definir el numero de zonas 

# Elbow Method

set.seed(7)
wss <- sapply(1:15,function(k){kmeans(datos_est, k, nstart=50,iter.max = 15 )$tot.withinss}) 
plot(1:15, wss, type="b", pch = 19, frame = FALSE, xlab="Number of clusters K", ylab="Total within-clusters sum of squares")

# Coeficiente de Particion, Entropia de Particion, Indice XieBeni

cant_zonas<-function(grupo) {
  
  MC_2 <- cmeans(datos_est,grupo,100,method="cmeans",m=1.1)
  
  I2CM <- fclustIndex(MC_2,datos_est, index=c("xie.beni", "fukuyama.sugeno",
                                              "partition.coefficient", 
                                              "partition.entropy", 
                                              "proportion.exponent" ))
  Indices0 <- cbind(I2CM)
  
  XieBeni <-Indices0[1,]
  FukSug <-Indices0[2,]
  CoefPart_1 <-Indices0[3,]
  CoefPart <- 1/CoefPart_1
  EntrPart <-Indices0[4,]
  ExpProp <-Indices0[5,]
  
  Indices <- as.data.frame(rbind(XieBeni,CoefPart,EntrPart))
  Indices
  
  return(Indices)
  
}

tabla_cant_zonas <- do.call("cbind",lapply (4:8,cant_zonas))
colnames(tabla_cant_zonas) <- c("4","5","6", "7", "8")
tabla_cant_zonas

## CLUSTERIZACION ##

# Se eligen 4 zonas para clusterizar - se aplica fuzzy c means

grupo = 4

set.seed (7)
zona <- cmeans(datos_est,grupo,100,method="cmeans",m=1.1)

radios_zona <- base_final[,c("link","vut" , "edif" , "sup" , "ctas" , "priv" , "est" , 
                             "personas_con_nbi_porc" , "hogares_con_nbi_porc" , 
                             "cons_sat" , "cons_bas" , "cons_insf" , "serv_sat" , 
                             "serv_bas" , "serv_insf" , "mat_1" , "mat_2" , "mat_3" ,
                             "mat_4" , "x" , "y" , "ocupada" , "desocupada" , "inactiva" ,
                             "pea" , "población_14_años_y_más" , "tasa_actividad" , 
                             "tasa_empleo" , "tasa_desocupacion")]

radios_zona$cluster <- zona$cluster

radios_zona$zona <- case_when(radios_zona$cluster==1 ~ "A",
                              radios_zona$cluster==4 ~ "B",
                              radios_zona$cluster==3 ~ "C",
                              radios_zona$cluster==2 ~ "D")

table(radios_zona$zona)

radios_zona$cluster <- NULL


library(mapview)
library(RColorBrewer)

# Definir paleta de colores 

col <- c("#d7191c", "#fdae61", "#ffffbf", "#a6d96a")

# Mapa de las zonas

mapview::mapview(radios_zona,  zcol="zona", col.regions = col, gl =TRUE,
                 alpha.region = 1 , lwd = 1, alpha = 0.3) 

# Guardar la base

getwd()
st_write(radios_zona, "radios_zona.gpkg", delete_dsn = T, delete_layer = T)

#####################################################################
#### ANALISIS DE COMPONENTES PRINCIPALES ############################
##################################################################

library(factoextra)
library(sf)

rm(list=ls())

setwd(choose.dir(getwd(), "Seleccione Dirección de Trabajo"))

radios_zona <- st_read("radios_zona.gpkg")

# Se elimina la geometria 

acp <- st_drop_geometry(radios_zona) 

# Se genera una semilla para que siempre surja el mismo valor

set.seed (7) 

# observo nombre de las variables

names(acp)

# Los componentes principales se encuentran escalados

res.pca <- prcomp(acp[,c(-1,-20,-21,-30)], scale = TRUE) # Se hace el análisis de CP y se los escala
eig.val <- get_eigenvalue(res.pca) # Se calculan los valores propios - Landa - que acopaña a cada CP
round(eig.val, digits = 2) # se los redondea a dos decimales

# Se obtienen los valores para las variables 

res.var <- get_pca_var(res.pca) # Calcula las componentes principales
round(res.var$coord, digits = 2)  # Coordinates
round(res.var$contrib, digits = 2) # Contributions to the PCs
round(res.var$cos2, digits = 2)    # Quality of representation 

round(res.var$coord[,1], digits = 2) 
round(res.var$coord[,2], digits = 2)

fviz_eig(res.pca, ylab= "% CP", xlab= "Comp. Principales", main = "Componentes Principales", font.tickslab = c(12, "bold", "black"), font.title= 20,font.y=15, font.x=15)
fviz_contrib(res.pca, choice = "var",  axes = 1, fill="#06623b", top = 10, font.tickslab = c(12, "bold", "black"), font.title= 20,font.y=15, title=" Contribución CP 1") 
fviz_contrib(res.pca, choice = "var", axes = 2,  fill="#6f0000", top = 10, font.tickslab = c(14, "bold", "black"), font.title= 20,font.y=15, title=" Contribución CP 2") 
fviz_contrib(res.pca, choice = "var", axes = 3,  fill="#00263b", top = 10,font.tickslab = c(14, "bold", "black"), font.title= 20,font.y=15,  title=" Contribución CP 3") 

# Graficar las variables y las CP


col<-c("#000000") # color hunt -https://colorhunt.co/ -

fviz_pca_var(res.pca,
             col.var = "contrib", # Color by contributions to the PC
             gradient.cols = col,
             axes=c(1, 2),
             title="Comp. Princ. Variables Economicos",
             repel = TRUE    # Avoid text overlapping
)

# Se obtienen los valores para las observaciones - individuos-

res.ind <- get_pca_ind(res.pca)
res.ind$coord          # Coordinates
res.ind$contrib        # Contributions to the PCs
res.ind$cos2           # Quality of representation 

groups <- as.factor(radios_zona$zona)

col <- c("#d7191c", "#a6d96a", "#ffffbf", "#fdae61")


fviz_pca_ind(res.pca,
             col.ind = groups, # color by groups
             palette = col,
             addEllipses = TRUE,
             legend.title = "Grupos",
             axes=c(1, 2),
             geom = c("point"),
             title="CP observaciones",
             alpha=1 ) # Concentration ellipses 

## Observar tanto varables como observaciones en las CP

fviz_pca_biplot(res.pca,
                col.ind = groups, # color by groups
                palette = col,
                col.var = "#000000",
                gradient.cols = "fff3af",
                addEllipses = TRUE,
                legend.title = "Grupos",
                axes=c(1, 2),
                geom = c("point"),
                jitter = list(what = "label", width = NULL, height = NULL),
                title="BI - Plot variables - Individuos",
                alpha=1 )



names(radios_zona)

col1 <- c("#d7191c", "#a6d96a", "#ffffbf", "#fdae61")

mapview::mapview(radios_zona,  zcol="zona", col.regions = col1, gl =TRUE,
                 alpha.region = 1 , lwd = 1, alpha = 0.3)