LibrerÃas requeridas
library(tidyverse)## -- Attaching packages ----------------------------------------------------------------------------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.2.1 v purrr 0.3.3
## v tibble 2.1.3 v dplyr 0.8.3
## v tidyr 1.0.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.4.0## -- Conflicts -------------------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(leaflet)
library(sf)## Warning: package 'sf' was built under R version 3.6.3## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3library(sp)
library(rgdal)## Warning: package 'rgdal' was built under R version 3.6.2## rgdal: version: 1.4-8, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/jose.ayala/Documents/R/win-library/3.6/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/jose.ayala/Documents/R/win-library/3.6/rgdal/proj
## Linking to sp version: 1.4-0Lectura de los archivos
df_trafos <- read.csv(file = "transformadores-pcb.csv", stringsAsFactors = FALSE)
comunas_archivo <- readOGR("comunas.shp")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\jose.ayala\Documents\mapatechos\comunas.shp", layer: "comunas"
## with 15 features
## It has 6 fieldsAjuste de coordenadas para el mapa de las comunas
crs_str <- sp::CRS('+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0')
comunas <- sp::spTransform(comunas_archivo,crs_str)Extracción de latitud y longitud para cada transformador
df_trafos$coord <- as.character(df_trafos$WKT)
df_trafos$coord <- df_trafos$coord %>% str_replace_all("POINT \\(", "") %>% str_replace_all("\\)", "") %>% str_split_fixed(" ", 2)
df_trafos$lon <- as.numeric(df_trafos$coord[,1])
df_trafos$lat <- as.numeric(df_trafos$coord[,2])
df_trafos$coord <- NULLVisualización de cada transformador en el mapa (agrupado, esto para que no muestre 7000 puntos juntos)
leaflet() %>%
addTiles()%>%
addMarkers(lng = df_trafos$lon, lat = df_trafos$lat, clusterOptions = markerClusterOptions())Intersercción de coordenadas de cada transformador con los polÃgonos de cada comuna para determinar a cual pertenece
df_traf_sp <- SpatialPointsDataFrame(df_trafos[,c("lon", "lat")], df_trafos[,1:6], proj4string = crs_str)
df_intersect <- over(df_traf_sp, comunas)
df_trafos$comuna <- df_intersect$COMUNAS
df_traf_comuna <- df_trafos %>% group_by(comuna) %>% summarize(transformadores = n()) %>% drop_na()
df_traf_comuna$transformadores <- as.numeric(df_traf_comuna$transformadores)Visualización de transformadores por cada comuna
df_traf_comuna## # A tibble: 15 x 2
## comuna transformadores
## <dbl> <dbl>
## 1 1 1825
## 2 2 542
## 3 3 466
## 4 4 565
## 5 5 273
## 6 6 318
## 7 7 312
## 8 8 250
## 9 9 282
## 10 10 204
## 11 11 232
## 12 12 273
## 13 13 547
## 14 14 711
## 15 15 365Extracción de la superficie de cada comuna
df_area_comuna <- data.frame(comunas$COMUNAS, comunas$AREA)
colnames(df_area_comuna) <- c("comuna", "area")Merge con la el df de transformadores y cálculo de transformadores por km cuadrado
df_traf_comuna <- merge(df_traf_comuna, df_area_comuna, by = "comuna") %>% mutate(area = round(area/1000000,1), trafporkm2 = transformadores/area)Visualización de transformadores por km2 según comuna
df_traf_comuna## comuna transformadores area trafporkm2
## 1 1 1825 17.8 102.52809
## 2 2 542 6.3 86.03175
## 3 3 466 6.4 72.81250
## 4 4 565 21.7 26.03687
## 5 5 273 6.7 40.74627
## 6 6 318 6.9 46.08696
## 7 7 312 12.4 25.16129
## 8 8 250 22.2 11.26126
## 9 9 282 16.5 17.09091
## 10 10 204 12.7 16.06299
## 11 11 232 14.1 16.45390
## 12 12 273 15.6 17.50000
## 13 13 547 14.7 37.21088
## 14 14 711 15.8 45.00000
## 15 15 365 14.3 25.52448Agregamos como dato adicional la cantidad de habitantes según censo 2010 y calculamos transformadores cada 100.000 habitantes
habitantescomunas <- c(205886, 157932, 187537, 218245, 179005, 176076, 220591, 187237, 161797, 166022, 189832, 200116, 231331, 225970, 182574)
df_traf_comuna$habitantes <- habitantescomunas
df_traf_comuna <- df_traf_comuna %>% mutate(trafpor100khab = round(transformadores/(habitantes/100000),2))Visualización de transformadores por 100.000 habitantes según comuna
df_traf_comuna## comuna transformadores area trafporkm2 habitantes trafpor100khab
## 1 1 1825 17.8 102.52809 205886 886.41
## 2 2 542 6.3 86.03175 157932 343.19
## 3 3 466 6.4 72.81250 187537 248.48
## 4 4 565 21.7 26.03687 218245 258.88
## 5 5 273 6.7 40.74627 179005 152.51
## 6 6 318 6.9 46.08696 176076 180.60
## 7 7 312 12.4 25.16129 220591 141.44
## 8 8 250 22.2 11.26126 187237 133.52
## 9 9 282 16.5 17.09091 161797 174.29
## 10 10 204 12.7 16.06299 166022 122.88
## 11 11 232 14.1 16.45390 189832 122.21
## 12 12 273 15.6 17.50000 200116 136.42
## 13 13 547 14.7 37.21088 231331 236.46
## 14 14 711 15.8 45.00000 225970 314.64
## 15 15 365 14.3 25.52448 182574 199.92Volviendo al dataframe de transformadores original queremos obtener más detalles sobre el tipo de aislante y cantidad de PCB por cada comuna
df_trafos <- df_trafos %>% filter(DIELECTRIC == "ACEITE" | DIELECTRIC == "RESINA EPOXI")
df_trafos$PCB <- df_trafos$PCB %>% str_remove_all("<") %>% str_replace_all(",",".") %>% as.numeric() ## Warning in function_list[[k]](value): NAs introduced by coercionPrimero cuento la proporción de transformadores con resina vs con aceite
prop_traf_comuna <- df_trafos %>% group_by(comuna, DIELECTRIC) %>% summarize(transformadores = n()) %>% drop_na() %>% spread(DIELECTRIC, transformadores)
colnames(prop_traf_comuna) <- c("comuna", "aceite", "resina")
prop_traf_comuna <- prop_traf_comuna %>% mutate(
total = aceite + resina,
propaceite = round((aceite/total)*100,1),
propresina = round((resina/total)*100,1)
)
prop_traf_comuna## # A tibble: 15 x 6
## # Groups: comuna [15]
## comuna aceite resina total propaceite propresina
## <dbl> <int> <int> <int> <dbl> <dbl>
## 1 1 1557 268 1825 85.3 14.7
## 2 2 484 46 530 91.3 8.7
## 3 3 417 49 466 89.5 10.5
## 4 4 531 33 564 94.1 5.9
## 5 5 251 22 273 91.9 8.1
## 6 6 287 31 318 90.3 9.7
## 7 7 291 21 312 93.3 6.7
## 8 8 244 6 250 97.6 2.4
## 9 9 271 11 282 96.1 3.9
## 10 10 201 3 204 98.5 1.5
## 11 11 226 6 232 97.4 2.6
## 12 12 257 16 273 94.1 5.9
## 13 13 507 39 546 92.9 7.1
## 14 14 653 57 710 92 8
## 15 15 318 47 365 87.1 12.9En general se considera que los transformadores con aceite son más seguros que los de resina. En base a esto elaboramos el mapa
pal <- colorBin("YlOrRd", domain = prop_traf_comuna$propaceite)
etiquetas <- sprintf(
"<strong>Comuna %s</strong><br/>Porcentaje con aceite %g <br/>Porcentaje con resina %g <br/>Total transformadores %g ",
prop_traf_comuna$comuna, prop_traf_comuna$propaceite, prop_traf_comuna$propresina, prop_traf_comuna$total
) %>% lapply(htmltools::HTML)
leaflet(comunas) %>%
addPolygons(
fillColor = ~pal(prop_traf_comuna$propaceite),
weight = 0.1,
fillOpacity = 1,
label = etiquetas,
highlight = highlightOptions(
weight = 5,
color = "#666",
fillOpacity = 0.7,
bringToFront = FALSE)
)%>% addLegend(pal = pal, values = ~prop_traf_comuna$propaceite, opacity = 0.7, title = "% traf. con aceite",
position = "bottomright")