Mapping Accessibility to Health Food in Brazilian Favelas and Urban Comunities

Marcio Antonio Pizzonia Junior

Mapping Accessibility to Health Food in Brazilian Favelas and Urban Comunities

Author

Marcio Antonio Pizzonia Junior

Published

September 23, 2026

Overview

This project investigates the spatial accessibility of healthy-food retail in peripheral areas of São Paulo, Brazil, where the lower availability of establishments selling healthy foods may reflect broader patterns of social inequality in large urban centers (Rocha et al., 2024). Areas characterized by the absence or limited availability of such establishments are commonly described as food deserts, representing not only a spatial constraint on access to food but also a potential barrier to the human right to adequate food and to public health.

The study aims to develop quantitative methods to assess accessibility to healthy-food retail, investigate how accessibility measures vary across different data sources, and provide evidence to support public policies addressing food insecurity and improving the identification of food deserts. Two official databases are compared: the Brazilian Federal Revenue Service (Receita Federal Brasileira, RFB) and the Annual Social Information Report (Relação Anual de Informações Sociais, RAIS).

Food-healthfulness classification follows the Locais-Nova methodology, which uses the Classificação Nacional de Atividades Econômicas (CNAE) codes reported in establishments CNPJ records to classify food-retail establishments according to the healthfulness of the foods they sell. Based on this classification, the AcessoSAN project proposed Group 0 (G0) as a subgroup of the original G1–G2 category, with the purpose of providing greater detail for establishments whose CNAE classification corresponds exclusively to healthy-food retail. Although the analysis focuses primarily on G0, the other groups defined by the Locais-Nova classification will also be mapped to provide a broader spatial characterization of food-retail establishments.

For instructions on how to run this work, see the repository README.

Problem

The AcessoSAN project aims to develop new methods for measuring and analyzing inequalities in access to healthy food, particularly in Brazilian favelas and urban communities. To achieve these objectives, it is necessary to use reliable and up-to-date datasets regarding establishments that sell healthy food.

The Brazilian Federal Revenue Service (RFB) and the Annual Social Information Report (RAIS) provides records of companies operating in Brazil, including food retailers. However, it does not provide geocoded information, which is necessary for spatial analysis.

Data Availability

The processed data are available in csv, rds and gpkg formats through a dedicated repository on the Open Science Framework (OSF). A metadata file is included alongside the validated datasets.

Only authorized personnel can access the processed files. They are protected with RSA 4096-bit encryption (OpenSSL) and a 32-byte password to ensure data security.

If you already have access to the OSF repository and the project keys, click here to access the data. You can also retrieve these files directly from R using the osfr package.

Methods

Source of Data

The data used for this this work comes from:

Brazilian Federal Revenue Service (RFB):
- Legal entities data from the National Register of Legal Entities (CNPJ) database, used for company identification (Brasil et al., n.d.).
Brazilian Ministry of Labor and Employment (MTE):
- Legal food establishments entities data from the Annual Social Information Report (RAIS) related to 2022 (provided by request) (Brasil & Ministério do Trabalho e Emprego, 2022), used for company identification.
- Codes of Brazilian agencies and municipalities (Tabela de Órgãos e Municípios, TOM) (Brasil & Sistema de Vigilância Alimentar e Nutricional, n.d.), used to identify municipalities in the CNPJ database.
Brazilian Institute of Geography and Statistics (IBGE):
- Codes and descriptions from the National Classification of Economic Activities (CNAE) (Instituto Brasileiro de Geografia e Estatística & Comissão Nacional de Classificação, n.d.), used to match establishments to the Locais-Nova groups.
- Brazilian municipalities data via the geobr R package (Pereira & Goncalves, n.d.), used for cross-referencing.
- Open spatial datasets of Brazilian addresses from the National Address Register for Statistical Purposes (CNEFE) via the geocodebr R package (Pereira & Herszenhut, n.d.), used for geocoding the establishments.
- Tract boundary data from the Brazilian 2022 Census and municipal boundary data via the geobr R package (Pereira & Goncalves, n.d.), used for spatial analysis and population density calculations.
- Tract population data from the Brazilian 2022 Census via the censobr R package (Pereira & Goncalves, n.d.), used to calculate population density.
Center for Metropolitan Studies (CEM):
- Lookup table providing the classification linking the establishments CNAEs codes to the Locais-Nova groups (Fernandes & Giannotti, 2026).
Lookup table providing address normalization for geocoding purposes (Vartanian et al., 2026).
MapBiomas
- Land cover and land use data for Brazil using Sentinel 2 images (10m resolution, collection 2) (MapBiomas, n.d.), used to show land use in the plots.
Mapzen Terrain Tile
- Terrain data (Mapzen, n.d.) via the elevatr R package (Hollister et al., 2025), used to calculate travel times and accessibility metrics.
OpenStreetMap (OSM):
- Road and footpaths spatial data (OpenStreetMap Contributors, n.d.) via the Geofabrik project (Geofabrik GmbH, n.d.), used to calculate travel times and accessibility metrics.
Uber Technologies:
- Hexagonal grid data structured using the hexagonal hierarchical geospatial indexing system (H3) (Uber Technologies, n.d.), processed with the h3jsr (O’Brien, 2023) R package, used to create hexagonal grids for spatial analysis and visualization.

Data Munging

The data munging follows the data science workflow outlined by Wickham et al. (2023), as illustrated in Figure 1. All processes were made using the Quarto publishing system and R (R Core Team, n.d.) programming languages, supported by several R packages.

For data manipulation and workflow, priority was given to packages from the tidyverse, rOpenSci and r-spatial frameworks, as well as other packages adhering to the tidy tools manifesto (Wickham, 2023).

Figure 1: Data science workflow created by Wickham, Çetinkaya-Runde, and Grolemund.

Travel Matrix

The travel matrix is computed using the r5r R package (Pereira et al., 2021), which interfaces with the R5 routing engine. This tool enables detailed travel time calculations based on multimodal transportation networks, accounting for road types, traffic conditions, and elevation changes. Travel times are derived from OpenStreetMap (OSM) road and footpath data (OpenStreetMap Contributors, n.d.) and elevation data from the Mapzen Terrain Tile dataset (Mapzen, n.d.).

Code Style

The Tidyverse tidy tools manifesto (Wickham, 2023), code style guide (Wickham, n.d.-a) and design principles (Wickham, n.d.-b) were followed to ensure consistency and enhance readability.

Reproducibility

The pipeline is fully reproducible and can be run again at any time. To ensure consistent results, the renv package (Ushey & Wickham, 2025) was used to manage and restore the R environment. See the README file in the code repository to learn how to run it.

Set the Environment

Load Packages

library(askpass) 
library(brandr)
library(checkmate)
library(cli)
library(cubelyr)
library(dplyr)
library(enderecobr)
library(fs)
library(geobr)
library(geocodebr)
library(ggnewscale)
library(ggplot2)
library(ggspatial)
library(googlesheets4)
library(groomr) # github.com/danielvartan/groomr
library(h3jsr)
library(here)
library(httr2)
library(labelled)
library(lobstr)
library(lockr) # github.com/danielvartan/lockr
library(lubridate)
library(magrittr)
library(nanoparquet)
library(numbersBR)
library(orbis) # github.com/danielvartan/orbis
library(osfr)
library(osmextract)
library(parallel)
library(purrr)
library(pmtiles) # github.com/walkerke/pmtiles
library(quarto)
library(ragg)
library(RColorBrewer)
library(readr)
library(rlang)
library(scales)
library(sf)
library(stars)
library(stringr)
library(tidyr)
library(vroom)
library(xml2)
library(zip)

Define Custom Functions

library(checkmate)
library(sf)

expand_bbox <- function(bbox, ratio = 0.02) {
  assert_class(bbox, "bbox")
  assert_number(ratio, lower = 0.01)

  bbox["xmin"] <- bbox["xmin"] - ((bbox["xmax"] - bbox["xmin"]) * ratio)
  bbox["xmax"] <- bbox["xmax"] + ((bbox["xmax"] - bbox["xmin"]) * ratio)
  bbox["ymin"] <- bbox["ymin"] - ((bbox["ymax"] - bbox["ymin"]) * ratio)
  bbox["ymax"] <- bbox["ymax"] + ((bbox["ymax"] - bbox["ymin"]) * ratio)

  bbox
}

library(checkmate)

get_plot_ratio <- function(plot) {
  assert_class(plot, "ggplot")

  plot_build <- plot |> ggplot_build()

  panel_scales <- plot_build$layout$panel_scales_y[[1]]
  panel_range <- panel_scales$range$range
  panel_height <- diff(panel_range)

  panel_scales <- plot_build$layout$panel_scales_x[[1]]
  panel_range <- panel_scales$range$range
  panel_width <- diff(panel_range)

  panel_height / panel_width
}

map_labs <- function(municipality_code, title, source, note = NULL) {
  assert_string(title)
  assert_string(source)

  municipality_row <-
    municipalities_data |>
    filter(municipality_code == .env$municipality_code)

  labs(
    title = paste0(
      title, " ",
      municipality_row |> pull(municipality), "-",
      municipality_row |> pull(federal_unit)
    ),
    caption = if (is.null(note)) {
      str_glue("Source: {source}.")
    } else {
      str_glue("Source: {source}, {note}.")
    }
  )
}

Set Keys

See the Data Availability section for more information.

osf_auth(Sys.getenv("OSF_PAT")) # askpass())

gs4_auth(
  cache = here(".secrets"),
  email = TRUE,
  use_oob = FALSE
)

public_key <- here(".ssh", "id_rsa.pub")

private_key <- here(".ssh", "id_rsa")

password <- Sys.getenv("ACESSOSAN_PASSWORD") # askpass()

Set and Check Parameters

The pipeline main parameters are defined in the YAML header of its Quarto document.

The year and month variables represent the year and month of the consolidated RFB CNPJ dataset.

The municipality_codes variable identifies the IBGE codes of the municipalities included in the analysis. Click here to see a list of all Brazilian municipalities and their codes.

params |> lobstr::tree()
#> <list>
#> ├─year: 2022
#> ├─month: 1
#> ├─municipality_codes<int [5]>: 2507507, 3106200, 3550308, 4314902, 1501402
#> ├─state_codes<int [5]>: 31, 15, 25, 43, 35
#> └─h3_resolution: 9

The plot_style variable sets the color scheme used in the plots. The default colors are based on the Voyager CARTO color scheme, which is the default color scheme used in the mapgl R package (Walker, 2026) for mapping geospatial data.

plot_style <- list(
  land = "#fbf8f3",
  water_bodies = "#cbe1e4",
  green_areas = "#e8f0dd",
  highways = "#fbdb98"
)

The plot_expansion variable sets the expansion ratio for the bounding box of the municipality shape when plotting the background shapes. The default value of 0.02 means that the bounding box will be expanded by 2% in all directions.

plot_expansion <- 0.02

The group_colors variable sets the sets the color scheme used in the plots.

group_colors <- c(
  "G0" = "#ff087f",  
  "G1-G2" = "#8bc34a",
  "G3" = "#ffc107",  
  "G4" = "#d32f2f"   
)

label_style <- list(
  fcu_label = "#111111" 
)

