Mapeando Ingresos en las Ciudades de Uruguay

Rafa · 2019/04/13 · 5 minute read

Todas las ciudades tienen zonas ricas y zonas pobres. En Montevideo, las zonas costeras concentran los hogares de mayores ingresos (Carrasco, Punta Gorda, Malvín, Pocitos, etc.).

En este post armo un mapa para visualizar cuáles son las zonas con mayores ingresos en las ciudades de Uruguay linkeando los mapas que publica el INE con las formas de los segmentos censales y los datos de ingresos de la Encuesta Continua de Hogares:

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/rafa/blog_data/mapas_vectoriales_ine", layer: "ine_seg_11"
## with 4313 features
## It has 13 fields

Las Zonas Censales

La información geográfica puede venir en formato raster o en formato shapefile. Un objeto raster es una grilla de números, y cada celda de la grilla está georeferenciada a una parte de la superficie de la tierra. Los shapefiles son una forma de almacenar datos de formas, que pueden representar entidades geográficas físicas o imaginarias (ríos, países, etc.)

El INE publica shapefiles con las formas de los segmentos censales que usa para muestrear la Encuesta de Hogares. En Uruguay hay nrow(ine_seg_11_sf) segmentos censales.

La librería sf tiene un montón de herramientas para trabajar con datos geográficos. La mejor fuente que conozco para aprender más es Geocomputation with R.

library(sf)

# Cargar mapa
ine_seg_11_sf <- st_read('~/blog_data/mapas_vectoriales_ine/ine_seg_11.shp')
## Reading layer `ine_seg_11' from data source `/Users/rafa/blog_data/mapas_vectoriales_ine/ine_seg_11.shp' using driver `ESRI Shapefile'
## Simple feature collection with 4313 features and 13 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 366582.2 ymin: 6127919 xmax: 858252.1 ymax: 6671738
## epsg (SRID):    NA
## proj4string:    NA
# asignar crs
st_crs(ine_seg_11_sf) <- "+proj=utm +zone=21 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
ine_seg_11_sf <- st_transform(ine_seg_11_sf, crs=4326)
class(ine_seg_11_sf)
## [1] "sf"         "data.frame"
plot(st_geometry(ine_seg_11_sf))

Calculamos las medianas

Para calcular el ingreso en cada segmento censal, usamos los datos de la Encuesta Continua de Hogares. Los microdatos incluyen el segmento censal de cada hogar y una estimación de su ingreso. dplyr hace muy fácil calcular la mediana del ingreso de los hogares encuestados:

Sys.setlocale("LC_ALL", "UTF-8")
## [1] "en_US.UTF-8/UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8"
medianas <- read_tsv("~/blog_data/H_2018_TERCEROS.dat") %>% 
  transmute(id       =as.character(numero),
            depto_id = dpto, 
            depto = nomdpto,
            localidad_id = as.numeric(locagr),
            localidad = nom_locagr,
            seccion_id = as.numeric(secc),
            segmento_id = as.numeric(segm),
            ingreso = YSVL) %>% 
  group_by(depto_id, seccion_id, segmento_id) %>% 
  summarize(mediana = median(ingreso))

Linkear los datos al mapa

Para visualizar todo en un mapa, hay que linkear los datos con el dataframe de los segmentos censales. Agregarle atributos no espaciales a un objeto sf es simple, ya que se comportan como data.frames y soportan las funciones del tidyverse, por lo que podemos usar left_join con depto_id, seccion_id y segmento_id como identificadores:

ine_seg_11_sf <- ine_seg_11_sf %>% 
  left_join(medianas, by=c("DEPTO" = "depto_id",
                           "SECCION" = "seccion_id",
                           "SEGMENTO" = "segmento_id"))

Elegir una localidad

