Entender o impacto regional de determinadas políticas públicas, de iniciativas privadas, de mudanças em indicadores de todo ordem, é determinante para a tomada de decisões adequadas no campo da gestão pública. Nestes contextos a manipulação de dados espaciais é uma ferramenta fundamental.

Imagine, por exemplo, que nos últimos seis meses observou-se aumento médio do uso da capacidade instalada da indústria do município de Piracicaba. Sabendo disso, queremos investigar o impacto regional deste crescimento no que diz respeito à absorção de mão de obra em termos regionais. A ideia seria analisar como a distância afeta a absorção de mão de obra e por isso vamos desagregar os resultados em de raios de 20km, 50km e 100km.

Neste caso é preciso identificar as cidades localizadas dentro destes raios. Este tutorial ensina a manipular dados espaciais no R com a finalidade de realizar tal tarefa.

Problema

Dada a escolha de uma cidade determinada, precisamos identificar quais são os municípos vizinhos sob seu raio de influência.

Preparar ambiente

Realizaremos a análise calculando os raios de duas formas diferentes. A primeira delas considerará os raios a partir das bordas do município de Piracicaba. A segunda considerará os raios a partir de um ponto específico da cidade, no caso o centroide, que representa o ponto central do território.

Pacotes utilizados

O procedimento que realizaremos para solucionar o problema depende do uso de alguns pacotes do R, sendo:

sf: estabelece classes e estruturas de dados espaciais, além de funções de manipulação. É uma forma aplicação da linguagem espacial no R que vem ganhando popularidade nos últimos anos, entre outras coisas, por incorporar a análise espacial na estrutura de sintaxe derivada do tidyverse.
dplyr: é um dos principais pacotes do tidyverse e um dos principais pacotes de manipulação de dados no R. Pode ser usado com dados espaciais de classe sf pois esta estrutura de dados é idêntica a um data.frame, diferentemente das classes de dados espaciais mais tradicional, derivadas do pacote sp.
ggplot2: também é um pacote do tidyverse para visualização de dados. É um dos mais usados e integra com o pacote sf.

Importar e preparar os dados

Para exemplificar, utilizaremos base de dados de municípos do estado de São Paulo do IBGE de 2017, a mais atual. Ela pode ser baixada diretamente do GeoFTP do IBGE.

# Baixar malha de municípios do estado de São Paulo
if (!file.exists(paste0("data/sp_municipios_2017"))){
  url_base <- "http://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/malhas_municipais/municipio_2017/UFs/SP/"
  arquivo <- "sp_municipios.zip"
  tmp <- tempfile(fileext = ".zip")
  download.file(paste0(url_base, arquivo), tmp)
  unzip(tmp, exdir = "data/sp_municipios_2017")
  unlink(tmp)
}

O código acima procura pela existência de um diretório chamado data/sp_municipios_2017. Se não existe, baixa o arquivo .zip do site do IBGE, o descompactado e depois armazena no diretório criado. O conjunto de arquivos descomptados são dados espaciais no formato shape (.shp).

Agora já é possível importar os dados no R, mas antes é preciso carregar os pacotes.

# carregar pacotes
library(dplyr)
library(sf)
library(ggplot2)

A função sf:read_sf() nos permite importar os dados de arquivos shape (.shp).

# importar dados
munic <- read_sf(dsn = "data/sp_municipios_2017",
                         layer = "35MUE250GC_SIR",
                         stringsAsFactors = FALSE)

Os dados foram importados e atribuídos a um objeto que chamamos munic. Agora já podemos checar os dados.

dim(munic)

## [1] 645   3

head(munic)

## Simple feature collection with 6 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -51.17924 ymin: -23.03648 xmax: -46.54881 ymax: -21.19701
## epsg (SRID):    4674
## proj4string:    +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
## # A tibble: 6 x 3
##   NM_MUNICIP       CD_GEOCMU                                       geometry
##   <chr>            <chr>                                 <MULTIPOLYGON [°]>
## 1 ADAMANTINA       3500105   (((-51.05425 -21.40465, -51.053 -21.40486, -5…
## 2 ADOLFO           3500204   (((-49.65795 -21.20333, -49.65645 -21.20549, …
## 3 AGUAÍ            3500303   (((-46.9764 -21.96818, -46.97599 -21.96832, -…
## 4 ÁGUAS DA PRATA   3500402   (((-46.73501 -21.81891, -46.73431 -21.81904, …
## 5 ÁGUAS DE LINDÓIA 3500501   (((-46.60614 -22.44173, -46.60347 -22.44307, …
## 6 ÁGUAS DE SANTA … 3500550   (((-49.19076 -22.72584, -49.19027 -22.72589, …