assert_int(params$year, lower = 1980, upper = 2050)
assert_integerish(params$municipality_codes, lower = 1100015, upper = 5300108)
assert_integerish(params$state_codes, lower = 11, upper = 53)
assert_list(plot_style)
assert_number(plot_expansion, lower = 0)

Set Data Directories

accessibility_data_dir <- here("data-input/accessibility")
data_dir <- here("data")
input_data_dir <- here("data-input")
plot_dir <- here("plots")
rais_data_dir <- here("data-input/rais")
rfb_data_dir <- here("data-input/rfb")
raw_data_dir <- here("data-raw")

for (i in c(raw_data_dir, input_data_dir, data_dir, rfb_data_dir, rais_data_dir, accessibility_data_dir, plot_dir)) {
  if (!dir_exists(i)) {
    dir_create(i, recurse = TRUE)
  }
}

Download, Import, Process, and Plot Municipalities Data

See the Source of Data section for more information.

Download and Save Data

read_municipalities_file <-
  input_data_dir |>
  file.path(
    paste0(
      paste(
        "read_municipality",
        "code",
        "all",
        "year",
        params$year |> closest_geobr_year(type = "municipality"),
        "simplified",
        FALSE,
        sep = "-"
      ),
      ".rds"
    )
  )

if (!file_exists(read_municipalities_file)) {
  municipalities_data <-
    read_municipality(
      year = params$year |> closest_geobr_year(type = "municipality"),
      simplified = FALSE,
      showProgress = FALSE
    )

  municipalities_data |> write_rds(read_municipalities_file)
}

Visualize Municipalities and Background Shapes

Municipalities Shapes

Import Data

municipalities_data <- read_municipalities_file |> read_rds()

municipalities_data |> glimpse()
#> Rows: 5,570
#> Columns: 9
#> $ code_muni    <dbl> 1100015, 1100023, 1100031, 1100049, 1100056, 1100064, …
#> $ name_muni    <chr> "Alta Floresta d'Oeste", "Ariquemes", "Cabixi", "Cacoa…
#> $ code_state   <dbl> 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11…
#> $ abbrev_state <chr> "RO", "RO", "RO", "RO", "RO", "RO", "RO", "RO", "RO", …
#> $ name_state   <chr> "Rondônia", "Rondônia", "Rondônia", "Rondônia", "Rondô…
#> $ code_region  <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
#> $ name_region  <chr> "Norte", "Norte", "Norte", "Norte", "Norte", "Norte", …
#> $ year         <dbl> 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, …
#> $ geometry     <MULTIPOLYGON [°]> MULTIPOLYGON (((-62.0080637..., MULTIPOLY…

Filter and Transform Data

municipalities_data <-
  municipalities_data |>
  filter(code_muni %in% params$municipality_codes) |> 
  st_transform(st_crs(4326)) |>
  select(
    municipality = name_muni,
    municipality_code = code_muni,
    state = name_state,
    state_code = code_state,
    federal_unit = abbrev_state,
    region = name_region,
    region_code = code_region,
    geometry
  ) |>
  mutate(
    municipality_code = as.integer(municipality_code),
    state_code = as.integer(state_code),
    region_code = as.integer(region_code)
  )

municipalities_data |> glimpse()
#> Rows: 5
#> Columns: 8
#> $ municipality      <chr> "Belém", "João Pessoa", "Belo Horizonte", "São Pa…
#> $ municipality_code <int> 1501402, 2507507, 3106200, 3550308, 4314902
#> $ state             <chr> "Pará", "Paraíba", "Minas Gerais", "São Paulo", "…
#> $ state_code        <int> 15, 25, 31, 35, 43
#> $ federal_unit      <chr> "PA", "PB", "MG", "SP", "RS"
#> $ region            <chr> "Norte", "Nordeste", "Sudeste", "Sudeste", "Sul"
#> $ region_code       <int> 1, 2, 3, 3, 4
#> $ geometry          <MULTIPOLYGON [°]> MULTIPOLYGON (((-48.5231116..., MULTIPOLYGON (((-…

municipalities_data_shapes <- list()

for (i in params$municipality_codes) {
  municipalities_data_shapes[[as.character(i)]] <- 
    municipalities_data |> 
    filter(municipality_code == i) |>
    st_transform(st_crs(4326))
}

for (i in params$municipality_codes){
  municipalities_data_shapes[[as.character(i)]] |> 
    glimpse()
}
#> Rows: 1
#> Columns: 8
#> $ municipality      <chr> "João Pessoa"
#> $ municipality_code <int> 2507507
#> $ state             <chr> "Paraíba"
#> $ state_code        <int> 25
#> $ federal_unit      <chr> "PB"
#> $ region            <chr> "Nordeste"
#> $ region_code       <int> 2
#> $ geometry          <MULTIPOLYGON [°]> MULTIPOLYGON (((-34.8558491...
#> Rows: 1
#> Columns: 8
#> $ municipality      <chr> "Belo Horizonte"
#> $ municipality_code <int> 3106200
#> $ state             <chr> "Minas Gerais"
#> $ state_code        <int> 31
#> $ federal_unit      <chr> "MG"
#> $ region            <chr> "Sudeste"
#> $ region_code       <int> 3
#> $ geometry          <MULTIPOLYGON [°]> MULTIPOLYGON (((-43.869603 ...
#> Rows: 1
#> Columns: 8
#> $ municipality      <chr> "São Paulo"
#> $ municipality_code <int> 3550308
#> $ state             <chr> "São Paulo"
#> $ state_code        <int> 35
#> $ federal_unit      <chr> "SP"
#> $ region            <chr> "Sudeste"
#> $ region_code       <int> 3
#> $ geometry          <MULTIPOLYGON [°]> MULTIPOLYGON (((-46.6351157...
#> Rows: 1
#> Columns: 8
#> $ municipality      <chr> "Porto Alegre"
#> $ municipality_code <int> 4314902
#> $ state             <chr> "Rio Grande do Sul"
#> $ state_code        <int> 43
#> $ federal_unit      <chr> "RS"
#> $ region            <chr> "Sul"
#> $ region_code       <int> 4
#> $ geometry          <MULTIPOLYGON [°]> MULTIPOLYGON (((-51.1629204...
#> Rows: 1
#> Columns: 8
#> $ municipality      <chr> "Belém"
#> $ municipality_code <int> 1501402
#> $ state             <chr> "Pará"
#> $ state_code        <int> 15
#> $ federal_unit      <chr> "PA"
#> $ region            <chr> "Norte"
#> $ region_code       <int> 1
#> $ geometry          <MULTIPOLYGON [°]> MULTIPOLYGON (((-48.5231116...

Plot Data

for (i in params$municipality_codes) {
  plot <-  
    ggplot() +
      geom_sf(data = municipalities_data_shapes[[as.character(i)]]) +
      scale_x_continuous(labels = NULL, breaks = NULL) +
      scale_y_continuous(labels = NULL, breaks = NULL) +
      annotation_scale(
        aes(),
        location = "br",
        style = "tick",
        height = unit(0.5, "lines")
      ) +
      annotation_north_arrow(
        location = "br",
        height = unit(2, "lines"),
        width = unit(2, "lines"),
        pad_x = unit(0.25, "lines"),
        pad_y = unit(1.25, "lines"),
        style = north_arrow_fancy_orienteering
      ) +
      labs(
        title = paste0(
          municipalities_data |> 
          filter(municipality_code == i) |>
          pull(municipality),
          "-",
          municipalities_data |>
            filter(municipality_code == i) |>
            pull(federal_unit),
          " (",
          params$year |>
            closest_geobr_year(type = "municipality"),
          ")"
        ),
        caption = "Source: geobr R package."
      )
  print(plot)
}

Background Municipalities Shape

Import Data

background_municipalities_data <- read_municipalities_file |> read_rds()

background_municipalities_data |> glimpse()
#> Rows: 5,570
#> Columns: 9
#> $ code_muni    <dbl> 1100015, 1100023, 1100031, 1100049, 1100056, 1100064, …
#> $ name_muni    <chr> "Alta Floresta d'Oeste", "Ariquemes", "Cabixi", "Cacoa…
#> $ code_state   <dbl> 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11…
#> $ abbrev_state <chr> "RO", "RO", "RO", "RO", "RO", "RO", "RO", "RO", "RO", …
#> $ name_state   <chr> "Rondônia", "Rondônia", "Rondônia", "Rondônia", "Rondô…
#> $ code_region  <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
#> $ name_region  <chr> "Norte", "Norte", "Norte", "Norte", "Norte", "Norte", …
#> $ year         <dbl> 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, …
#> $ geometry     <MULTIPOLYGON [°]> MULTIPOLYGON (((-62.0080637..., MULTIPOLY…

Filter and Transform Data

background_municipalities_data <-
  background_municipalities_data |>
  filter(code_state %in% params$state_codes)|>
    select(
      municipality = name_muni,
      municipality_code = code_muni,
      state = name_state,
      state_code = code_state,
      federal_unit = abbrev_state,
      region = name_region,
      region_code = code_region,
      geometry
    ) |>
    mutate(
      municipality_code = as.integer(municipality_code),
      state_code = as.integer(state_code),
      region_code = as.integer(region_code)
    )

background_municipalities_shapes <- list()

for (i in params$municipality_codes) {
  background_municipalities_shapes[[as.character(i)]] <-
    background_municipalities_data |> 
    st_transform(st_crs(4326)) |>
    st_crop(
      y = municipalities_data_shapes[[as.character(i)]] |>
        st_bbox() |>
        expand_bbox(ratio = plot_expansion)
    ) |>
    suppressMessages()
}
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries

for (i in params$municipality_codes){
  background_municipalities_shapes[[as.character(i)]] |> 
    glimpse()
}
#> Rows: 6
#> Columns: 8
#> $ municipality      <chr> "Alhandra", "Bayeux", "Cabedelo", "Conde", "João …
#> $ municipality_code <int> 2500601, 2501807, 2503209, 2504603, 2507507, 2513…
#> $ state             <chr> "Paraíba", "Paraíba", "Paraíba", "Paraíba", "Para…
#> $ state_code        <int> 25, 25, 25, 25, 25, 25
#> $ federal_unit      <chr> "PB", "PB", "PB", "PB", "PB", "PB"
#> $ region            <chr> "Nordeste", "Nordeste", "Nordeste", "Nordeste", "…
#> $ region_code       <int> 2, 2, 2, 2, 2, 2
#> $ geometry          <POLYGON [°]> POLYGON ((-34.94059883 -7.2..., POLYGON ((-34.936…
#> Rows: 11
#> Columns: 8
#> $ municipality      <chr> "Belo Horizonte", "Betim", "Brumadinho", "Contage…
#> $ municipality_code <int> 3106200, 3106705, 3109006, 3118601, 3129806, 3144…
#> $ state             <chr> "Minas Gerais", "Minas Gerais", "Minas Gerais", "…
#> $ state_code        <int> 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31
#> $ federal_unit      <chr> "MG", "MG", "MG", "MG", "MG", "MG", "MG", "MG", "…
#> $ region            <chr> "Sudeste", "Sudeste", "Sudeste", "Sudeste", "Sude…
#> $ region_code       <int> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3
#> $ geometry          <GEOMETRY [°]> POLYGON ((-43.86970337 -19...., POLYGON ((-44.067…
#> Rows: 34
#> Columns: 8
#> $ municipality      <chr> "Arujá", "Barueri", "Caieiras", "Cajamar", "Carap…
#> $ municipality_code <int> 3503901, 3505708, 3509007, 3509205, 3510609, 3513…
#> $ state             <chr> "São Paulo", "São Paulo", "São Paulo", "São Paulo…
#> $ state_code        <int> 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 3…
#> $ federal_unit      <chr> "SP", "SP", "SP", "SP", "SP", "SP", "SP", "SP", "…
#> $ region            <chr> "Sudeste", "Sudeste", "Sudeste", "Sudeste", "Sude…
#> $ region_code       <int> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3…
#> $ geometry          <GEOMETRY [°]> POLYGON ((-46.35567725 -23...., POLYGON …
#> Rows: 11
#> Columns: 8
#> $ municipality      <chr> "Alvorada", "Barra do Ribeiro", "Cachoeirinha", "…
#> $ municipality_code <int> 4300604, 4301909, 4303103, 4304606, 4306767, 4309…
#> $ state             <chr> "Rio Grande do Sul", "Rio Grande do Sul", "Rio Gr…
#> $ state_code        <int> 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43
#> $ federal_unit      <chr> "RS", "RS", "RS", "RS", "RS", "RS", "RS", "RS", "…
#> $ region            <chr> "Sul", "Sul", "Sul", "Sul", "Sul", "Sul", "Sul", …
#> $ region_code       <int> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4
#> $ geometry          <GEOMETRY [°]> POLYGON ((-51.00599962 -29...., POLYGON ((-51.309…
#> Rows: 12
#> Columns: 8
#> $ municipality      <chr> "Acará", "Ananindeua", "Barcarena", "Belém", "Ben…
#> $ municipality_code <int> 1500206, 1500800, 1501303, 1501402, 1501501, 1501…
#> $ state             <chr> "Pará", "Pará", "Pará", "Pará", "Pará", "Pará", "…
#> $ state_code        <int> 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15
#> $ federal_unit      <chr> "PA", "PA", "PA", "PA", "PA", "PA", "PA", "PA", "…
#> $ region            <chr> "Norte", "Norte", "Norte", "Norte", "Norte", "Nor…
#> $ region_code       <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
#> $ geometry          <GEOMETRY [°]> POLYGON ((-48.28942722 -1.5..., POLYGON ((-48.395…

Plot Background Shapes

for (i in params$municipality_codes) {
  plot <-
    ggplot() +
      geom_sf(data = background_municipalities_shapes[[as.character(i)]]) +
      coord_sf(expand = FALSE) +
      scale_x_continuous(labels = NULL, breaks = NULL) +
      scale_y_continuous(labels = NULL, breaks = NULL) +
      annotation_scale(
        location = "br",
        style = "tick",
        height = unit(0.5, "lines")
      ) +
      annotation_north_arrow(
        location = "br",
        height = unit(2, "lines"),
        width = unit(2, "lines"),
        pad_x = unit(0.25, "lines"),
        pad_y = unit(1.25, "lines"),
        style = north_arrow_fancy_orienteering
      ) +
      labs(
        title = paste0(
          municipalities_data |> 
            filter(municipality_code == i) |>
            pull(municipality),
          "-",
          municipalities_data |>
            filter(municipality_code == i) |>
            pull(federal_unit),
          " (",
          params$year |>
            closest_geobr_year(type = "municipality"),
          ")"
        ),
        caption = "Source: geobr R package."
      ) +
      theme(
        panel.background = element_rect(fill = "#e5e5e5")
      )
  print(plot)
}

Download, Import, Process, and Visualize Land Cover Shapes

See the Source of Data section for more information.

Download Data

mapbiomas_file <- here(
  input_data_dir,
  paste0(
    "mapbiomas_10m_collection2_integration_v1-classification_",
    params$year,
    ".tif"
  )
)

if (!file_exists(mapbiomas_file)) {
  request("https://storage.googleapis.com") |>
    req_url_path_append("mapbiomas-public") |>
    req_url_path_append("initiatives") |>
    req_url_path_append("brasil") |>
    req_url_path_append("lulc_10m") |>
    req_url_path_append("collection_2") |>
    req_url_path_append("integration") |>
    req_url_path_append(basename(mapbiomas_file)) |>
    req_user_agent("https://cem-usp.github.io/kernel-accessibility") |>
    req_options(connecttimeout = 10) |>
    req_retry(
      max_tries = 3,
      retry_on_failure = TRUE,
      backoff = \(attempt) 5^attempt
    ) |>
    req_progress() |>
    req_perform(path = mapbiomas_file)
}

Import Data

mapbiomas_data <- mapbiomas_file |> read_stars(proxy = TRUE)

Crop Data

mapbiomas_background_municipalities_shapes <- list()

crop_bboxes <- list()

for (i in params$municipality_codes) {
  crop_bboxes[[as.character(i)]] <-
    background_municipalities_data |>
    st_transform(st_crs(4326)) |>
    filter(municipality_code == i) |>
    st_bbox() |>
    expand_bbox(ratio = plot_expansion)

  mapbiomas_background_municipalities_shapes[[as.character(i)]] <-
    mapbiomas_data |>
    st_crop(crop_bboxes[[as.character(i)]])
}

Transform Data

Click here to see the legend codes for MapBiomas Collection 2.0.

green_areas_shapes <- list()

for (i in params$municipality_codes) {
  green_areas_shapes[[as.character(i)]] <-
    mapbiomas_background_municipalities_shapes[[as.character(i)]] |>
    set_names("value") |>
    mutate(
      value = if_else(
        value %in% c(1:8, 49),
        value,
        NA_integer_
      )
    ) |>
    st_as_sf(merge = TRUE, na.rm = TRUE) |>
    mutate(category = "Green Areas") |>
    st_transform(st_crs(4326))
}

water_bodies_shapes <- list()

for (i in params$municipality_codes) {
  water_bodies_shapes[[as.character(i)]] <-
    mapbiomas_background_municipalities_shapes[[as.character(i)]] |>
    set_names("value") |>
    mutate(
      value = if_else(
        value %in% c(26, 31, 27, 33),
        value,
        NA_integer_
      )
    ) |>
    st_as_sf(merge = TRUE, na.rm = TRUE) |>
    mutate(category = "Water Bodies") |>
    st_transform(st_crs(4326))
}

land_cover_shapes <- list()

for (i in params$municipality_codes){
  land_cover_shapes[[as.character(i)]] <-
  green_areas_shapes[[as.character(i)]] |>
  bind_rows(water_bodies_shapes[[as.character(i)]]) |>
  mutate(category = category |> factor())
}

Plot Data

for (i in params$municipality_codes){
  plot <-
    ggplot() +
      geom_sf(
        data = background_municipalities_shapes[[as.character(i)]],
        color = NA,
        fill = plot_style$land
      ) +
      geom_sf(
        data = land_cover_shapes[[as.character(i)]],
        aes(fill = category),
        color = NA
      ) +
      geom_sf(
        data = municipalities_data_shapes[[as.character(i)]],
        color = NA,
        fill = "black",
        alpha = 0.1
      ) +
      scale_fill_manual(
        name = NULL,
        values = c(
          "Green Areas" = plot_style$green_areas,
          "Water Bodies" = plot_style$water_bodies
        )
      ) +
      coord_sf(expand = FALSE) +
      scale_x_continuous(labels = NULL, breaks = NULL) +
      scale_y_continuous(labels = NULL, breaks = NULL) +
      annotation_scale(
        location = "br",
        style = "tick",
        height = unit(0.5, "lines")
      ) +
      annotation_north_arrow(
        location = "br",
        height = unit(2, "lines"),
        width = unit(2, "lines"),
        pad_x = unit(0.25, "lines"),
        pad_y = unit(1.25, "lines"),
        style = north_arrow_fancy_orienteering
      ) +
      labs(
        title = paste0(
          municipalities_data_shapes[[as.character(i)]] |> pull(municipality),
          "-",
          municipalities_data_shapes[[as.character(i)]] |> pull(federal_unit),
          " (",
          params$year |>
            closest_geobr_year(type = "municipality"),
          ")"
        ),
        caption = "Source: Mapbiomas Collection 2.0."
      )
    print(plot)
}

Download, Import, Process, and Plot Favela’s and Urban Communitie’s Data

See the Source of Data section for more information.

Download and Save Data

read_fcu_file <-
  input_data_dir |>
  file.path(
    paste0(
      paste(
        "read_fcu",
        "code",
        "all",
        "year",
        params$year |> closest_geobr_year(type = "favela"),
        "simplified",
        FALSE,
        sep = "-"
      ),
      ".rds"
    )
  )

if (!file_exists(read_fcu_file)) {
  fcu_data <-
    read_favela(
      year = params$year |> closest_geobr_year(type = "favela"),
      simplified = FALSE,
      showProgress = FALSE
    )

  fcu_data |> write_rds(read_fcu_file)
}

Import Data

fcu_data <- read_fcu_file |> read_rds()

fcu_data |> glimpse()
#> Rows: 12,348
#> Columns: 11
#> $ code_muni    <dbl> 1100049, 1100049, 1100106, 1100106, 1100106, 1100106, …
#> $ name_muni    <chr> "Cacoal", "Cacoal", "Guajará-Mirim", "Guajará-Mirim", …
#> $ code_favela  <dbl> 11000490002, 11000490006, 11001060003, 11001060004, 11…
#> $ name_favela  <chr> "Avenida 07 de Setembro", "Vitória", "Santa Luzia", "N…
#> $ code_state   <dbl> 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11…
#> $ abbrev_state <chr> "RO", "RO", "RO", "RO", "RO", "RO", "RO", "RO", "RO", …
#> $ name_state   <chr> "Rondônia", "Rondônia", "Rondônia", "Rondônia", "Rondô…
#> $ code_region  <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
#> $ name_region  <chr> "Norte", "Norte", "Norte", "Norte", "Norte", "Norte", …
#> $ year         <dbl> 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, …
#> $ geometry     <MULTIPOLYGON [°]> MULTIPOLYGON (((-61.4446486..., MULTIPOLY…

Filter and Transform Data

fcu_data <-
  fcu_data |>
  filter(code_muni %in% params$municipality_codes)

fcu_data <-
  fcu_data |>
  st_transform(st_crs(4326)) |>
  select(
    municipality = name_muni,
    municipality_code = code_muni,
    favela_code = code_favela,
    favela = name_favela,
    state = name_state,
    state_code = code_state,
    federal_unit = abbrev_state,
    region = name_region,
    region_code = code_region
  )

fcu_data |> glimpse()
#> Rows: 2,073
#> Columns: 10
#> $ municipality      <chr> "Belém", "Belém", "Belém", "Belém", "Belém", "Bel…
#> $ municipality_code <dbl> 1501402, 1501402, 1501402, 1501402, 1501402, 1501…
#> $ favela_code       <dbl> 15014020001, 15014020002, 15014020003, 1501402000…
#> $ favela            <chr> "Bacia do Una-Telégrafo", "Bacia do Una-Sacrament…
#> $ state             <chr> "Pará", "Pará", "Pará", "Pará", "Pará", "Pará", "…
#> $ state_code        <dbl> 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 1…
#> $ federal_unit      <chr> "PA", "PA", "PA", "PA", "PA", "PA", "PA", "PA", "…
#> $ region            <chr> "Norte", "Norte", "Norte", "Norte", "Norte", "Nor…
#> $ region_code       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ geometry          <MULTIPOLYGON [°]> MULTIPOLYGON (((-48.481316 ..., MULT…

