Mapa espacios verdes

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.3

library(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-0

library(rgeos)

## Warning: package 'rgeos' was built under R version 3.6.3

## rgeos version: 0.5-3, (SVN revision 634)
##  GEOS runtime version: 3.8.0-CAPI-1.13.1 
##  Linking to sp version: 1.4-2 
##  Polygon checking: TRUE

library(ggplot2)

Lectura de los archivos

Lectura de los archivos

verde_pub_archivo <- readOGR("espacio-verde-publico.shp", encoding = "UTF-8")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\jose.ayala\Documents\mapatechos\espacio-verde-publico.shp", layer: "espacio-verde-publico"
## with 1736 features
## It has 36 fields
## Integer64 fields read as strings:  SUP_TOTAL

verde_priv_archivo <- readOGR("espacio-verde-privado.shp", encoding = "UTF-8")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\jose.ayala\Documents\mapatechos\espacio-verde-privado.shp", layer: "espacio-verde-privado"
## with 443 features
## It has 19 fields
## Integer64 fields read as strings:  id_ev_priv

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 fields

Ajuste de proyección

crs_str <- sp::CRS('+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0')
comunas <- sp::spTransform(comunas_archivo,crs_str)
verde_pub <- sp::spTransform(verde_pub_archivo,crs_str)
verde_priv <- sp::spTransform(verde_priv_archivo,crs_str)

Agregamos 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)

Centroides para cada comuna

centros_comunas <- data.frame(gCentroid(comunas, byid = TRUE))
centros_comunas$comuna <- paste0("Comuna ", comunas$COMUNAS)

df_instalaciones_pub <- data.frame(verde_pub$patio_de_j, verde_pub$Canil, verde_pub$Posta_aero, verde_pub$REJA)
colnames(df_instalaciones_pub) <- c("patiojuegos", "canil", "postaaero", "reja")
df_instalaciones_pub <- df_instalaciones_pub %>% mutate(
  patiojuegos = ifelse(patiojuegos == "SI", "Patio de juegos ", ""),
  canil = ifelse(canil == "SI", "Canil ", ""),
  postaaero = ifelse(postaaero == "SI", "Posta aeróbica ", ""),
  reja = ifelse(!is.na(reja), "Reja ", ""),
  instalaciones = paste0(patiojuegos, canil, postaaero, reja, sep = " ")
)
instalaciones <- str_remove_all(df_instalaciones_pub$instalaciones, "NA")
instalaciones <- ifelse(instalaciones == " ", "No", instalaciones)

observaciones <- ifelse(is.na(verde_priv$Observacio),"No",verde_priv$Observacio)



etiquetas_pub <- sprintf(
  "<strong>Espacio público %s</strong><br/>Comuna: %g <br/>Area: %g m2<br/>Instalaciones: %s ",
  verde_pub$nombre, verde_pub$COMUNA, verde_pub$area, instalaciones
) %>% lapply(htmltools::HTML)

etiquetas_priv <- sprintf(
  "<strong>Espacio privado %s</strong><br/>Comuna: %g <br/>Area: %g m2 <br/>Observaciones: %s ",
  verde_priv$nombre, verde_priv$COMUNA, verde_priv$area, observaciones
) %>% lapply(htmltools::HTML)

leaflet() %>%
  addProviderTiles(providers$CartoDB.PositronNoLabels)%>%
  addLabelOnlyMarkers(data = centros_comunas,
                      lng = ~x,
                      lat = ~y,
                      label = ~comuna,
                      labelOptions = labelOptions(noHide = TRUE, direction = 'top', textOnly = TRUE))%>%
 addPolygons(data = comunas,
              fillOpacity = 0,
              weight = 0.5) %>%
 addPolygons(data = verde_priv,
              fillColor = "blue",
              weight = 0.1,
              fillOpacity = 1,
              label = etiquetas_priv)%>%
 addPolygons(data = verde_pub,
              fillColor = "green",
              weight = 0.1,
              fillOpacity = 1,
              label = etiquetas_pub)

Cálculo de espacios verdes públicos y privados por comuna

df_priv <- data.frame(verde_priv$COMUNA, verde_priv$area)
df_pub <- data.frame(verde_pub$COMUNA, verde_pub$area)
colnames(df_priv) <- c("comuna", "area")
colnames(df_pub) <- c("comuna", "area")
df_pub <- df_pub %>% filter(comuna != 0)
df_priv$tipo <- rep("privado", length(df_priv$area))
df_pub$tipo <- rep("publico", length(df_pub$area))
df_area_verde <- rbind(df_pub, df_priv) %>% group_by(comuna, tipo) %>% summarize(area = sum(round(area)))

ggplot(df_area_verde, aes(fill = tipo, x = comuna, y = area)) + 
  geom_bar(position = "stack", stat = "identity") +
  scale_y_continuous(labels = scales::comma)+
  scale_x_discrete(limits = 1:15, labels = 1:15)

Cálculo de proporción de espacio verde público por habitante

df_area_verde_pub <- df_area_verde %>% filter(tipo == "publico") %>% mutate(comuna_txt = paste0("Comuna ", comuna))
df_area_verde_pub$habitantes <- habitantescomunas
df_area_verde_pub <- df_area_verde_pub %>% mutate(espverdehab = round(area/habitantes,1))

ggplot(df_area_verde_pub, aes(x = comuna, y = espverdehab)) + 
  geom_bar(stat = "identity", color = "green", fill = "green") +
  scale_y_continuous(labels = scales::comma) +
  scale_x_discrete(limits = 1:15, labels = 1:15)

Volvemos a generar el mapa para mostrar estas diferencias

listacomunas <- comunas$COMUNAS

df_area_verde_pub <- df_area_verde_pub[match(listacomunas, df_area_verde_pub$comuna),]

pal <- colorBin("RdYlGn", domain = df_area_verde_pub$espverdehab, bins = c(0, 1, 2, 4, 6, 8, 10, 25))

etiquetas <- sprintf(
  "<strong>Comuna %s</strong><br/>Habitantes (censo 2010) %g <br/>Espacios verdes públicos %d m2 <br/>Proporción m2 espacio verde por habitante %g ",
  df_area_verde_pub$comuna, df_area_verde_pub$habitantes, df_area_verde_pub$area, df_area_verde_pub$espverdehab
) %>% lapply(htmltools::HTML)

leaflet(comunas) %>%
  addProviderTiles(providers$CartoDB.PositronNoLabels)%>%
  addPolygons(
    fillColor = ~pal(df_area_verde_pub$espverdehab),
    weight = 0.1,
    fillOpacity = 1,
    label = etiquetas,
    highlight = highlightOptions(
      weight = 5,
      color = "#666",
      fillOpacity = 0.7,
      bringToFront = FALSE)
  )%>% addLegend(pal = pal, values = ~df_area_verde_pub$espverdehab, opacity = 1, title = "Espacio verde público por habitante",
  position = "bottomright")