O dataset tem 645 linhas, sendo cada uma correspondente a um município do estado de São Paulo, e 2 variáveis tabulares e uma coluna com as coordenadas espaciais. NM_MUNICIP guarda o nome das cidades e CD_GEOCMU o código de identificação do município. A outra columa, geometry, guarda as coordenadas dos polígonos de cada cidade.

Podemos ver que a coordenada de referência dos dados está em +proj=longlat +ellps=GRS80 +no_defs, com as coordenadas em latitude e logitude. As funções de manipulação geométrica do sf assumem uma projeção diferente, baseada em DATUM

Método 1: raios a partir das bordas da cidade

Etapas de manipulação

A manipulação dos dados, efetivamente, envolve as seguintes etapas:

Criação dos raios de 10km, 50km e 100km a partir dos centroides
Identificação das cidades cujo território faz intersecção com cada um dos raios
Criação de variáveis binárias identificando quais cidades estão nos raios de cada uma das seis cidades (1 = a cidade está no raio, 0 = a cidade não está no raio)
Mescla das variáveis binárias na base de dados do excel

1. Criar os raios

A função sf::st_buffer() é a que usaremos para criar os raios. Ela pode ser aplicada tanto para pontos como para polígonos, e o que faz é criar um novo polígono baseado no argumento de distância indicado na função.

raio20_borda <- munic %>% 
  filter(NM_MUNICIP == "PIRACICABA") %>% 
  st_buffer(20000)

raio50_borda <- munic %>% 
  filter(NM_MUNICIP == "PIRACICABA") %>% 
  st_buffer(50000)

raio100_borda <- munic %>% 
  filter(NM_MUNICIP == "PIRACICABA") %>% 
  st_buffer(100000)

Vamos visualizar os raios:

ggplot() +
  geom_sf(data = munic %>% filter(NM_MUNICIP == "PIRACICABA")) +
  geom_sf(data = raio20_borda, alpha = 0, aes(col = "020 km")) +
  geom_sf(data = raio50_borda, alpha = 0, aes(col = "050 km")) +
  geom_sf(data = raio100_borda, alpha = 0, aes(col = "100 km"))

2. Identificar as cidades com interseção nos raios

Queremos saber quais as cidades vizinhas a Piracicaba que estão dentro dos limites dos raios de 20km, 50km e 100km. Isto é uma tarefa de intersecção de geometrias. Assim como a escolha de entre criar os raios a partir das bordas ou dos centroides é uma decisão a ser tomada, dependendo das características da análise a ser realizada, é preciso definir também o critério de intersecção. Ou seja, o que define a intersecção? Podemos considerar o polígono das cidades ou os seus centroides, por exemplo.

Neste caso, consideraremos sob os raios de influência qualquer cidade com ao menos uma parte do seu território (o polígono do seu território) dentro do raio.

Para realizar as interseções, usaremos em combinação as funções sf::st_join() e sf::st_intersects.

# raio de 20km
cidades_borda20 <- st_join(
  select(raio20_borda, NM_MUNICIP),
  select(munic, CD_GEOCMU),
  join = st_intersects
  ) %>% 
  mutate(raioborda20 = 1)

# raio de 50km
cidades_borda50 <- st_join(
  select(raio50_borda, NM_MUNICIP),
  select(munic, CD_GEOCMU ),
  join = st_intersects
  ) %>% 
  mutate(raioborda50 = 1)

# raio de 100km
cidades_borda100 <- st_join(
  select(raio100_borda, NM_MUNICIP),
  select(munic, CD_GEOCMU ),
  join = st_intersects
  ) %>% 
  mutate(raioborda100 = 1)

