Mapas no R com pacote ‘sp’
Pontos, linhas, polígonos e classes sp
Em nossos exemplos utilizamos data frames para armazenar as infomações espaciais. Isso foi possível por que, ao trabalharmos com latitude e longitude de pontos, ficamos limitados a um único sistema de coordenadas e tampouco precisamos de mais de uma linha do data frame para armazenar o “desenho” de um objeto (ponto) em um mapa.
Há vários sistemas de coordenadas utilizados em diferentes bases de dados cartográficas. Assim, por um lado, precisamos de métodos que nos permitam transformar um sistema em outro para combinar diferentes fontes de dados. Dados especiais armazenados em data frames não contêm referência alguma ao sistema de coordenadas ao qual pertencem.
Por outro lado, se precisarmos de mapas que incluem outras representações geométricas, como linhas e pontos, somos obrigados a multiplicar o número de linhas em um data frame pelo número de vértices de tais figuras (2 vezes para uma linha e n vezes para um polígono, onde n é o número de vértices do polígono). Esse método, ainda que possível, é muito ineficiente.
A solução é trabalharmos com classes de objetos mais complexas (e completas) do que data frames. Vamos falar sobre isso ao explorarmos as principais classes do pacote sp. Em nosso tutorial vamos trabalhar com as classes de objetos do pacote sp e com as funções disponíveis para sua manipulação.
Antes disso, poré, vamos conhecer o formato mais comum para armazenar dados espaciais – shapefile – e como abrí-los em R.
Shapefiles
library(sp)
library(ggplot2)## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlanglibrary(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionAdemais de precisarmos de outras classes de objetos em R para armazenarmos informação espacial, precisamos também de formatos de arquivos diferentes para compartilhar dados espaciais. Shapefiles são os mais populares e diversos repositórios de dados espaciais, como a Prefeitura de São Paulo e o Centro de Estudos da Metrópole (CEM).
Nosso primeiro exemplo foi retirado do site do CEM e ligeiramente modificado. Sugiro que você dê uma olhada nas bases cartográficas do CEM, pois há bastante coisa interessante lá.
Vamos começar fazendo o download do arquivo .zip onde estão os dados:
url_shp_eleicoes_sp <- "https://github.com/leobarone/FLS6397/raw/master/data/rmsp.zip"
download.file(url_shp_eleicoes_sp, "temp.zip")
unzip("temp.zip")Veja que o arquivo descompactado é uma pasta (“rmsp”) com arquivos de diversas extensões. Não vamos entrar nos detalhes dos arquivos. Quando falamos de “shapefile”, estamos nos referindo a arquivos .shp. Mas raramente ele vem desacompanhado. A informação que nos importa é: todos eles, exceto, obviamente, o arquivo em .pdf, farão parte do objeto que criaremos.
list.files("rmsp")## [1] "MunRM07.dbf" "MunRM07.pdf" "MunRM07.shp" "MunRM07.shx"É fundamental que os arquivos de diferentes extensões tenham o mesmo nome, inclusive maiúsculas e minúsculas, e que os nomes das extensões estejam em letras minúsculas. Por esta razão, tome cuidado ao usar os arquivos do CEM e da PMSP, pois, em geral, algum arquivo vem com o nome e extensão em letras maiúscula.
A biblioteca para abertura de dados espaciais em R é rgdal. readORG, função de rgdal que utilizaremos para abrir os dados, recebe dois argumentos: a pasta em que estão os arquivos – no nosso caso “rmsp”, e o nome dos arquivos sem extensão (e por esta razão é fundamental que tenham o mesmo nome):
library(rgdal)## rgdal: version: 1.4-3, (SVN revision 828)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/jonny/OneDrive/Documents/R/win-library/3.6/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/jonny/OneDrive/Documents/R/win-library/3.6/rgdal/proj
## Linking to sp version: 1.3-1rmsp <- readOGR("rmsp", "MunRM07", stringsAsFactors = F)## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\jonny\Google Drive\Academic\USP\Teaching\FLS6397 - Intro to Programming\2019\FLS6397_2019\tutorials\rmsp", layer: "MunRM07"
## with 39 features
## It has 8 fields
## Integer64 fields read as strings: ID COD_IBGE POP_2000 DENS_DEMOVamos observar a classe do objeto importado:
class(rmsp)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"“SpatialPolygonsDataFrame” é, como é fácil deduzir, um objeto espacial que contém polígonos (municípios da Região Metropolitana de São Paulo) e que acompanha um data frame. Falaremos sobre essa classe de objetos a seguir.
Tal classe de objetos pode ser rapidamente visualizada utilizando o comando plot:
plot(rmsp)
Podemos utilizar o pacote ggplot2 e sua gramática para plotar objetos da classe “SpatialPolygonsDataFrame”:
library(mapproj)## Loading required package: mapsggplot(data = rmsp,
aes(x = long, y = lat, group = group)) +
geom_polygon(color = 'white', fill = 'darkgrey') +
coord_map()## Regions defined for each Polygons
Exercício:
Vá a uma das duas fontes de mapas indicadas – Prefeitura de São Paulo ou Centro de Estudos da Metrópole (CEM) – importe os dados e produza um mapa. Dependendo do tema que você escolher, você produzirá mapas com polígonos (por exemplo, mapas de divisões políticas ou administrativas), linhas (ruas, ciclovias, etc) ou pontos (unidades de saúde, escolas, etc).
Classes do pacote sp e data frames
Se nosso objetivo fosse apenas importar os dados e produzir um mapa, poderíamos para na sessão anterior. Contudo, a regra é tentarmos combinar bases de dados espaciais entre si ou com data frames que contém informações não presentes no arquivos que acompanham o “shapefile”. No Desafio 4, adicionaremos informações geradas no Desafio 3 ao mapa que acabamos de produzir com municípios da região metropolitana, por exemplo.
Como combinar dois “SpatialPolygonsDataFrame” com sistemas coordenadas diferentes? Como relacionar um data frame a um objeto da classe “SpatialPolygonsDataFrame”?
Vamos conhecer o que há dentro de um “SpatialPolygonsDataFrame” com a função str. Para não poluir o console, vamos adicionar o argumento “max.level = 2”:
str(rmsp, max.level = 2)## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 39 obs. of 8 variables:
## ..@ polygons :List of 39
## ..@ plotOrder : int [1:39] 12 37 10 5 11 1 16 6 3 4 ...
## ..@ bbox : num [1:2, 1:2] -47.2 -24.1 -45.7 -23.2
## .. ..- attr(*, "dimnames")=List of 2
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slotUm “SpatialPolygonsDataFrame” contém vários elementos. O principal deles são, obviamente, os polígonos. Para o objeto “rmsp”, são 39 polígonos. Já entraremos no detalhe de cada um. “plotOrder” e “bbox” são elementos do objeto que definem, respectivamente, a ordem de “plotagem” e as dimensões de um retângulo que contém todos os polígonos e raramente nos importaremos cm ambos.
“proj4string” armazena a informação sobre qual é sistema de coordenadas ao qual as coordenadas dos polígonos pertencem.
Um objeto da classe “SpatialPolygons” contém estes quatro elementos. Os objetos “SpatialPolygonsDataFrame” são iguais aos “SpatialPolygons” com a adição do elemento @data.
“SpatialLines”, “SpatialLinesDataFrame”, “SpatialMultiPoints” e “SpatialMultiPointsDataFrame” são as classes correspondentes para linhas e pontos e tem estrutura semelhante.
Podemos selecionar um elemento da estrutura de um objeto das classes do pacote sp usando o símbolo @. Vamos separar o objeto “rmsp” em 3 novos objetos:
rmsp_data <- rmsp@data
rmsp_poligonos <- rmsp@polygons
rmsp_projecao <- rmsp@proj4string“rmsp_data” é um data frame com os dados dos municípios:
head(rmsp_data)## ID COD_IBGE SIGLA NOME NOMECAPS POP_2000 DENS_DEMO
## 0 74 3546801 SIS Santa Isabel SANTA ISABEL 43740 120
## 1 76 3518305 GMA Guararema GUARAREMA 21904 81
## 2 75 3506607 BBM Biritiba-Mirim BIRITIBA MIRIM 24653 77
## 3 211 3518800 GRU Guarulhos GUARULHOS 1072717 3361
## 4 238 3545001 SPS Salesópolis SALESOPOLIS 14357 34
## 5 246 3528502 MRP Mairiporã MAIRIPORA 60111 187
## AREA_KM2
## 0 364.1687
## 1 271.8909
## 2 318.7029
## 3 319.0683
## 4 424.5088
## 5 321.4058“rmsp_projecao” é um
Conhecendo a estrutura de tais classes, já temos uma pista do que precisamos para (1) colocar duas projeções no mesmo sistema de coordenadas (alterando a informação sobre projeção) e (2) adicionar dados proveninentes de um data frame (combinando-o com o que estiverm em @data). Faremos este último a seguir e não cobriremos o primeiro. Antes disso, vamos entender a complexa estrutura dos elementos em @polygons.
Vamos criar o objeto “poligono1” com o primeiro elemento de “rmsp_poligonos” e examiná-lo.
poligono1 <- rmsp_poligonos[[1]]
str(poligono1)## Formal class 'Polygons' [package "sp"] with 5 slots
## ..@ Polygons :List of 1
## .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. ..@ labpt : num [1:2] -46.2 -23.3
## .. .. .. ..@ area : num 0.0321
## .. .. .. ..@ hole : logi FALSE
## .. .. .. ..@ ringDir: int 1
## .. .. .. ..@ coords : num [1:436, 1:2] -46.2 -46.2 -46.2 -46.2 -46.2 ...
## ..@ plotOrder: int 1
## ..@ labpt : num [1:2] -46.2 -23.3
## ..@ ID : chr "0"
## ..@ area : num 0.0321Note que ele é da classe “Polygons”, também do pacote sp. Um polígono contém várias informações: os vértices do polígono; as conexões entre os vértices; a existências de “buracos” no polígono (que o buraco é também um polígono e polígonos com buracos lembram rosquinhas); etc. Não vamos examinar o que há em cada polígono, mas você já deve ter percebido que “SpatialPolygons” e “SpatialPolygonsDataFrame” são conjuntos de “Polygons” com algumas informações adicionais (dentre as quais um data frame com características de cada polígono no caso do último). “SpatialPoints” e “Lines” são classes análogas.
A não se que você deseje criar novos (ou redesenhar) polígonos, não convém examinar essa classe de objetos.
Particionando um objeto de dados espaciais
Eventualmente, não são todos os polígonos importados de um “shapefile” que interessam. As classes do pacote sp aceitam seleção de polígonos da mesma forma que os data frames. No objeto “rmsp”, São Paulo é o 12o município, como vemos pelo vetor de nomes dos municípios em “rmsp_data”:
rmsp_data$NOME## [1] "Santa Isabel" "Guararema"
## [3] "Biritiba-Mirim" "Guarulhos"
## [5] "Salesópolis" "Mairiporã"
## [7] "Franco da Rocha" "Santana de Parnaíba"
## [9] "Pirapora do Bom Jesus" "Juquitiba"
## [11] "São Bernardo do Campo" "São Paulo"
## [13] "Cajamar" "Itapevi"
## [15] "Santo André" "Cotia"
## [17] "São Lourenço da Serra" "Osasco"
## [19] "Carapicuíba" "Jandira"
## [21] "Francisco Morato" "Rio Grande da Serra"
## [23] "Itaquaquecetuba" "Suzano"
## [25] "Mauá" "Itapecerica da Serra"
## [27] "Embu" "Taboão da Serra"
## [29] "Diadema" "São Caetano do Sul"
## [31] "Arujá" "Poá"
## [33] "Embu-Guaçu" "Barueri"
## [35] "Caieiras" "Ferraz de Vasconcelos"
## [37] "Moji das Cruzes" "Ribeirão Pires"
## [39] "Vargem Grande Paulista"Podemos separar o polígono de São Paulo e criar seu mapa da seguinte forma:
saopaulo <- rmsp[12,]
plot(saopaulo)
Repetindo o exemplo, mas para os 3 municípios do ABC:
posicoes_abc <- c(11, 15, 30)
abc <- rmsp[posicoes_abc,]
plot(abc)
Examinando a estrutura de “abc” notamos que o conjunto de polígonos e o data frame contém agora apenas 3 observações, como esperávamos.
str(abc, max.level = 2)## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 3 obs. of 8 variables:
## ..@ polygons :List of 3
## ..@ plotOrder : int [1:3] 1 2 3
## ..@ bbox : num [1:2, 1:2] -46.7 -24 -46.3 -23.6
## .. ..- attr(*, "dimnames")=List of 2
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slotAgregando e visualizando características de objetos espaciais
Vamos voltar ao nosso mapa da Região Metropolitana de São Paulo:
plot(rmsp)
Bastante sem graça, certo? Não estamos aprensentando nenhuma característica dos municipios no mapa, apenas os polígonos que os definem.
Os dados que acompanham o “shapefile”, porém, nos apresentam uma característica interessante dos munípios: a densidade demográfica (“DENS_DEMO”)
head(rmsp@data)## ID COD_IBGE SIGLA NOME NOMECAPS POP_2000 DENS_DEMO
## 0 74 3546801 SIS Santa Isabel SANTA ISABEL 43740 120
## 1 76 3518305 GMA Guararema GUARAREMA 21904 81
## 2 75 3506607 BBM Biritiba-Mirim BIRITIBA MIRIM 24653 77
## 3 211 3518800 GRU Guarulhos GUARULHOS 1072717 3361
## 4 238 3545001 SPS Salesópolis SALESOPOLIS 14357 34
## 5 246 3528502 MRP Mairiporã MAIRIPORA 60111 187
## AREA_KM2
## 0 364.1687
## 1 271.8909
## 2 318.7029
## 3 319.0683
## 4 424.5088
## 5 321.4058Vamos repetir o mapa, colorindo os municípios de acordo com a densidade demográfica:
plot(rmsp, col = log(as.numeric(rmsp@data$DENS_DEMO)))
Não muito bonito, mas bastante mais interessante. Vamos trazer uma variável externa aos dados para tornar nosso exemplo mais atraente. Utilizaremos os dados de planejamento urbano da MUNIC-IBGE 2015. Note que, por parcimônia, não estamos buscando os dados diretamente na fonte, mas em uma cópia dos dados no repositório do curso.
No nosso exemplo buscaremos o ano de elaboração do plano diretor de cada município da RMSP nas informações sobre planejamento urbano. Separaremos esta variável de interesse e o código do município, ambas variáveis jpa renomeadas:
url_munic_15 <- "https://raw.githubusercontent.com/leobarone/FLS6397/master/data/planejamento_munic_2015.csv"
munic_15 <- read.table(url_munic_15, header = T, sep = ";") %>%
rename(COD_IBGE = A1, ano_pd = A18) %>%
select(COD_IBGE, ano_pd) %>%
mutate(ano_pd = as.numeric(as.character(ano_pd)))## Warning: NAs introduced by coercionPara adicioanrmos essa informação ao objeto de dados espaciais, basta fazermos um left_join com os dados em @data:
rmsp@data <- rmsp@data %>%
mutate(ID = as.numeric(ID), COD_IBGE = as.numeric(COD_IBGE),
DENS_DEMO = as.numeric(DENS_DEMO)) %>%
left_join(munic_15, by = "COD_IBGE")Nada de novo, exceto o fato de que o data frame está dentro de um objeto mais complexo. Vamos agora plotar o ano de elaboração dos Planos Diretores de cada município como cores nos mapas:
plot(rmsp, col = factor(rmsp@data$ano_pd))
De novo, bastante feio, mas conseguimos trazer rapidamente dados de uma fonte externa para nosso mapa.
Deixando o mapa bonito com ggplot2
Vimos anterioremente que podemos plotar mapas usando o pacote ggplot2. Vejamos o mpa básico:
ggplot(data = rmsp,
aes(x = long, y = lat, group = group)) +
geom_polygon(color = 'white', fill = 'darkgrey') +
coord_map()## Regions defined for each Polygons
Para adicionarmos novas variáveis utilizando o ggplot2, porém, precisamos reorganizar nossos dados. Infelizmente, se simplemente adicionarmos uma variável que esteja em @data ao argumento “fill”, por exemplo, obteremos um mensagem de erro.
A solução será reorganizar os dados de polígonos em um data frame novo com a função fortify, perdendo justamente a eficiência que as classes de objetos do pacote sp forneceram.
Na sequência abaixo, criaremos “rmsp_df” com a função fortify. “rmsp_df” terá apenas variáveis “espaciais”. Para adicionar as demais variáveis que estão em @data, teremos de criar um “id” comum a ambos os data frames e realizar um left_join:
rmsp_df <- fortify(rmsp)## Regions defined for each Polygonsrmsp_df$id <- as.numeric(rmsp_df$id)
rmsp@data$id <- 0:(nrow(rmsp@data)-1)
rmsp_df <- left_join(rmsp_df, rmsp@data, by = "id")Observemos o objeto “rmsp_df”:
glimpse(rmsp_df)## Observations: 19,673
## Variables: 16
## $ long <dbl> -46.15734, -46.16225, -46.16271, -46.16266, -46.1622...
## $ lat <dbl> -23.34258, -23.34526, -23.34622, -23.34684, -23.3475...
## $ order <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1...
## $ hole <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FAL...
## $ piece <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ id <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ group <fct> 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0....
## $ ID <dbl> 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, ...
## $ COD_IBGE <dbl> 3546801, 3546801, 3546801, 3546801, 3546801, 3546801...
## $ SIGLA <chr> "SIS", "SIS", "SIS", "SIS", "SIS", "SIS", "SIS", "SI...
## $ NOME <chr> "Santa Isabel", "Santa Isabel", "Santa Isabel", "San...
## $ NOMECAPS <chr> "SANTA ISABEL", "SANTA ISABEL", "SANTA ISABEL", "SAN...
## $ POP_2000 <chr> "43740", "43740", "43740", "43740", "43740", "43740"...
## $ DENS_DEMO <dbl> 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 12...
## $ AREA_KM2 <dbl> 364.1687, 364.1687, 364.1687, 364.1687, 364.1687, 36...
## $ ano_pd <dbl> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007...Note que todas as variáveis que precisamos estão aí. Reorganizando os argumentos do nosso mapa anterior, produzido com ggplot, podemos gerar mapas bastante mais elegantes. O primeiro, usando a informação sobre densidade demográfica:
ggplot(data = rmsp_df,
aes(x = long, y = lat, group = group, fill = DENS_DEMO)) +
geom_polygon() +
coord_map()
E agora o ano dos planos diretores, variável que contém um número grande de missing values para os municípios da RMSP:
ggplot(data = rmsp_df,
aes(x = long, y = lat, group = group, fill = ano_pd)) +
geom_polygon() +
coord_map()
A grande vantagem de trabalharmos com a gramática do ggplot2 é que podemos editar os mapas da mesma forma que editamos gráficos e as possibilidades são inúmeras.