Plot Data

Prepare Data

fcu_data_shapes <- list()

for (i in params$municipality_codes) {
  fcu_data_shapes[[as.character(i)]] <- 
    fcu_data |> 
    filter(municipality_code == i) |>
    st_transform(st_crs(4326))
}

for (i in params$municipality_codes){
  fcu_data_shapes[[as.character(i)]] |> 
    glimpse()
}
#> Rows: 157
#> Columns: 10
#> $ municipality      <chr> "João Pessoa", "João Pessoa", "João Pessoa", "Joã…
#> $ municipality_code <dbl> 2507507, 2507507, 2507507, 2507507, 2507507, 2507…
#> $ favela_code       <dbl> 25075070001, 25075070002, 25075070003, 2507507000…
#> $ favela            <chr> "Maria de Nazaré", "Boa Esperança - Pedra Branca"…
#> $ state             <chr> "Paraíba", "Paraíba", "Paraíba", "Paraíba", "Para…
#> $ state_code        <dbl> 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 2…
#> $ federal_unit      <chr> "PB", "PB", "PB", "PB", "PB", "PB", "PB", "PB", "…
#> $ region            <chr> "Nordeste", "Nordeste", "Nordeste", "Nordeste", "…
#> $ region_code       <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
#> $ geometry          <MULTIPOLYGON [°]> MULTIPOLYGON (((-34.8837568..., MULT…
#> Rows: 218
#> Columns: 10
#> $ municipality      <chr> "Belo Horizonte", "Belo Horizonte", "Belo Horizon…
#> $ municipality_code <dbl> 3106200, 3106200, 3106200, 3106200, 3106200, 3106…
#> $ favela_code       <dbl> 31062000001, 31062000002, 31062000004, 3106200000…
#> $ favela            <chr> "Vila Paquetá", "Vila Nova dos Milionários", "Vil…
#> $ state             <chr> "Minas Gerais", "Minas Gerais", "Minas Gerais", "…
#> $ state_code        <dbl> 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 3…
#> $ federal_unit      <chr> "MG", "MG", "MG", "MG", "MG", "MG", "MG", "MG", "…
#> $ region            <chr> "Sudeste", "Sudeste", "Sudeste", "Sudeste", "Sude…
#> $ region_code       <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3…
#> $ geometry          <MULTIPOLYGON [°]> MULTIPOLYGON (((-43.9925751..., MULT…
#> Rows: 1,359
#> Columns: 10
#> $ municipality      <chr> "São Paulo", "São Paulo", "São Paulo", "São Paulo…
#> $ municipality_code <dbl> 3550308, 3550308, 3550308, 3550308, 3550308, 3550…
#> $ favela_code       <dbl> 35503080005, 35503080007, 35503080009, 3550308001…
#> $ favela            <chr> "Mirangoaba Ii", "Dezenove", "Jardim Cotinha", "S…
#> $ state             <chr> "São Paulo", "São Paulo", "São Paulo", "São Paulo…
#> $ state_code        <dbl> 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 3…
#> $ federal_unit      <chr> "SP", "SP", "SP", "SP", "SP", "SP", "SP", "SP", "…
#> $ region            <chr> "Sudeste", "Sudeste", "Sudeste", "Sudeste", "Sude…
#> $ region_code       <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3…
#> $ geometry          <MULTIPOLYGON [°]> MULTIPOLYGON (((-46.680735 ..., MULT…
#> Rows: 125
#> Columns: 10
#> $ municipality      <chr> "Porto Alegre", "Porto Alegre", "Porto Alegre", "…
#> $ municipality_code <dbl> 4314902, 4314902, 4314902, 4314902, 4314902, 4314…
#> $ favela_code       <dbl> 43149020004, 43149020009, 43149020015, 4314902001…
#> $ favela            <chr> "Vila Dona Teodora", "Vila Mariano de Matos Ii", …
#> $ state             <chr> "Rio Grande do Sul", "Rio Grande do Sul", "Rio Gr…
#> $ state_code        <dbl> 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 4…
#> $ federal_unit      <chr> "RS", "RS", "RS", "RS", "RS", "RS", "RS", "RS", "…
#> $ region            <chr> "Sul", "Sul", "Sul", "Sul", "Sul", "Sul", "Sul", …
#> $ region_code       <dbl> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4…
#> $ geometry          <MULTIPOLYGON [°]> MULTIPOLYGON (((-51.1935989..., MULT…
#> Rows: 214
#> Columns: 10
#> $ municipality      <chr> "Belém", "Belém", "Belém", "Belém", "Belém", "Bel…
#> $ municipality_code <dbl> 1501402, 1501402, 1501402, 1501402, 1501402, 1501…
#> $ favela_code       <dbl> 15014020001, 15014020002, 15014020003, 1501402000…
#> $ favela            <chr> "Bacia do Una-Telégrafo", "Bacia do Una-Sacrament…
#> $ state             <chr> "Pará", "Pará", "Pará", "Pará", "Pará", "Pará", "…
#> $ state_code        <dbl> 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 1…
#> $ federal_unit      <chr> "PA", "PA", "PA", "PA", "PA", "PA", "PA", "PA", "…
#> $ region            <chr> "Norte", "Norte", "Norte", "Norte", "Norte", "Nor…
#> $ region_code       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ geometry          <MULTIPOLYGON [°]> MULTIPOLYGON (((-48.481316 ..., MULT…

Layer on `municipailities_data_shapes`

for (i in params$municipality_codes) {
  plot <-  
    ggplot() +
      geom_sf(data = municipalities_data_shapes[[as.character(i)]]) +
      geom_sf(
        data = fcu_data_shapes[[as.character(i)]],
        color = "black",
        linewidth = 0.1
        ) + 
      scale_fill_manual(
        name = NULL,
        values = c(
          "Green Areas" = plot_style$green_areas,
          "Water Bodies" = plot_style$water_bodies
        )
      ) +
      scale_x_continuous(labels = NULL, breaks = NULL) +
      scale_y_continuous(labels = NULL, breaks = NULL) +
      annotation_scale(
        aes(),
        location = "br",
        style = "tick",
        height = unit(0.5, "lines")
      ) +
      annotation_north_arrow(
        location = "br",
        height = unit(2, "lines"),
        width = unit(2, "lines"),
        pad_x = unit(0.25, "lines"),
        pad_y = unit(1.25, "lines"),
        style = north_arrow_fancy_orienteering
      ) +
      labs(
        title = paste0(
          "Favelas and Urban Comunities in ",
          municipalities_data |> 
          filter(municipality_code == i) |>
          pull(municipality),
          "-",
          municipalities_data |>
            filter(municipality_code == i) |>
            pull(federal_unit),
          " (",
          params$year |>
            closest_geobr_year(type = "favela"),
          ")"
        ),
        caption = "Source: geobr R package."
      ) + 
      theme(plot.title = element_text(hjust = 0.5))
  print(plot)
}
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.

#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.

#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.

#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.

#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.

Layer on `background_municipalities_shapes`

for (i in params$municipality_codes) {
  plot <-
    ggplot() +
      geom_sf(data = background_municipalities_shapes[[as.character(i)]]) +
      geom_sf(
        data = fcu_data_shapes[[as.character(i)]],
        color = "black",
        linewidth = 0.1
        ) +
      coord_sf(expand = FALSE) +
      scale_x_continuous(labels = NULL, breaks = NULL) +
      scale_y_continuous(labels = NULL, breaks = NULL) +
      annotation_scale(
        location = "br",
        style = "tick",
        height = unit(0.5, "lines")
      ) +
      annotation_north_arrow(
        location = "br",
        height = unit(2, "lines"),
        width = unit(2, "lines"),
        pad_x = unit(0.25, "lines"),
        pad_y = unit(1.25, "lines"),
        style = north_arrow_fancy_orienteering
      ) +
      labs(
        title = paste0(
          "Favelas and Urban Comuunities in ",
          municipalities_data |> 
            filter(municipality_code == i) |>
            pull(municipality),
          "-",
          municipalities_data |>
            filter(municipality_code == i) |>
            pull(federal_unit),
          " (",
          params$year |>
            closest_geobr_year(type = "favela"),
          ")"
        ),
        caption = "Source: geobr R package."
      ) +
      theme(
        plot.title = element_text(hjust = 0.5),
        panel.background = element_rect(fill = "#e5e5e5")
      )
  print(plot)
}

Layer on `landcover_shapes`

for (i in params$municipality_codes){
  plot <-
    ggplot() +
      geom_sf(
        data = background_municipalities_shapes[[as.character(i)]],
        color = NA,
        fill = plot_style$land
      ) +
      geom_sf(
        data = land_cover_shapes[[as.character(i)]],
        aes(fill = category),
        color = NA
      ) +
      geom_sf(
        data = municipalities_data_shapes[[as.character(i)]],
        color = NA,
        fill = "black",
        alpha = 0.1
      ) +
      geom_sf(
        data = fcu_data_shapes[[as.character(i)]],
        color = NA,
        fill = "black",
        alpha = 0.75
      ) +
      scale_fill_manual(
        name = NULL,
        values = c(
          "Green Areas" = plot_style$green_areas,
          "Water Bodies" = plot_style$water_bodies
        )
      ) +
      coord_sf(expand = FALSE) +
      scale_x_continuous(labels = NULL, breaks = NULL) +
      scale_y_continuous(labels = NULL, breaks = NULL) +
      annotation_scale(
        location = "br",
        style = "tick",
        height = unit(0.5, "lines")
      ) +
      annotation_north_arrow(
        location = "br",
        height = unit(2, "lines"),
        width = unit(2, "lines"),
        pad_x = unit(0.25, "lines"),
        pad_y = unit(1.25, "lines"),
        style = north_arrow_fancy_orienteering
      ) +
      labs(
        title = paste0(
          "Favelas and Urban Comunities in ",
          municipalities_data_shapes[[as.character(i)]] |> pull(municipality),
          "-",
          municipalities_data_shapes[[as.character(i)]] |> pull(federal_unit),
          " (",
          params$year |>
            closest_geobr_year(type = "favela"),
          ")"
        ),
        caption = "Source: Mapbiomas Collection 2.0."
      ) + 
      theme(plot.title = element_text(hjust = 0.5))
    print(plot)
}

Set Plot Template

Define Template