Como queremos mapear los ingresos en cada ciudad, tenemos que encontrar una forma de seleccionar una zona del mapa. Para esto podemos usar filter y seleccionar los polígonos que pertenecen a cierta localidad`:

# Montevideo es la localidad 1020
cod_loc <- "1020"
mapa_loc <- filter(ine_seg_11_sf, CODLOC==cod_loc)
plot(st_geometry(mapa_loc))

Agregar mapa de fondo

Este mapa funciona bien si conocemos la ciudad, pero para zonas geográficas no tan familiares puede ser conveniente agregar referencias como calles importantes o cuerpos de agua.

Para eso podemos usar google maps, que publica una API con mapas basadas en imágenes satelitales.Para consumir la API, es necesario registrarse y darles una tarjeta de crédito 👎 para obtener una clave. Acá. hay más info sobre como trabajar con API keys en R.

library(ggmap)
# Leer la API key de variable de environment
API_KEY <- Sys.getenv("GOOGLE_API")
# Loguearse a la API
register_google(API_KEY)

Para bajar el mapa de Google Maps, necesitamos las coordenadas del mapa que vamos a dibujar. Para obtener estas coordenadas, subseteo el mapa original con el código de la localidad que quiero mapear y obtengo las coordenadas del mapa subseteado.

cod_loc <- "1020"
mapa_loc <- filter(ine_seg_11_sf, CODLOC==cod_loc)

bbox <- st_bbox(mapa_loc) %>% 
   setNames(c("left", "bottom", "right", "top"))

mapa <- get_map(location=bbox)

Para poder dibujar este mapa con tmap, tengo que convertirlo a un Raster con la función ggmap_rast.

# https://gis.stackexchange.com/questions/155334/ggmap-clip-a-map
library(raster)
library(tmap)

ggmap_rast <- function(map) {
  map_bbox <- attr(map, 'bb') 
  .extent <- extent(as.numeric(map_bbox[c(2,4,1,3)]))
  my_map <- raster(.extent, nrow= nrow(map), ncol = ncol(map))
  rgb_cols <- setNames(as.data.frame(t(col2rgb(map))), c('red','green','blue'))
  red <- my_map
  values(red) <- rgb_cols[['red']]
  green <- my_map
  values(green) <- rgb_cols[['green']]
  blue <- my_map
  values(blue) <- rgb_cols[['blue']]
  stack(red,green,blue)
}

mapa_raster <- ggmap_rast(mapa)

Ahora que tengo el Raster, los polígonos de los segmentos censales y la estimación de las medianas, podemos hacer un mapa en capas con tmap:

tm_shape(mapa_raster) + 
  tm_rgb() + 
  tm_shape(mapa_loc) +
  tm_fill(col="mediana", 
          style="quantile",
          title = "Ingreso",
          textNA = "S/D",
          alpha = 0.75) + 
  tm_layout("Montevideo", 
            frame = FALSE,
            legend.position = c("right", "bottom"),
            legend.format = list(
              "text.separator" = "-"
            )) +
  tm_legend(show=FALSE)

Una función para dibujar le mapa de una localidad

tmap nos permite dibujar un mapa con todos los elementos que vinimos desarrollando en distintas capas. También tiene una función que permite armar grillas de mapas. Para usarla, refactoreamos todo el código para armar un mapa en una función que toma el código de la localidad y el título y la llamamos para Salto y Durazno.

mapa_ingresos <- function(cod_loc, nom_loc, legend_pos = c("right","top")) {
  


  mapa_loc <- filter(ine_seg_11_sf, CODLOC==cod_loc)

  bbox <- st_bbox(mapa_loc) %>% 
    setNames(c("left", "bottom", "right", "top"))

  mapa_raster <- get_map(location=bbox) %>% 
    ggmap_rast

  tm_shape(mapa_raster) + 
    tm_rgb() + 
    tm_shape(mapa_loc) +
    tm_fill(col="mediana", 
          style="quantile",
          title = "Ingreso",
          textNA = "S/D",
          alpha = 0.75) + 
  tm_layout(nom_loc, 
            frame = FALSE,
            legend.position = legend_pos,
            legend.format = list(
              "text.separator" = "-"
            ))
}

Para sacar varios mapas a la vez podemos combinar tmap_arrange con esta función:

tmap_arrange(
            mapa_ingresos("15120", "Salto"),
            mapa_ingresos("6220", "Durazno"),
            ncol = 1)