Visualisando mapas

Lembra que nossos gráficos em ggplot foram conectados diretamente com um data.frame? Como um objeto de ‘simple features’ é um data.frame, podemos usar ggplot para criar mapas! Para visualizar o nosso mapa, vamos usar uma camada de geometria especial (e espacial) se chama geom_sf, e o formato dos parâmetros são os mais simples possíveis: branco!

emef %>% ggplot() +
  geom_sf()

Já criou o seu primeiro mapa! O eixo x é o longitude e o eixo y é o latitude. Debaixo, vamos ver como a representar outros aspectos de nosso data.frame em mapas. Agora, tem um outro abordagem para criar um mapa interativo, com um mapa de terra no fundo. Usamos o pacote ‘mapview’ e, de novo, não precisamos especificar nenhum parâmetro.

library(mapview)
emef %>% mapview()

Georeferencing

Um outro fonte de dados espaciais é georeferenciar um endereço usando uma ferramenta de busca como o open street maps (mas fácil a acessar do que o google maps), por exemplo com a biblioteca tmaptools do pacote tmap.

Vamos agora criar um novo data frame a partir dos dados do cadastro da PMSP que contém apenas os CEUs. Vamos juntar algumas informações de endereço e agregar a elas o texto “, Sao Paulo, Brazil”

ceu <- escolas  %>%
  filter(tipo == "CEU") %>%
  mutate(endereco = paste(ENDERECO, NUMERO, ", Sao Paulo, Brazil"))

Com a função geocode_OSM, procuraremos a latitude e longitude dos 46 CEUs. Vamos ver o exemplo do primeiro CEU, criando um objeto tipo ‘sf’ com o nosso sistema de referência de coordenadas (WGS84 = epsg 4326) para ficar consistente:

library(tmaptools)
ceu_geocoded <- geocode_OSM(ceu$endereco, projection=4326, as.sf=T)

Podemos mapear os resultados:

ceu_geocoded %>% mapview()

Agora, podemos comparar as localizações no banco de dados original da PMSP com os nossos resultados de geocodificação, aplicando o nosso conhecimento do ggplot2.

ceu %>% ggplot() + 
  geom_sf(color="blue") +
  geom_sf(data=ceu_geocoded, color="red",
          alpha=0.5)

Simples, não! Alguns erros, mas a maioria dos conjuntos de coordenados parecem bem pertos. O principal problema da função geocode_OSM - e por que você pode ver menos de 46 pontos geocodificados em seu mapa - é que ele depende muito na qualidade de endereços. ### Exercício Use a função geocode_OSM para geolocalizar três endereços com os quais você está familiarizado em São Paulo e criar um mapa com ggplot e um mapa interativa com mapview.

Trabalhando com Polígonos

Áreas administrativas são geralmente representadas como polígonos em mapas. Em geral, obtemos esses polígonos como ‘shapefiles’ produzidos por uma agência oficial. Podemos abrir qualquer tipo de shapefile (pontos, linhas ou polígonos) com a função read_sf. Vamos abrir um shapefile (simplificado) de IBGE dos municipios do Brasil.

download.file("https://github.com/JonnyPhillips/Curso_R/raw/master/Brazil_s.zip",destfile="Brazil_s.zip")
unzip("Brazil_s.zip")
municipios <- read_sf("Brazil_s.shp")

Como podemos visualizar este mapa? Exatamente o mesmo que antes (pode demorar para abrir):

municipios %>% ggplot() +
  geom_sf()

Enquanto trabalhando com dados espaciais de diversas fontes (ex. as escolas e os municipios), é essencial ter certeza que estamos trabalhando com a mesma projeção para todos os objetos - caso contrário nossos mapas podem não aparecer corretamente e nossas medidas espaciais serão imprecisas. Dar uma olhada na CRS/projeção de ‘municipios’ - está escrito “+proj=longlat +ellps=GRS80 +no_defs” - em contraste de CRS/projeção de ‘escolas’ - “+proj=longlat +datum=WGS84 +no_defs”. O parte ‘WGS84’ é o mais comum então vamos transformar os dados dos municipios para a mesma projeção usando a função st_transform e o shortcut 4326 para WGS84. Agora estamos tudo pronto para fazer análises ou mapas incorporando as duas camadas espaciais.

municipios <- municipios %>% st_transform(4326)

Joins Não-Espaciais

Se abre o objeto municipios vai ver os detalhes de CRS/projeção e o data.frame dele mesmo, incluindo a coluna de ‘geometry’ e também todas as colunas normais, incluindo o código municipal (CD_GEOCODM). Isso é uma oportunidade para nós - se tivermos dados não-espaciais de todos os municípios, podemos simplesmente ‘join’ estes dados com o nosso shapefile e, em seguida, podemos visualizar mapas dessas variáveis.

Por exemplo, vamos baixar os dados eleitorais de 2010 para cada município com o percentagem de voto da Dilma no segundo turno.

download.file("https://github.com/JonnyPhillips/FLS6397_2019/raw/master/data/Dilma_2010_segundo_turno.csv",destfile="Dilma_2010_segundo_turno.csv")

Dilma_2010 <- read_csv("Dilma_2010_segundo_turno.csv")

Agora, se quisermos tornar esses dados eleitorais espaciais podemos fazer um left_join com o nosso shapefile, contanto que, como normal, os nomes e tipo de colunas chaves - o código do IBGE - nos dois bancos de dados são os mesmos.