plot_template <- function(coord_sf = TRUE) {
  assert_flag(coord_sf)

  out <- list()

  for (i in params$municipality_codes) {
    out[[as.character(i)]] <-
      ggplot() +
      geom_sf(
        data = background_municipalities_shapes[[as.character(i)]],
        color = NA,
        fill = plot_style$land
      ) +
      geom_sf(
        data = land_cover_shapes[[as.character(i)]],
        aes(fill = category),
        color = NA
      ) +
      scale_fill_manual(
        name = NULL,
        values = c(
          "Green Areas" = plot_style$green_areas,
          "Water Bodies" = plot_style$water_bodies
        ),
        guide = "none"
      ) +
      new_scale_fill() +
      geom_sf(
        data = municipalities_data_shapes[[as.character(i)]],
        color = NA,
        fill = "black",
        alpha = 0.1
      ) + 
      geom_sf(
          data = fcu_data_shapes[[as.character(i)]],
          color = NA,
          fill = "black",
          alpha = 0.75
      ) +
      scale_x_continuous(labels = NULL, breaks = NULL) +
      scale_y_continuous(labels = NULL, breaks = NULL) +
      annotation_scale(
        location = "br",
        style = "tick",
        height = unit(0.5, "lines")
      ) +
      annotation_north_arrow(
        location = "br",
        height = unit(2, "lines"),
        width = unit(2, "lines"),
        pad_x = unit(0.25, "lines"),
        pad_y = unit(1.25, "lines"),
        style = north_arrow_fancy_orienteering
      ) +
      theme(
        legend.box = "horizontal",
        legend.box.just = "left",
        legend.direction = "horizontal",
        legend.text = element_text(size = 8),
        plot.title = element_text(hjust = 0.5)
      )

    if (isTRUE(coord_sf)) {
      out[[as.character(i)]] <- out[[as.character(i)]] + coord_sf(expand = FALSE)
    }
  }

  out
}

Plot Template

plot_template()
#> $`2507507`

#> 
#> $`3106200`

#> 
#> $`3550308`

#> 
#> $`4314902`

#> 
#> $`1501402`

Download, Import, and Visualize Establishments Data from Bazilian Federal Revenue (RFB)

Download Data

osf_rfb_data_id <- "2x6jb"

osf_rfb_data_file <-
  osf_rfb_data_id |>
  osf_retrieve_node() |>
  osf_ls_files(
    type = "file",
    pattern = rfb_file_pattern,
    n_max = Inf
  )

osf_rfb_data_file

rfb_file <- dir_ls(
  rfb_data_dir, 
  regexp = "\\.(csv|lockr)$"
  ) |>
  pluck(
    1, 
    .default = file.path(
      rfb_data_dir, 
      "placeholder")
      )

rfb_file <-
  osf_rfb_data_file |>
  osf_download(
    path = rfb_data_dir,
    conflicts = "overwrite"
  ) |>
  pull("local_path")

Unlock Files

rfb_file <- dir_ls(here(rfb_data_dir), regexp = "\\.lockr$") |>
  map_chr(
    \(x) unlock_file(
      file = x,
      private_key = here(".ssh", "id_rsa"),
      remove_file = TRUE,
      password = password
    )
  )

Import Data

rfb_data <- 
  dir_ls(
    here(rfb_data_dir), 
    regexp = "2026-01\\.rds$"
  ) |>
  read_rds()

rfb_data |> glimpse()
#> Rows: 286,759
#> Columns: 24
#> $ cnpj                       <chr> "07.755.352/0002-87", "13.121.990/0001-6…
#> $ legal_name                 <chr> "S M COMERCIO DE CARNES LTDA", "ACOUGUE …
#> $ trade_name                 <chr> "BIG BOI", "ACOUGUE CURTICAO DAS CARNES"…
#> $ legal_nature_code          <int> 2062, 2062, 2062, 2062, 2062, 2135, 2062…
#> $ share_capital              <dbl> 20000, 110000, 100000, 20000, 300000, 10…
#> $ opening_date               <date> 2007-06-19, 2011-01-18, 2011-05-03, 201…
#> $ simples_nacional           <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE…
#> $ primary_cnae               <chr> "47.22-9/01", "47.22-9/01", "47.22-9/01"…
#> $ secondary_cnae             <chr> "53.20-2/02,47.23-7/00", "47.21-1/03,47.…
#> $ street_type                <chr> "TRAVESSA", "TRAVESSA", "RUA", "TRAVESSA…
#> $ street_name                <chr> "TEOFILO CONDURU", "SIQUEIRA MENDES", "B…
#> $ number                     <chr> "940", "32", "140", "91", "582", "20 A",…
#> $ complement                 <chr> NA, NA, "LOJA", "BOX 57", NA, NA, NA, NA…
#> $ neighborhood               <chr> "CANUDOS", "VILA (MOSQUEIRO)", "BENGUI",…
#> $ postal_code                <chr> "66073-570", "66910-110", "66630-140", "…
#> $ municipality               <chr> "BELÉM", "BELÉM", "BELÉM", "BELÉM", "BEL…
#> $ municipality_code          <int> 1501402, 1501402, 1501402, 1501402, 1501…
#> $ latitude                   <dbl> -1.455280487, -1.165107000, -1.376883091…
#> $ longitude                  <dbl> -48.46080877, -48.47249900, -48.45089664…
#> $ geocodebr_precision        <fct> numero_aproximado, numero, numero, numer…
#> $ geocodebr_type_of_result   <chr> "da01", "dn02", "dn01", "da02", "da01", …
#> $ geocodebr_deviation_meters <int> 7, 6, 88, 6, 10, 63, 6, 6, 276, 15, 9, 6…
#> $ geocodebr_address_found    <chr> "TRAVESSA TEOFILO CONDURU, 940 (aprox) -…
#> $ locais_nova_group          <ord> G0, G0, G0, G0, G0, G0, G0, G0, G0, G0, …

Transform and Filter Data

rfb_data_fcu <- rfb_data |> 
  st_as_sf(
    crs = 4326,
    coords = c("longitude", "latitude")
    ) |> 
  st_transform(st_crs(fcu_data)) |>
  st_join(
    fcu_data,
    join = st_intersects,
    left = TRUE
    )

rfb_data_fcu <- rfb_data_fcu |>
  st_transform(4326) |>
  mutate(
    longitude = st_coordinates(geometry)[, 1], 
    latitude  = st_coordinates(geometry)[, 2]  
  ) |>
  st_drop_geometry() |>
  select(
    cnpj,
    legal_name,
    trade_name,
    legal_nature_code,
    share_capital,
    opening_date,
    simples_nacional,
    primary_cnae,
    secondary_cnae,
    street_type,
    street_name,
    number,
    complement,
    neighborhood,
    postal_code,
    municipality.x,
    municipality_code.x,
    favela_code,
    favela,
    latitude,
    longitude,
    geocodebr_precision,
    geocodebr_type_of_result,
    geocodebr_deviation_meters,
    geocodebr_address_found,
    locais_nova_group
  ) |>
  rename(
    municipality = municipality.x,
    municipality_code = municipality_code.x
  ) |> 
  st_as_sf(
    coords = c("longitude", "latitude"),
    crs = 4326,
    remove = FALSE  
  )

rfb_data_fcu <- rfb_data_fcu |>
  filter(!is.na(favela))

rfb_data_fcu |> glimpse()
#> Rows: 30,001
#> Columns: 27
#> $ cnpj                       <chr> "07.755.352/0002-87", "20.430.696/0001-0…
#> $ legal_name                 <chr> "S M COMERCIO DE CARNES LTDA", "DISTRIBU…
#> $ trade_name                 <chr> "BIG BOI", "DISTRIBUIDORA DE CARNES GEMA…
#> $ legal_nature_code          <int> 2062, 2062, 2062, 2135, 2062, 2135, 2062…
#> $ share_capital              <dbl> 20000, 20000, 300000, 10000, 2023000, 10…
#> $ opening_date               <date> 2007-06-19, 2014-06-11, 2012-04-24, 201…
#> $ simples_nacional           <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE…
#> $ primary_cnae               <chr> "47.22-9/01", "47.22-9/01", "47.22-9/01"…
#> $ secondary_cnae             <chr> "53.20-2/02,47.23-7/00", "47.21-1/03", "…
#> $ street_type                <chr> "TRAVESSA", "TRAVESSA", "TRAVESSA", "PAS…
#> $ street_name                <chr> "TEOFILO CONDURU", "CRISTOVAO COLOMBO", …
#> $ number                     <chr> "940", "91", "582", "20 A", "159", "S/N"…
#> $ complement                 <chr> NA, "BOX 57", NA, NA, NA, "BOX 02", NA, …
#> $ neighborhood               <chr> "CANUDOS", "CRUZEIRO (ICOARACI)", "TELEG…
#> $ postal_code                <chr> "66073-570", "66810-000", "66083-000", "…
#> $ municipality               <chr> "BELÉM", "BELÉM", "BELÉM", "BELÉM", "BEL…
#> $ municipality_code          <int> 1501402, 1501402, 1501402, 1501402, 1501…
#> $ favela_code                <dbl> 15014020050, 15014020035, 15014020001, 1…
#> $ favela                     <chr> "Bacia do Tucunduba Canudos", "Parque Vi…
#> $ latitude                   <dbl> -1.455280487, -1.363515971, -1.424346988…
#> $ longitude                  <dbl> -48.46080877, -48.46399571, -48.47769878…
#> $ geocodebr_precision        <fct> numero_aproximado, numero_aproximado, nu…
#> $ geocodebr_type_of_result   <chr> "da01", "da02", "da01", "dc01", "dn02", …
#> $ geocodebr_deviation_meters <int> 7, 6, 10, 63, 6, 276, 15, 6, 7, 6, 623, …
#> $ geocodebr_address_found    <chr> "TRAVESSA TEOFILO CONDURU, 940 (aprox) -…
#> $ locais_nova_group          <ord> G0, G0, G0, G0, G0, G0, G0, G0, G0, G0, …
#> $ geometry                   <POINT [°]> POINT (-48.46080877 -1.4552..., PO…

Visualize Data

Prepare Data

rfb_data_plot_fcu <- list()

for (m in params$municipality_codes) {
  rfb_data_plot_fcu[[as.character(m)]] <-
    rfb_data_fcu |>
    filter(municipality_code == m)
}

locais_nova_groups <- c("g0", "g1_g2", "g3", "g4")

Build Base Maps

i_maps_rfb <- plot_template(coord_sf = FALSE)

Combine Data and Maps

i_plots_rfb <- list()
i_data_rfb <- list()

for (m in params$municipality_codes) {
  for (i in locais_nova_groups) {
    i_data_rfb[[as.character(m)]][[i]] <-
      rfb_data_plot_fcu[[as.character(m)]] |>
      filter(
        locais_nova_group |>
          equals(i |> str_replace("_", "-") |> str_to_upper())
      )

    i_plots_rfb[[as.character(m)]][[i]] <-
      i_maps_rfb[[as.character(m)]] +
      geom_sf(
        data = i_data_rfb[[as.character(m)]][[i]],
        aes(color = locais_nova_group),
        size = 0.3
      ) +
      scale_color_manual(
        values = group_colors,
        name = NULL,
        guide = guide_legend(order = 3, override.aes = list(size = 2))
      ) +
      coord_sf(expand = FALSE) +
      map_labs(
        municipality_code = m,
        title = "Locais-Nova Establishments in Favelas and Urban Comunities in",
        source = "Brazilian Federal Revenue Service (RFB) CNPJ database",
        note = str_glue("as of {month.name[params$month]} 2026")
      )
  }
}

Compute Plot Ratio