Fizemos a seleção das variáveis tabulares no procedimento para evitar redundância no resultado (se quiser entender, retire as seleções, rode o código e veja os resultados). Além disso, criamos uma variável chamada Raio, que identifica com o valor 1 os municípios extraídos. Esta variável será transformada na variável binária de identificação dos raios, na base munic.

3. Mescla na base de municípios do estado de São Paulo

munic <- munic %>%
  merge(cidades_borda20 %>% as.data.frame %>% select(-NM_MUNICIP, -geometry),
        by = "CD_GEOCMU", all.x = TRUE) %>% 
  merge(cidades_borda50 %>% as.data.frame %>% select(-NM_MUNICIP, -geometry),
        by = "CD_GEOCMU", all.x = TRUE) %>% 
  merge(cidades_borda100 %>% as.data.frame %>% select(-NM_MUNICIP, -geometry),
        by = "CD_GEOCMU", all.x = TRUE) %>% 
  mutate(raioborda20 = recode(raioborda20, .missing = 0),
         raioborda50 = recode(raioborda50, .missing = 0),
         raioborda100 = recode(raioborda100, .missing = 0))

Vamos ver agora os resultados:

head(munic %>% data.frame)

##   CD_GEOCMU             NM_MUNICIP raioborda20 raioborda50 raioborda100
## 1   3500105             ADAMANTINA           1           1            1
## 2   3500204                 ADOLFO           1           1            1
## 3   3500303                  AGUAÍ           1           1            1
## 4   3500402         ÁGUAS DA PRATA           1           1            1
## 5   3500501       ÁGUAS DE LINDÓIA           1           1            1
## 6   3500550 ÁGUAS DE SANTA BÁRBARA           1           1            1
##                         geometry
## 1 MULTIPOLYGON (((-51.05425 -...
## 2 MULTIPOLYGON (((-49.65795 -...
## 3 MULTIPOLYGON (((-46.9764 -2...
## 4 MULTIPOLYGON (((-46.73501 -...
## 5 MULTIPOLYGON (((-46.60614 -...
## 6 MULTIPOLYGON (((-49.19076 -...

Temos então as variáveis binárias que identificam com 0 e 1 se a cidade está ou não no respectivo raio.

Visualizar os resultados

Vamos visualizar agora o mapa do estado de São Paulo, com os raios de influência da cidade de Piracicaba. Para facilitar, vamos criar uma variável categórica, a partir das binárias, identificando os raios de influência.

munic <- munic %>% 
  mutate(
    Raio_borda = factor(if_else(NM_MUNICIP == "PIRACICABA", "Piracicaba",
                                if_else(raioborda20 == 1, "20km",
                                        if_else(raioborda50 == 1, "50km",
                                                if_else(raioborda100 == 1, "100km", "Fora dos raios")))),
                  levels = c("20km", "50km", "100km", "Fora dos raios", "Piracicaba"))
    )

Agora vamos criar o mapa:

p1 <- ggplot(munic) +
  geom_sf(aes(fill = Raio_borda), lwd = .1) +
  labs(title = "Raios de influência de Piracicaba",
       subtitle = "A partir das bordas") +
  seda::light_map_theme() +
  theme(legend.position = "bottom")
p1

Método 2: raios a partir das bordas da cidade

1. Gerar centroide

centroide <- munic %>% 
  filter(NM_MUNICIP == "PIRACICABA") %>% 
  st_centroid()

raio20_centroide <- centroide %>% 
  filter(NM_MUNICIP == "PIRACICABA") %>% 
  st_buffer(20000)

raio50_centroide <- centroide %>% 
  filter(NM_MUNICIP == "PIRACICABA") %>% 
  st_buffer(50000)

raio100_centroide <- centroide %>% 
  filter(NM_MUNICIP == "PIRACICABA") %>% 
  st_buffer(100000)

ggplot() +
  geom_sf(data = munic %>% filter(NM_MUNICIP == "PIRACICABA")) +
  geom_sf(data = centroide, size = 1) +
  geom_sf(data = raio20_centroide, alpha = 0, aes(col = "020 km")) +
  geom_sf(data = raio50_centroide, alpha = 0, aes(col = "050 km")) +
  geom_sf(data = raio100_centroide, alpha = 0, aes(col = "100 km"))