municipios <- municipios %>% rename("COD_MUN_IBGE"="CD_GEOCODM") %>% 
  mutate(COD_MUN_IBGE=as.numeric(COD_MUN_IBGE)) %>%
  left_join(Dilma_2010,by="COD_MUN_IBGE")

Quantas colunas o nosso objeto municipios tem agora? Mais! - os dados eletorais também.

municipios

Para visualizar a coluna ‘pct_voto’ num mapa, podemos trabalhar em ggplot como normal. Para polígonos, colocamos o nome de coluna para o parâmetro fill ( color com pontos e linhas). Um mapa de todo o Brasil pode ser esmagador, então vamos nos concentrar em São Paulo com um filter.

municipios %>% filter(UF=="SP") %>% 
  ggplot() +
  geom_sf(aes(fill=pct_voto))

Vamos destacar as fronteiras em branco, aplicar um tema e tirar as linhas de longitude/latitude.

municipios %>% filter(UF=="SP") %>% 
  ggplot() +
  geom_sf(aes(fill=pct_voto),color="white") +
  theme_classic() +
  coord_sf(datum=NA)

Além disso, podemos alterar as escalas de cores como normal em ggplot.

municipios %>% filter(UF=="SP") %>% 
  ggplot() +
  geom_sf(aes(fill=pct_voto),color="white") +
  theme_classic() +
  coord_sf(datum=NA) +
  scale_fill_gradient(low="white",high="red")

Podemos colocar os pontos dos CEUs no mesmo mapa? Sim - apenas especificamos o parâmetro de ‘data’ para uma nova camada dentro de geom_sf. (Isto assume que nós confirmamos que ambas as camadas têm a mesma projeção, como confirmamos acima).

municipios %>% filter(UF=="SP") %>% 
  ggplot() +
  geom_sf(aes(fill=pct_voto),color="white") +
  theme_classic() +
  scale_fill_gradient(low="white",high="red") +
  geom_sf(data=ceu,color="dark green") +
  coord_sf(datum=NA)

Joins Espaciais

O mundo espacial abre um novo tipo de join entre diversas bancos de dados - joins espaciais que são definido pela localização semelhante e não por uma chave comum nas tabelas de dados. Existe diversas tipos de joins espaciais mas vamos focar sobre um join entre uma camada de polígonos e uma camada de pontos. Queremos pegar os dados de polígono (município) em que fica cada ponto (escola). Vamos usar st_join para fazer um join espacial e usar o tipo de comparação de st_contains (que simplesmente significa que um join existe quando um ponto é dentro de um polígono). (Lembra-se que é importante confirmar que ambas as camadas têm a mesma projeção, como confirmamos acima).

joined <- ceu %>% st_join(municipios,st_intersects)

Agora, o objeto ‘joined’ contém os dados de escolas e colunas adicionais com dados dos polígonos em que cada escola existe. Isso facilita comparações entre cada escola e os dados eleitorais da sua região, mas também podemos simplesmente contar o número de escolas em cada polígono (município). Claro que neste caso, a maioria das escolas estão localizadas na PMSP.

joined %>% group_by(NM_MUNICIP) %>% count()

Outras operações espaciais

O pacote simple features inclui muitas metodologias espaciais. Como exemplo, podemos calcular a distancia entre cada um das escolas CEU (em metros) com st_distance. O resultado é um matriz (simétrico) com distancias entre escolas.

ceu %>% st_distance()

Rasters

Existe um outro formato para dados espaciais que não é baseado em formas geométricas (polígonos, pontos e linhas), mas em uma grade regular com valores específicos em cada célula x, y - isto é um ‘raster’ e para trabalhar com ele usamos o pacote ‘raster’. Vamos usar o código debaixo para abrir um arquivo raster de densidade populacional no Camboja, que é simplesmente um imagem com extensão .tif.

library(raster)
download.file("https://github.com/JonnyPhillips/Curso_R/raw/master/khm_popdenr_landscan_2011.zip",destfile="khm_popdenr_landscan_2011.zip")
unzip("khm_popdenr_landscan_2011.zip")
cambodia <- raster("khm_popdenr_landscan_2011.tif")

Para visualizar o nosso raster, precisamos transformar ele em um data.frame simples e usar o ggplot com a geometria de geom_tile.

cambodia %>% as("SpatialPixelsDataFrame") %>% 
  as.data.frame() %>% 
  ggplot() + 
  geom_tile(aes(x=x,y=y,fill=khm_popdenr_landscan_2011))

Este mapa parece bem chato porque os dados são altamente ‘skewed’, com grandes outliers apenas na capital. Frequentemente com rasters é útil transformá-los em uma escala de log para visualizar. Vamos também limpar o fundo e adicionar uma escala de cores.

cambodia %>% as("SpatialPixelsDataFrame") %>% 
  as.data.frame() %>% 
  ggplot() + 
  geom_tile(aes(x=x,y=y,fill=khm_popdenr_landscan_2011)) +
  coord_equal() +
  theme_void() +
  scale_fill_gradient(low="white",high="red",na.value="white", trans="log")

Também podemos criar mapas de raster interativos com mapview. (Adicionamos uma pequena quantidade ao log para evitar valores infinitos e erros).

log(cambodia+0.00001) %>% 
  brick() %>% 
  mapview(layer.name="khm_popdenr_landscan_2011")