plot_ratio_rfb <-
  i_plots_rfb[[as.character(params$municipality_codes[1])]] |>
  pluck("g0") |>
  get_plot_ratio()

Save Plots

for (m in params$municipality_codes) {
  for (i in locais_nova_groups) {
    plot_dir |>
      file.path(
        paste0(
          "rfb-", m, "-",
          i |> str_replace_all("_", "-"),
          ".png"
        )
      ) |>
      ggsave(
        plot = i_plots_rfb[[as.character(m)]][[i]],
        width = 21,
        height = (21 * plot_ratio_rfb) - 1.75,
        units = "cm",
        dpi = 300,
        device = agg_png
      )
  }
}

Download, Import, and Visualize Establishmnts Data from Annual Social Information Report (RAIS)

Download Data

osf_rais_data_id <- "p5ufj"

osf_rais_data_file <-
  osf_rais_data_id |>
  osf_retrieve_node() |>
  osf_ls_files(
    type = "file",
    pattern = rais_file_pattern,
    n_max = Inf
  )

osf_rais_data_file

rais_file <- dir_ls(
  rais_data_dir, 
  regexp = "\\.(csv|lockr)$"
  ) |>
  pluck(
    1, 
    .default = file.path(
      rais_data_dir, 
      "placeholder")
      )

rais_file <-
  osf_rais_data_file |>
  osf_download(
    path = rais_data_dir,
    conflicts = "overwrite"
  ) |>
  pull("local_path")

Unlock Files

rais_file <- dir_ls(here(rais_data_dir), regexp = "\\.lockr$") |>
  map_chr(
    \(x) unlock_file(
      file = x,
      private_key = here(".ssh", "id_rsa"),
      remove_file = TRUE,
      password = password
    )
  )

Import Data

rais_data <- 
  dir_ls(
    here(rais_data_dir), 
    regexp = "2022\\.rds$"
  ) |>
  read_rds()

rais_data |> glimpse()
#> Rows: 56,159
#> Columns: 25
#> $ cnpj                       <chr> "45.097.533/0001-71", "14.703.079/0001-2…
#> $ legal_name                 <chr> "NÃO INFORMADO", "NÃO INFORMADO", "NÃO I…
#> $ legal_nature_code          <int> 2062, 2135, 2135, 2062, 2135, 2062, 2062…
#> $ opening_date               <date> 2022-02-01, 2011-11-30, 2010-09-23, 200…
#> $ simples_nacional           <lgl> FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, TR…
#> $ primary_cnae               <chr> "47.22-9/01", "47.22-9/01", "47.22-9/01"…
#> $ active_employees           <int> 2, 4, 0, 9, 0, 2, 7, 0, 2, 23, 0, 1, 2, …
#> $ employee_size_class        <fct> 1-4 EMPLOYEES, 1-4 EMPLOYEES, 0 EMPLOYEE…
#> $ clt_employees              <int> 2, 4, 0, 9, 0, 2, 7, 0, 2, 23, 0, 1, 2, …
#> $ statutory_employees        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ street_type                <chr> "AVENIDA", "TRAVESSA", "AVENIDA", "AVENI…
#> $ street_name                <chr> "ROMULO MAIORANA", "SIQUEIRA MENDES", "E…
#> $ number                     <chr> "231", "602", "606", "1146", "2172", "18…
#> $ complement                 <chr> "LOJA 1", NA, NA, NA, NA, NA, "CASA B", …
#> $ neighborhood               <chr> "MARCO", "MOSQUEIRO", "JURUNAS", "MARAMB…
#> $ postal_code                <chr> "66093-005", "66910-110", "66030-250", "…
#> $ municipality               <chr> "BELÉM", "BELÉM", "BELÉM", "BELÉM", "BEL…
#> $ municipality_code          <int> 1501402, 1501402, 1501402, 1501402, 1501…
#> $ latitude                   <dbl> -1.442252123, -1.165174137, -1.471696890…
#> $ longitude                  <dbl> -48.46634257, -48.47179267, -48.49118097…
#> $ geocodebr_precision        <fct> numero_aproximado, numero_aproximado, nu…
#> $ geocodebr_type_of_result   <chr> "da01", "da01", "da01", "dn01", "da04", …
#> $ geocodebr_deviation_meters <int> 6, 6, 12, 6, 11, 13, 7, 6, 453, 12, 6, 1…
#> $ geocodebr_address_found    <chr> "AVENIDA ROMULO MAIORANA, 231 (aprox) - …
#> $ locais_nova_group          <ord> G0, G0, G0, G0, G0, G0, G0, G0, G0, G0, …

Transform and Filter Data

rais_data_fcu <- rais_data |> 
  st_as_sf(
    crs = 4326,
    coords = c("longitude", "latitude")
    ) |> 
  st_transform(st_crs(fcu_data)) |>
  st_join(
    fcu_data,
    join = st_intersects,
    left = TRUE
    )

rais_data_fcu <- rais_data_fcu |>
  st_transform(4326) |>
  mutate(
    longitude = st_coordinates(geometry)[, 1], 
    latitude  = st_coordinates(geometry)[, 2]  
  ) |>
  st_drop_geometry() |>
  select(
    cnpj,
    legal_name,
    legal_nature_code,
    opening_date,
    simples_nacional,
    primary_cnae,
    active_employees,
    employee_size_class,
    clt_employees,
    statutory_employees,
    street_type,
    street_name,
    number,
    complement,
    neighborhood,
    postal_code,
    municipality.x,
    municipality_code.x,
    favela_code,
    favela,
    latitude,
    longitude,
    geocodebr_precision,
    geocodebr_type_of_result,
    geocodebr_deviation_meters,
    geocodebr_address_found,
    locais_nova_group
  ) |>
  rename(
    municipality = municipality.x,
    municipality_code = municipality_code.x
  ) |> 
  st_as_sf(
    coords = c("longitude", "latitude"),
    crs = 4326,
    remove = FALSE  
  )

rais_data_fcu |> glimpse()
#> Rows: 56,159
#> Columns: 28
#> $ cnpj                       <chr> "45.097.533/0001-71", "14.703.079/0001-2…
#> $ legal_name                 <chr> "NÃO INFORMADO", "NÃO INFORMADO", "NÃO I…
#> $ legal_nature_code          <int> 2062, 2135, 2135, 2062, 2135, 2062, 2062…
#> $ opening_date               <date> 2022-02-01, 2011-11-30, 2010-09-23, 200…
#> $ simples_nacional           <lgl> FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, TR…
#> $ primary_cnae               <chr> "47.22-9/01", "47.22-9/01", "47.22-9/01"…
#> $ active_employees           <int> 2, 4, 0, 9, 0, 2, 7, 0, 2, 23, 0, 1, 2, …
#> $ employee_size_class        <fct> 1-4 EMPLOYEES, 1-4 EMPLOYEES, 0 EMPLOYEE…
#> $ clt_employees              <int> 2, 4, 0, 9, 0, 2, 7, 0, 2, 23, 0, 1, 2, …
#> $ statutory_employees        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ street_type                <chr> "AVENIDA", "TRAVESSA", "AVENIDA", "AVENI…
#> $ street_name                <chr> "ROMULO MAIORANA", "SIQUEIRA MENDES", "E…
#> $ number                     <chr> "231", "602", "606", "1146", "2172", "18…
#> $ complement                 <chr> "LOJA 1", NA, NA, NA, NA, NA, "CASA B", …
#> $ neighborhood               <chr> "MARCO", "MOSQUEIRO", "JURUNAS", "MARAMB…
#> $ postal_code                <chr> "66093-005", "66910-110", "66030-250", "…
#> $ municipality               <chr> "BELÉM", "BELÉM", "BELÉM", "BELÉM", "BEL…
#> $ municipality_code          <int> 1501402, 1501402, 1501402, 1501402, 1501…
#> $ favela_code                <dbl> NA, NA, 15014020007, NA, NA, NA, NA, NA,…
#> $ favela                     <chr> NA, NA, "Baixadas da Estrada Nova Juruna…
#> $ latitude                   <dbl> -1.442252123, -1.165174137, -1.471696890…
#> $ longitude                  <dbl> -48.46634257, -48.47179267, -48.49118097…
#> $ geocodebr_precision        <fct> numero_aproximado, numero_aproximado, nu…
#> $ geocodebr_type_of_result   <chr> "da01", "da01", "da01", "dn01", "da04", …
#> $ geocodebr_deviation_meters <int> 6, 6, 12, 6, 11, 13, 7, 6, 453, 12, 6, 1…
#> $ geocodebr_address_found    <chr> "AVENIDA ROMULO MAIORANA, 231 (aprox) - …
#> $ locais_nova_group          <ord> G0, G0, G0, G0, G0, G0, G0, G0, G0, G0, …
#> $ geometry                   <POINT [°]> POINT (-48.46634257 -1.4422..., PO…

rais_data_fcu <- rais_data_fcu |>
  filter(!is.na(favela))

rais_data_fcu |> glimpse()
#> Rows: 2,055
#> Columns: 28
#> $ cnpj                       <chr> "12.573.555/0001-01", "04.422.772/0035-1…
#> $ legal_name                 <chr> "NÃO INFORMADO", "FRIGORIFICO SANTA CRUZ…
#> $ legal_nature_code          <int> 2135, 2062, 2135, 2135, 2135, 2135, 2062…
#> $ opening_date               <date> 2010-09-23, 2017-12-21, 2017-12-19, 201…
#> $ simples_nacional           <lgl> FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, FA…
#> $ primary_cnae               <chr> "47.22-9/01", "47.22-9/01", "47.22-9/01"…
#> $ active_employees           <int> 0, 2, 23, 0, 1, 8, 8, 1, 10, 12, 0, 5, 4…
#> $ employee_size_class        <fct> 0 EMPLOYEES, 1-4 EMPLOYEES, 20-49 EMPLOY…
#> $ clt_employees              <int> 0, 2, 23, 0, 1, 8, 8, 1, 10, 12, 0, 5, 4…
#> $ statutory_employees        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ street_type                <chr> "AVENIDA", "RUA", "AVENIDA", "AVENIDA", …
#> $ street_name                <chr> "ENGENHEIRO FERNANDO GUILHON", "RODOLFO …
#> $ number                     <chr> "606", "1081", "40", "2118", "527", "333…
#> $ complement                 <chr> NA, "CONJUNTO MEDICI I", "COM A", NA, "B…
#> $ neighborhood               <chr> "JURUNAS", "MARAMBAIA", "JURUNAS", "TELE…
#> $ postal_code                <chr> "66030-250", "66620-000", "66030-250", "…
#> $ municipality               <chr> "BELÉM", "BELÉM", "BELÉM", "BELÉM", "BEL…
#> $ municipality_code          <int> 1501402, 1501402, 1501402, 1501402, 1501…
#> $ favela_code                <dbl> 15014020007, 15014020037, 15014020007, 1…
#> $ favela                     <chr> "Baixadas da Estrada Nova Jurunas", "Ent…
#> $ latitude                   <dbl> -1.471696890, -1.398502701, -1.474153138…
#> $ longitude                  <dbl> -48.49118097, -48.45362036, -48.49310939…
#> $ geocodebr_precision        <fct> numero_aproximado, cep, numero_aproximad…
#> $ geocodebr_type_of_result   <chr> "da01", "dc01", "da01", "da01", "dn01", …
#> $ geocodebr_deviation_meters <int> 12, 453, 12, 6, 14, 6, 12, 12, 18, 6, 12…
#> $ geocodebr_address_found    <chr> "AVENIDA ENGENHEIRO FERNANDO GUILHON, 60…
#> $ locais_nova_group          <ord> G0, G0, G0, G0, G0, G0, G0, G0, G0, G0, …
#> $ geometry                   <POINT [°]> POINT (-48.49118097 -1.4716..., PO…

Visualize Data

Prepare Data

rais_data_plot_fcu <- list()

for (m in params$municipality_codes) {
  rais_data_plot_fcu[[as.character(m)]] <-
    rais_data_fcu |>
    filter(municipality_code == m)
}

Build Base Maps

i_maps_rais <- plot_template(coord_sf = FALSE)

Combine Data and Maps

i_plots_rais <- list()
i_data_rais <- list()

for (m in params$municipality_codes) {
  for (i in locais_nova_groups) {
    i_data_rais[[as.character(m)]][[i]] <-
      rais_data_plot_fcu[[as.character(m)]] |>
      filter(
        locais_nova_group |>
          equals(i |> str_replace("_", "-") |> str_to_upper())
      )

    i_plots_rais[[as.character(m)]][[i]] <-
      i_maps_rais[[as.character(m)]] +
      geom_sf(
        data = i_data_rais[[as.character(m)]][[i]],
        aes(color = locais_nova_group),
        size = 0.3
      ) +
      scale_color_manual(
        values = group_colors,
        name = NULL,
        guide = guide_legend(order = 3, override.aes = list(size = 2))
      ) +
      coord_sf(expand = FALSE) +
      map_labs(
        municipality_code = m,
        title = "Locais-Nova Establishments in Favelas and Urban Comunities in",
        source = "Annual Social Information Report (RAIS) CNPJ Database",
        note = str_glue("as of {params$year}")
      )
  }
}

Compute Plot Ratio

plot_ratio_rais <-
  i_plots_rais[[as.character(params$municipality_codes[1])]] |>
  pluck("g0") |>
  get_plot_ratio()

Save Plots

for (m in params$municipality_codes) {
  for (i in locais_nova_groups) {
    plot_dir |>
      file.path(
        paste0(
          "rais-", m, "-",
          i |> str_replace_all("_", "-"),
          ".png"
        )
      ) |>
      ggsave(
        plot = i_plots_rfb[[as.character(m)]][[i]],
        width = 21,
        height = (21 * plot_ratio_rais) - 1.75,
        units = "cm",
        dpi = 300,
        device = agg_png
      )
  }
}

Download, Import, and Visualize Accesibilitty Data from Bazilian Federal Revenue (RFB) and Annual Social Information Report (RAIS)

Annual Social Information Report (RAIS)

Download Data

osf_rais_accessibility_data_id <- "fb5dv"

osf_rais_accessibility_data_file <-
  osf_rais_accessibility_data_id |>
  osf_retrieve_node() |>
  osf_ls_files(
    type = "file",
    pattern = rais_accessibility_file_pattern,
    n_max = Inf
  )

osf_rais_accessibility_data_file

rais_accessibility_file <- dir_ls(
  accessibility_data_dir, 
  regexp = "\\.(csv|lockr)$"
  ) |>
  pluck(
    1, 
    .default = file.path(
      accessibility_data_dir, 
      "placeholder")
      )

rais_accessibility_file <-
  osf_rais_accessibility_data_file |>
  osf_download(
    path = accessibility_data_dir,
    conflicts = "overwrite"
  ) |>
  pull("local_path")

Unlock Files

rais_accessibility_file <- dir_ls(here(accessibility_data_dir), regexp = "\\.lockr$") |>
  map_chr(
    \(x) unlock_file(
      file = x,
      private_key = here(".ssh", "id_rsa"),
      remove_file = TRUE,
      password = password
    )
  )

Import Data

rais_accessibility_data <- 
  dir_ls(
    here(accessibility_data_dir), 
    regexp = "2022_3550308.rds"
  ) |>
  read_rds()

rais_accessibility_data |> glimpse()
#> Rows: 34,168
#> Columns: 7
#> $ h3_address             <chr> "89a81000003ffff", "89a81000003ffff", "89a81…
#> $ population             <dbl> 379.4540780, 379.4540780, 379.4540780, 379.4…
#> $ accessibility          <dbl> 1, 2, 0, 6, 2, 3, 1, 10, 1, 2, 3, 8, 1, 2, 0…
#> $ accessibility_per_1000 <dbl> 2.6353650095, 5.2707300191, 0.0000000000, 15…
#> $ locais_nova_group      <chr> "G0", "G1-G2", "G3", "G4", "G0", "G1-G2", "G…
#> $ percentile             <int> 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, …
#> $ cutoff                 <int> 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, …

Filter and Transform Data

rais_accessibility_data <-
  rais_accessibility_data |>
  cell_to_polygon(simple = FALSE) |>
  st_transform(st_crs(4326))

Visualize Data

Define Plot Template

rais_accessibility_plot_template <- function(coord_sf = TRUE) {
  assert_flag(coord_sf)

  out <-
    ggplot() +
    geom_sf(
      data = background_municipalities_shapes[["3550308"]],
      color = NA,
      fill = plot_style$land
    ) +
    geom_sf(
      data = land_cover_shapes[["3550308"]],
      aes(fill = category),
      color = NA
    ) +
    scale_fill_manual(
      name = NULL,
      values = c(
        "Green Areas" = plot_style$green_areas,
        "Water Bodies" = plot_style$water_bodies
      ),
      guide = "none"
    ) +
    new_scale_fill() +
    geom_sf(
      data = municipalities_data_shapes[["3550308"]],
      color = NA,
      fill = "black",
      alpha = 0.1
    ) +
    geom_sf(
        data = fcu_data_shapes[["3550308"]],
        color = NA,
        fill = "black",
        alpha = 0.75
    ) +
    scale_x_continuous(labels = NULL, breaks = NULL) +
    scale_y_continuous(labels = NULL, breaks = NULL) +
    annotation_scale(
      location = "br",
      style = "tick",
      height = unit(0.5, "lines"),
      text_cex = 2.4
    ) +
    annotation_north_arrow(
      location = "br",
      height = unit(2.5, "lines"),
      width = unit(2.5, "lines"),
      pad_x = unit(0.25, "lines"),
      pad_y = unit(1.25, "lines"),
      style = north_arrow_fancy_orienteering(text_size = 20)
    ) +
    theme(
      legend.box = "horizontal",
      legend.box.just = "left",
      legend.direction = "horizontal",
      legend.text = element_text(size = 25)
    )

  if (isTRUE(coord_sf)) {
    out <- out + coord_sf(expand = FALSE)
  }

  out
}

Prepare Data

breaks <- c(-Inf, 0, 3, 5, 10, 15, Inf)

# fmt: skip
category_labels <- c(
  "0", "1 - 3", "4 - 5", "6 - 10", "11 - 15", "Above 15"
)

fill_scale <- list(
  g0 = brewer.pal(6, "PuRd") |> rev(),
  g1_g2 = brewer.pal(6, "YlGn") |> rev(),
  g3 = brewer.pal(6, "YlGnBu"),
  g4 = brewer.pal(6, "YlOrRd")
)

category_cut <- function(x) {
  assert_numeric(x)

  x |>
    cut(
      breaks = breaks,
      labels = category_labels,
      include.lowest = TRUE,
      right = TRUE
    ) |>
    factor(
      ordered = TRUE,
      levels = category_labels
    )
}

rais_accessibility_plot_data <-
  rais_accessibility_data |>
  mutate(
    accessibility_category = accessibility_per_1000 |> category_cut()
  ) |>
  relocate(h3_address) |>
  relocate(
    accessibility_category,
    .after = accessibility_per_1000
  )

Plot Accessibility

rais_accessibility_plots <- list()

for (i in locais_nova_groups) {
  i_rais_accessibility_data <-
    rais_accessibility_plot_data |>
    filter(
      locais_nova_group |>
        equals(
          i |>
            str_replace("_", "-") |>
            str_to_upper()
        )
    )

  i_rais_accessibility_plot <-
    rais_accessibility_plot_template(coord_sf = FALSE) +
    geom_sf(
      data = i_rais_accessibility_data,
      aes(fill = accessibility_category),
      color = NA
    ) +
    scale_fill_manual(
      values = fill_scale[[i]],
      name = NULL,
      guide = guide_legend(
        order = 3,
        ncol = 1,
        override.aes = list(color = "black", linewidth = 0.25)
      )
    ) +
    coord_sf(expand = FALSE)

  rais_accessibility_plots <-
    rais_accessibility_plots |>
    inset2(i, i_rais_accessibility_plot)
}

Compute Plot Ratio

rais_accessibility_plot_ratio <-
  rais_accessibility_plots |>
  pluck("g0") |>
  get_plot_ratio()

Save Plots

for (i in locais_nova_groups) {
  plot_dir |>
    file.path(
      i |>
        str_replace_all("_", "-") %>%
        paste0("rais-accessibility-", ., ".png")
    ) |>
    ggsave(
      plot = rais_accessibility_plots[[i]],
      width = 21,
      height = (21 * rais_accessibility_plot_ratio) - 1.75,
      units = "cm",
      dpi = 300,
      device = agg_png
    )
}

Bazilian Federal Revenue (RFB)

Download Data

Unlock Files

Import Data

rfb_accessibility_data <- 
  dir_ls(
    here(accessibility_data_dir), 
    regexp = "rfb_2026_3550308.rds"
  ) |>
  read_rds()

rfb_accessibility_data |> glimpse()
#> Rows: 34,148
#> Columns: 7
#> $ h3_address             <chr> "89a81000003ffff", "89a81000003ffff", "89a81…
#> $ population             <dbl> 379.4540780, 379.4540780, 379.4540780, 379.4…
#> $ locais_nova_group      <chr> "G0", "G1-G2", "G3", "G4", "G0", "G1-G2", "G…
#> $ accessibility          <dbl> 8, 18, 3, 22, 11, 31, 13, 58, 4, 27, 15, 70,…
#> $ accessibility_per_1000 <dbl> 21.082920076, 47.436570172, 7.906095029, 57.…
#> $ percentile             <int> 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, …
#> $ cutoff                 <int> 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, …

Filter and Transform Data

rfb_accessibility_data <-
  rfb_accessibility_data |>
  cell_to_polygon(simple = FALSE) |>
  st_transform(st_crs(4326))

Visualize Data

Define Plot Template

rfb_accessibility_plot_template <- function(coord_sf = TRUE) {
  assert_flag(coord_sf)

  out <-
    ggplot() +
    geom_sf(
      data = background_municipalities_shapes[["3550308"]],
      color = NA,
      fill = plot_style$land
    ) +
    geom_sf(
      data = land_cover_shapes[["3550308"]],
      aes(fill = category),
      color = NA
    ) +
    scale_fill_manual(
      name = NULL,
      values = c(
        "Green Areas" = plot_style$green_areas,
        "Water Bodies" = plot_style$water_bodies
      ),
      guide = "none"
    ) +
    new_scale_fill() +
    geom_sf(
      data = municipalities_data_shapes[["3550308"]],
      color = NA,
      fill = "black",
      alpha = 0.1
    ) +
    geom_sf(
        data = fcu_data_shapes[["3550308"]],
        color = NA,
        fill = "black",
        alpha = 0.75
    ) +
    scale_x_continuous(labels = NULL, breaks = NULL) +
    scale_y_continuous(labels = NULL, breaks = NULL) +
    annotation_scale(
      location = "br",
      style = "tick",
      height = unit(0.5, "lines"),
      text_cex = 2.4
    ) +
    annotation_north_arrow(
      location = "br",
      height = unit(2.5, "lines"),
      width = unit(2.5, "lines"),
      pad_x = unit(0.25, "lines"),
      pad_y = unit(1.25, "lines"),
      style = north_arrow_fancy_orienteering(text_size = 20)
    ) +
    theme(
      legend.box = "horizontal",
      legend.box.just = "left",
      legend.direction = "horizontal",
      legend.text = element_text(size = 25)
    )

  if (isTRUE(coord_sf)) {
    out <- out + coord_sf(expand = FALSE)
  }

  out
}

Prepare Data

rfb_accessibility_plot_data <-
  rfb_accessibility_data |>
  mutate(
    accessibility_category = accessibility_per_1000 |> category_cut()
  ) |>
  relocate(h3_address) |>
  relocate(
    accessibility_category,
    .after = accessibility_per_1000
  )

Plot Accessibility

rfb_accessibility_plots <- list()

for (i in locais_nova_groups) {
  i_rfb_accessibility_data <-
    rfb_accessibility_plot_data |>
    filter(
      locais_nova_group |>
        equals(
          i |>
            str_replace("_", "-") |>
            str_to_upper()
        )
    )

  i_rfb_accessibility_plot <-
    rfb_accessibility_plot_template(coord_sf = FALSE) +
    geom_sf(
      data = i_rfb_accessibility_data,
      aes(fill = accessibility_category),
      color = NA
    ) +
    scale_fill_manual(
      values = fill_scale[[i]],
      name = NULL,
      guide = guide_legend(
        order = 3,
        ncol = 1,
        override.aes = list(color = "black", linewidth = 0.25)
      )
    ) +
    coord_sf(expand = FALSE)

  rfb_accessibility_plots <-
    rfb_accessibility_plots |>
    inset2(i, i_rfb_accessibility_plot)
}

Compute Plot Ratio

rfb_accessibility_plot_ratio <-
  rfb_accessibility_plots |>
  pluck("g0") |>
  get_plot_ratio()

Save Plots

for (i in locais_nova_groups) {
  plot_dir |>
    file.path(
      i |>
        str_replace_all("_", "-") %>%
        paste0("rfb-accessibility-", ., ".png")
    ) |>
    ggsave(
      plot = rfb_accessibility_plots[[i]],
      width = 21,
      height = (21 * rfb_accessibility_plot_ratio) - 1.75,
      units = "cm",
      dpi = 300,
      device = agg_png
    )
}

Outcomes

Analysis and Discussion for Physical Accessibility Outcomes

Citation

When using this data, you must also cite the original data sources.

To cite this work, please use the following format:

Pizzonia, M. A. Jr. (2026). Mapping Accessibility to Health Food in Brazilian Favelas and Urban Comunities [Computer software]. Center for Metropolitan Studies, University of São Paulo. https://cem-usp.github.io/kernel-accessibility/

A BibTeX entry for LaTeX users is

@software{pizzonia2026,
  title = {Mapping Healthy Food Accessibility in Brazilian Slums and Urban Communities},
  author = {{Marcio Antonio Pizzonia Junior}},
  year = {2026},
  address = {São Paulo},
  institution = {Center for Metropolitan Studies, University of São Paulo},
  langid = {en},
  url = {https://cem-usp.github.io/kernel-accessibility/}
}

License

The original data sources may be subject to their own licensing terms and conditions.

The code in this repository is licensed under the GNU General Public License Version 3, while the report is available under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International license.

Copyright (C) 2026 Center for Metropolitan Studies

The code in this report is free software: you can redistribute it and/or
modify it under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your option)
any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with
this program. If not, see <https://www.gnu.org/licenses/>.

Acknowledgments

This work is part of a research project by the Polytechnic School (Poli) of the University of São Paulo (USP), in partnership with the National Secretariat for Food and Nutrition Security (SESAN) of the Brazilian Ministry of Social Development and Assistance, Family and Fight Against Hunger (MDS), titled: AcessoSAN: Mapping Food Access to Support Public Policies on Food and Nutrition Security and Hunger Reduction in Brazilian Cities.

This work was developed with support from the Center for Metropolitan Studies (CEM) based at the School of Philosophy, Letters and Human Sciences (FFLCH) of the University of São Paulo (USP) and at the Brazilian Center for Analysis and Planning (CEBRAP).

This study was financed, in part, by the São Paulo Research Foundation (FAPESP), Brazil. Process Number 2025/17879-2.

References

Brasil, Ministério da Saúde, & Departamento de Informática do SUS. (n.d.). Banco de dados do Sistema Único de Saúde - DATASUS. Retrieved July 17, 2024, from https://datasus.saude.gov.br/

Brasil, & Ministério do Trabalho e Emprego. (2022). Base de dados da Relação Anual de Informações Sociais (RAIS) [Annual Social Information Report (RAIS) database] [Data set].

Brasil, & Sistema de Vigilância Alimentar e Nutricional. (n.d.). Módulo gerador de relatórios públicos: estado nutricional dos indivíduos acompanhados por período, fase do ciclo da vida e índice. Retrieved September 23, 2024, from https://sisaps.saude.gov.br/sisvan/relatoriopublico

Geofabrik GmbH. (n.d.). OpenStreetMap for Brazil [Data set]. https://download.geofabrik.de/south-america/brazil.html

Hollister, J., Shah, T., Nowosad, J., Robitaille, A. L., Beck, M. W., & Johnson, M. (2025). elevatr: Access elevation data from various APIs [Computer software]. https://doi.org/10.32614/CRAN.package.elevatr

Instituto Brasileiro de Geografia e Estatística, & Comissão Nacional de Classificação. (n.d.). CNAE 2.3 Subclasses [Data set]. SIDRA. Retrieved November 16, 2023, from https://concla.ibge.gov.br/classificacoes/download-concla.html

MapBiomas. (n.d.). Collection 2 (beta) of the land cover and land use annual maps of Brazil using Sentinel 2 images [Data set]. https://brasil.mapbiomas.org/en/mapbiomas-cobertura-10m

Mapzen. (n.d.). Terrain tiles [Data set]. https://registry.opendata.aws/terrain-tiles

O’Brien, L. (2023). h3jsr: Access Uber’s H3 library [Computer software]. https://doi.org/10.32614/CRAN.package.h3jsr

OpenStreetMap Contributors. (n.d.). OpenStreetMap [Computer software]. OpenStreetMap Foundation. https://www.openstreetmap.org

Pereira, R. H. M., & Goncalves, C. N. (n.d.). geobr: Download official spatial data sets of Brazil [Computer software]. https://doi.org/10.32614/CRAN.package.geobr

Pereira, R. H. M., & Herszenhut, D. (n.d.). geocodebr: Geolocalização de endereços brasileiros [geocodebr: Geolocation of Brazilian addresses] [Computer software]. https://doi.org/10.32614/CRAN.package.geocodebr

Pereira, R. H. M., Saraiva, M., Herszenhut, D., Braga, C. K. V., & Conway, M. W. (2021). r5r: Rapid Realistic Routing on Multimodal Transport Networks with R⁵ in R. Findings. https://doi.org/10.32866/001c.21262

R Core Team. (n.d.). R: A language and environment for statistical computing [Computer software]. R Foundation for Statistical Computing. https://www.r-project.org

Rocha, L. L., Lima Friche, A. A. de, Melo, G. B. V. de, Cordeiro, N. G., Honório, O. S., Oliveira Cardoso, L. de, & Mendes, L. L. (2024). Food retail in favelas of a brazilian metropolis. Food Security, 16(1), 277–292. https://doi.org/10.1007/s12571-023-01425-w

Uber Technologies. (n.d.). H3: A hexagonal hierarchical geospatial indexing system [Data set]. https://h3geo.org

Ushey, K., & Wickham, H. (2025). renv: Project environments [Computer software]. https://doi.org/10.32614/CRAN.package.renv

Walker, K. (2026). mapgl: Interactive Maps with ’Mapbox GL JS’ and ’MapLibre GL JS’ [Computer software]. https://walker-data.com/mapgl

Wickham, H. (n.d.-a). The tidyverse style guide. https://style.tidyverse.org

Wickham, H. (n.d.-b). Tidy design principles. https://design.tidyverse.org

Wickham, H. (2023). The tidy tools manifesto. Tidyverse. https://tidyverse.tidyverse.org/articles/manifesto.html

Wickham, H., Çetinkaya-Rundel, M., & Grolemund, G. (2023). R for data science: Import, tidy, transform, visualize, and model data (2nd ed.). O’Reilly Media. https://r4ds.hadley.nz

Overview

Problem

Data Availability

Methods

Source of Data

Data Munging

Travel Matrix

Code Style

Reproducibility

Set the Environment

Load Packages

Define Custom Functions

Set Keys

Set and Check Parameters

Set Data Directories

Download, Import, Process, and Plot Municipalities Data

Download and Save Data

Visualize Municipalities and Background Shapes

Municipalities Shapes

Import Data

Filter and Transform Data

Plot Data

Background Municipalities Shape

Import Data

Filter and Transform Data

Plot Background Shapes

Download, Import, Process, and Visualize Land Cover Shapes

Download Data

Import Data

Crop Data

Transform Data

Plot Data

Download, Import, Process, and Plot Favela’s and Urban Communitie’s Data

Download and Save Data

Import Data

Filter and Transform Data

Plot Data

Prepare Data

Layer on municipailities_data_shapes

Layer on background_municipalities_shapes

Layer on landcover_shapes

Set Plot Template

Define Template

Plot Template

Download, Import, and Visualize Establishments Data from Bazilian Federal Revenue (RFB)

Download Data

Unlock Files

Import Data

Transform and Filter Data

Visualize Data

Prepare Data

Build Base Maps

Combine Data and Maps

Compute Plot Ratio

Save Plots

Download, Import, and Visualize Establishmnts Data from Annual Social Information Report (RAIS)

Download Data

Unlock Files

Import Data

Transform and Filter Data

Visualize Data

Prepare Data

Build Base Maps

Combine Data and Maps

Compute Plot Ratio

Save Plots

Download, Import, and Visualize Accesibilitty Data from Bazilian Federal Revenue (RFB) and Annual Social Information Report (RAIS)

Annual Social Information Report (RAIS)

Download Data

Unlock Files

Import Data

Filter and Transform Data

Visualize Data

Define Plot Template

Prepare Data

Plot Accessibility

Compute Plot Ratio

Save Plots

Bazilian Federal Revenue (RFB)

Download Data

Layer on `municipailities_data_shapes`

Layer on `background_municipalities_shapes`

Layer on `landcover_shapes`