invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL","pt_BR.UTF-8"))
options(warn=-1)
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(kableExtra, warn.conflicts=FALSE))
suppressMessages(library(foreign, warn.conflicts=FALSE))
suppressMessages(library(geobr, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(miniCRAN, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
# suppressMessages(library(rgdal, warn.conflicts=FALSE))
suppressMessages(library(rnaturalearth, warn.conflicts=FALSE))
# se nao estiver disponivel para sua versão de R:
# remotes::install_github("ropensci/rnaturalearthhires")
# devtools::install_github("ropensci/rnaturalearthhires")
suppressMessages(library(rnaturalearthhires, warn.conflicts=FALSE))
suppressMessages(library(sf, warn.conflicts=FALSE))
suppressMessages(library(sp, warn.conflicts=FALSE))
suppressMessages(library(dplyr, warn.conflicts=FALSE))
suppressMessages(library(openxlsx, warn.conflicts=FALSE))
suppressMessages(library(KernSmooth, warn.conflicts=FALSE))
suppressMessages(library(foreign, warn.conflicts=FALSE))
source("eiras.viridis.hot.R")
source("eiras.friendlycolor.R")
options(warn=0)

Material

  • HTML de R Markdown em RPubs

Plano da aula

Alguns exemplos

Mapas podem modelar a elevação do terreno, como neste link e neste link.

Este vídeo tem um exemplo de construção com raster, usando ggplot, para mostrar a distribuição da chuva ao longo do ano: o código R é exposto à medida que cada comando vai sendo utilizado. Também é possível criar modelos de terrenos, como este mapa da Tanzânia — tem código fonte usando GMT, mas há componentes adicionais que precisam ser encontrados para sua execução. Detalhes em:

Lemenkova, Polina. 2022. Tanzania Craton, Serengeti Plain and Eastern Rift Valley: Mapping of Geospatial Data by Scripting Techniques. Estonian Journal of Earth Sciences 71 (2)..

GMT (Generic Mapping Tools) é uma coleção de ferramentas de linha de comando usada para criar mapas e realizar análises de dados geoespaciais, amplamente utilizado em geociências, geofísica, oceanografia e outras áreas que envolvem a manipulação e visualização de dados geográficos. GMT é um software open-source e pode ser utilizado em Linux, macOS e Windows.

Um script básico do GMT geralmente segue esta estrutura:

  • Inicialização: Configurações iniciais para o ambiente GMT, como definir o tamanho do mapa, a resolução e outras opções de configuração.
  • Comandos de plotagem: Uso de comandos GMT para adicionar elementos ao mapa, como contornos, linhas costeiras, rótulos, legendas, etc.
  • Finalização: Comandos para fechar o arquivo de saída e limpar qualquer configuração temporária.

Por exemplo, costa.gmt:

#!/bin/bash

# Configuração inicial do GMT
gmt begin simple_map png

# Configura o tamanho da região e o tipo de projeção
gmt basemap -R-180/180/-90/90 -JG0/0/6i -Baf

# Adiciona as linhas costeiras
gmt coast -W1p -Glightgray

# Finaliza e gera o arquivo PDF
gmt end show

gera

É possível monitorar o regime de ventos, exibindo em tempo real, como neste site — mas não verifiquei se tem código fonte disponível (este exemplo é mencionado em r-bloggers.com mas o link que poderia ensinar como criar os mapas vetoriais está quebrado.

Para quem tem interesse, possíveis formas estão em r-graphics.org ou milospopovic.net/ (eu não tentei implementar estas rotinas).

\[~\]

Uma versão antiga, que fizemos no início da pandemia COVID_cases_UF.mp4, plotava os casos de acordo com a localização dos municípios, o que também é interessante ver.

Observe como foi feita a formatação do eixo das abscissas, exibindo as datas em um formato humanamente mais legível.

Os dados continuam sendo atualizados pelo Governo Federal em https://covid.saude.gov.br/.

Procurando algo similar sobre Dengue, encontrei o Plano de Dados Abertos do Ministério da Saúde. Encontrei o “Conjunto de dados” (um arquivo csv) onde há uma lista do que está disponível.

Destacamos alguns momentos da pandemia:

  • 15 de abril de 2020:

  • 27 de maio de 2020:

  • 27 de julho de 2020:

  • 22 de março de 2021:

A sobreposição dos mapas com a localização de aeroportos, rodovias, ferrovias e hidrovias mostra o seguinte:

  • Aeroportos

  • Rodovias

  • Ferrovias

  • Hidrovias

Pacotes principais (implementação em R)

sp e sf

Há uma grande quantidade de pacotes com funções interrelacionadas para a manipulação de mapas e sua associação com procedimentos estatísticos. Os mais integradores são sp para desenhar objetos de classes espaciais (Pebesma & Bivand, 2018) e spatstat (Baddeley, 2021; Kopczewska, 2021). O pacote sp, que cuida dos mapas (Mas et al., 2019), trabalhava com dependência de:

  • rgdal para ler e gravar objetos vetoriais (OGR), que foi removido do CRAN
  • raster (GDAL) (Bivand et al., 2017a) para ler e escrever dados
  • rgeos (Bivand & Rundel, 2017) para operacionalização espacial, removido do CRAN
  • maptools (Becker & Wilks, 2018) que lidava com conversões, removido do CRAN
  • PBSmapping (Schnute et al., 2003) para manipular os objetos espaciais e permitir a conversão para uso em outros pacotes (ainda existe).

Conexões dos pacotes com base em sp e spatstat (baseado em Kopczewska, 2021).

Existe, ainda, o pacote sf (Pebesma, 2017), que promete simplificar, mantendo compatibilidade, as diversas operações de análise espacial, por exemplo integrando a exibição de resultados com leitura e gravação dos dados[1]. Algumas de suas operações mais simples são mostradas adiante.

Conexões entre pacotes com base em sf.

Implementado com demo_miniCRAN.R.

Estas duas figuras mostram os principais pacotes utilizados em análise espacial.

A análise espacial em R começou com grande profusão e confusão de função. Foi o trabalho de Roger Bivand que integrou pacotes, organizando classes de objetos e permitindo a conversão dos respectivos formatos quando necessário passar dos recursos de um para outro[3].

O objetivo deste texto é auxiliar nos passos iniciais. Neste texto mostraremos como iniciar com os pacotes sp e sf; instale-os para replicar os exemplos.

Funções:

  • sf::st_point(): cria ponto (objeto sfg, de geometria) com 2 a 4 dimensões
  • sf::st_multipoint(): cria vários pontos (pode receber as coordenadas em matriz)
  • sf::st_linestring(): cria linha a partir das coordenadas
  • sf::st_multilinestring(): cria objeto com várias linhas
  • sf::st_polygon(): cria polígono a partir das coordenadas
  • sf::st_multilinestring(): cria objeto com várias linhas
  • sf::st_sfc(): cria objeto sfc, contendo uma coleção de geometrias
  • sf::st_sf(): adiciona atributos a uma geometria criada
  • sf::st_geometry(): retorna somente a geometria de um objeto sfc
  • sf::st_coordinate(): retorna as coordenadas de um objeto sfc
  • sf::st_read(): lê um shapefile

geobr

Dados existem em muitas fontes. Para dados oficiais do Brasil, por exemplo, existe o pacote geobr, que contém:

suppressMessages(library(geobr, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(kableExtra, warn.conflicts=FALSE))

dt_br <- geobr::list_geobr()
# dt_br %>%
#   knitr::kable("html") %>%
#   kableExtra::kable_styling("striped", full_width = FALSE) %>%
#   kableExtra::scroll_box(width = "100%", height = "400px")

table_html <- knitr::kable(dt_br, format = "html")
styled_table <- kableExtra::kable_styling(table_html,
                                          bootstrap_options = "striped", 
                                          full_width = FALSE)
final_output <- kableExtra::scroll_box(styled_table, width = "100%", height = "400px")
print(final_output)
function geography years source
read_country Country 1872, 1900, 1911, 1920, 1933, 1940, 1950, 1960, 1970, 1980, 1991, 2000, 2001, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 IBGE
read_region Region 2000, 2001, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 IBGE
read_state States 1872, 1900, 1911, 1920, 1933, 1940, 1950, 1960, 1970, 1980, 1991, 2000, 2001, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 IBGE
read_meso_region Meso region 2000, 2001, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 IBGE
read_micro_region Micro region 2000, 2001, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 IBGE
read_intermediate_region Intermediate region 2017, 2019, 2020 IBGE
read_immediate_region Immediate region 2017, 2019, 2020 IBGE
read_municipality Municipality 1872, 1900, 1911, 1920, 1933, 1940, 1950, 1960, 1970, 1980, 1991, 2000, 2001, 2005, 2007, 2010, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022 IBGE
read_municipal_seat Municipality seats (sedes municipais) 1872, 1900, 1911, 1920, 1933, 1940, 1950, 1960, 1970, 1980, 1991, 2010 IBGE
read_weighting_area Census weighting area (área de ponderação) 2010 IBGE
read_census_tract Census tract (setor censitário) 2000, 2010, 2017, 2019, 2020 IBGE
read_statistical_grid Statistical Grid of 200 x 200 meters 2010 IBGE
read_metro_area Metropolitan areas 1970, 2001, 2002, 2003, 2005, 2010, 2013, 2014, 2015, 2016, 2017, 2018 IBGE
read_urban_area Urban footprints 2005, 2015 IBGE
read_amazon Brazil’s Legal Amazon 2012 MMA
read_biomes Biomes 2004, 2019 IBGE
read_conservation_units Environmental Conservation Units 201909 MMA
read_disaster_risk_area Disaster risk areas 2010 CEMADEN and IBGE
read_indigenous_land Indigenous lands 201907, 202103 FUNAI
read_semiarid Semi Arid region 2005, 2017, 2021, 2022 IBGE
read_health_facilities Health facilities 201505, 202303 CNES, DataSUS
read_health_region Health regions and macro regions 1991, 1994, 1997, 2001, 2005, 2013 DataSUS
read_neighborhood Neighborhood limits 2010 IBGE
read_schools Schools 2020, 2023 INEP
read_comparable_areas Historically comparable municipalities, aka Áreas mínimas comparáveis (AMCs) 1872, 1900, 1911, 1920, 1933, 1940, 1950, 1960, 1970, 1980, 1991, 2000, 2010 IBGE
read_urban_concentrations Urban concentration areas (concentrações urbanas) 2015 IBGE
read_pop_arrangements Population arrangements (arranjos populacionais) 2015 IBGE

Muitos dos mapas disponíveis usam o modelo geométrico de recursos simples (simple feature, em inglês), um padrão estabelecido pelo Open Geospatial Consortium (OGC), o qual utilizamos com os procedimentos em R que veremos adiante. Pela Web é possível encontrar muitos mapas de diversos locais do mundo pronto para seu uso.

ggplot2

  • ggplot2::ggplot: função básica de gráfico do pacote ggplot2
  • ggplot2::geom_sf: exibição de shapefile
  • ggplot2::geom_point: pontos em mapa com longitude e latitude.
  • ggplot2::aes: controle de atributos dos elementos de ggplot()
  • ggplot2::labs: rótulos em eixos ou legendas
  • ggplot2::ylim: cortes de segmentos dos mapas em latitudes
  • ggplot2::xlim: cortes de segmentos dos mapas em longitudes
  • ggplot2::theme_XXX & ggplot2::theme(parameters): para controlar a aparência dos mapas com temas (XXX indica que há vários a escolher)
  • ggplot2::element_blank, ggplot2::element_line, ggplot2::element_rect: usadas com os parâmetros dos temas.

Tipos de dados espaciais

Podem ser tomados em três situações[3] quando:

  • o valor é de interesse e sua localização é assinalada em um ponto.
  • a localização é de interesse e o padrão de distribuição dos pontos é assinalado (obviamente, os pontos podem ter cores ou intensidade de acordo com seus valores).
  • dados agregados são de interesse e sua localização em área (polígono) de ocorrência é assinalada (também utilizando cores ou intensidade para dar noção dos valores).

Dados vetoriais

Um dos formatos mais difundidos para análise espacial são mapas distribuídos em simple feature[1]. Este tipo de objeto hierárquico permite que várias formas geográficas sejam armazenadas conjuntamente.

O pacote sf tem funções para criar e armazenar estas formas geográficas ou manipular objetos sf já existentes.

Dados raster

Além das formas vetoriais (pontos, linhas e polígonos) de expressar localização espacial, há pacotes que dividem a área em pequenas unidades discretas, formando um grid que pode receber cores ou intensidades de preto e branco, conhecido como formato raster. Como este formato preenche áreas com um padrão, é frequentemente sobreposto aos mapas vetoriais.

Formas geográficas elementares

Usando sf construiremos, passo-a-passo, um objeto no formato simple feature. Embora o mapa resultante seja muito elementar, facilitará o entendimento da estrutura dos mapas que encontramos prontos para baixar pela Internet.

pontos

Os pontos podem ter de 2 a 4 dimensões em suas coordenadas. As duas primeiras coordenadas servem para posicionar o ponto no plano (x,y). A terceira é utilizada para representar altura, em mapas com informação sobre relevo. A quarta pode servir para algum outro parâmetro que se deseje estudar (densidade, temperatura, etc.).

Embora não pareça em nada diferente de um gráfico comum, criar os pontos com o pacote sf permitirá, adiante, integrar os pontos em objetos geográficos. Veja o resultado de demo_ponto2d.R:

cat("\nPonto, duas coordenadas:\n")
pto <- sf::st_point(c(3,5)) # 2 dimensões (XY)
print(pto)
cat("\nclass:\n")
print(class(pto))
cat("\nstr:\n")
print(str(pto))
# mapa
plot(pto, axes = TRUE)

Gerando o seguinte:


Ponto, duas coordenadas:
POINT (3 5)

class:
[1] "XY"    "POINT" "sfg"  

str:
 'XY' num [1:2] 3 5
NULL

A função sf::st_point() cria um ponto, que é um objeto de geometria (ou, para simplificar a linguagem, uma geometria).

Note que usamos as funções nativas do R, class(), str() e plot(). As duas primeiras mostram que o objeto criado não é um ponto qualquer, pois já está preparado para uso em mapa. Observe, também, o uso de sf:: antes da função st_point. Indicar o pacote é um hábito, explicitando a qual pacote a função pertence e evitando ambiguidades entre funções de mesmo nome; é útil especialmente quando o R Script for longo e revisitado muito tempo depois.

Usando plot(), neste exemplo fizemos com que os eixos fossem exibidos, o que costuma ser inconveniente para mapas. Interessantemente, a função plot() exibe os eixos por default, mas com objetos deste tipo não o faz. Com demo_ponto2d_2.R:

pto <- sf::st_point(c(3,5)) # 2 dimensões (XY)
plot(pto)

obtém-se:

Com pontos de duas coordenadas, mapas são bidimensionais. No entanto, outras dimensões podem ser armazenadas em um ponto, por exemplo, Com demo_ponto4d.R monta um ponto com 4 coordenadas:

cat("\nPonto, quatro coordenadas:\n")
pto <- sf::st_point(c(3,5,18,25)) # 4 dimensões (XYZM)
print(pto)
cat("\nclass:\n")
print(class(pto))
cat("\nstr:\n")
print(str(pto))
# mapa
plot(pto)

que contém, como podemos verificar:


Ponto, quatro coordenadas:
POINT ZM (3 5 18 25)

class:
[1] "XYZM"  "POINT" "sfg"  

str:
 'XYZM' num [1:4] 3 5 18 25
NULL

O ponto tem quatro valores. Embora não haja efeito visível nesta imagem, observe que as duas coordenadas adicionais foram incorporadas ao ponto (observe o que class() informa), para uso em outras situações. A terceira coordenada pode ser usada para indicar profundidade; a quarta pode conter algum outro dado de interesse.

multiponto

Vários pontos podem ser armazenados em uma geometria. Com demo_multiponto.R

x <- c(10,1,6,4)
y <- c(6,6,4,11)
z <- c(12,15,11,8)
pontos <- cbind(x,y,z) # cria matriz com tres colunas
cat("\nmatrix com os pontos:\n")
print(pontos)
cat("\nobjeto multipoint:\n")
mtpto <- sf::st_multipoint(pontos) 
print(mtpto)
cat("\nclass:\n")
print(class(mtpto))
cat("\nstr:\n")
print(str(mtpto))
# mapa
plot(mtpto)

obtém-se:


matrix com os pontos:
      x  y  z
[1,] 10  6 12
[2,]  1  6 15
[3,]  6  4 11
[4,]  4 11  8

objeto multipoint:
MULTIPOINT Z ((10 6 12), (1 6 15), (6 4 11), (4 11 8))

class:
[1] "XYZ"        "MULTIPOINT" "sfg"       

str:
 'XYZ' num [1:4, 1:3] 10 1 6 4 6 6 4 11 12 15 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:3] "x" "y" "z"
NULL

Note que class(), agora, informa que há três dimensões e que o objeto é MULTIPOINT.

coleção de pontos

Os objetos definidos acima podem ser agrupados em uma simple feature collection (sfc) usando a função st_sfc(), como em demo_geometria.R.

x <- 3
y <- 5
pto <- sf::st_point(c(x,y))
x <- c(10,1,6,4)
y <- c(6,6,4,11)
z <- c(12,15,11,8)
ptos <- cbind(x,y,z) # cria matriz com tres colunas
mtpto <- sf::st_multipoint(ptos) 
# objeto composto
cat("\nobjeto sfc (colecao de geometrias):\n")
geometria <- sf::st_sfc(pto,mtpto)
print(geometria)
cat("\nclass:\n")
print(class(geometria))
cat("\nstr:\n")
print(str(geometria))
# mapa
plot(geometria)

que gera:


objeto sfc (colecao de geometrias):
Geometry set for 2 features 
Geometry type: GEOMETRY
Dimension:     XY, XYZ
Bounding box:  xmin: 1 ymin: 4 xmax: 10 ymax: 11
z_range:       zmin: 8 zmax: 15
CRS:           NA
POINT (3 5)
MULTIPOINT Z ((10 6 12), (1 6 15), (6 4 11), (4...

class:
[1] "sfc_GEOMETRY" "sfc"         

str:
sfc_GEOMETRY of length 2; first list element:  'XY' num [1:2] 3 5
NULL

Observe a saída, que informa ser “Geometry set for 2 features”. Os dois componentes são POINT e MULTIPOINT. Podemos verificar como foram armazenados na geometria que acabamos de criar.

print(class(geometria[1]))
print(geometria[1])
[1] "sfc_POINT" "sfc"      
Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3 ymin: 5 xmax: 3 ymax: 5
CRS:           NA
POINT (3 5)
print(class(geometria[2]))
print(geometria[2])
[1] "sfc_MULTIPOINT" "sfc"           
Geometry set for 1 feature 
Geometry type: MULTIPOINT
Dimension:     XYZ
Bounding box:  xmin: 1 ymin: 4 xmax: 10 ymax: 11
z_range:       zmin: 8 zmax: 15
CRS:           NA
MULTIPOINT Z ((10 6 12), (1 6 15), (6 4 11), (4...

As duas geometrias que incorporamos com sf::st_sfc() preservam, respectivamente, suas classes como sfc_POINT e sfc_MULTIPOINT. O ponto isolado tem duas e os quatro multipontos têm três dimensões cada um.

adicionando atributos à geometria

Esta geometria pode ser sofisticada, incluindo atributos para cada ponto. Suponha que o primeiro ponto seja uma unidade de tratamento de água, e que os multipontos definidos sejam poços artesianos, executando-se demo_geometria2.R.

# cria objeto sfc
x <- c(3,10,1,6,4)
y <- c(5,6,6,4,11)
z <- c(20,12,15,11,8)
coords <- cbind(x,y,z) # cria matriz com tres colunas
ptos <- c(NA,rep(nrow(coords)))
for (i in 1:nrow(coords))
{
  ptos[i] <- list(sf::st_point(coords[i,]))
}
geometria2 <- sf::st_sfc(ptos)
# cria atributos para os pontos
id <- 1:length(geometria2)
tipo <- c("Unidade de tratamento",rep("Poço",4))
dt <- data.frame(id,tipo)
# monta novo objeto com os pontos
geometria2 <- sf::st_sf(dt,geometry=geometria2)
# exibe o resultado
cat("\nobjeto:\n")
print(class(geometria2))
print(geometria2)
cat("\nmapas:\n")
plot(geometria2)

objeto:
[1] "sf"         "data.frame"
Simple feature collection with 5 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 1 ymin: 4 xmax: 10 ymax: 11
z_range:       zmin: 8 zmax: 20
CRS:           NA
  id                  tipo          geometry
1  1 Unidade de tratamento  POINT Z (3 5 20)
2  2                  Poço POINT Z (10 6 12)
3  3                  Poço  POINT Z (1 6 15)
4  4                  Poço  POINT Z (6 4 11)
5  5                  Poço  POINT Z (4 11 8)

mapas:

A função sf::st_sf() criou uma simple feature de classe híbrida, contendo um sf e um data frame. Do dataframe reconhecem-se dois atributos, que são id (o número de identificação) e o tipo (unidade de tratamento ou poço).

O comportamento da função plot() mudou, exibindo um mapa para cada atributo (e usando uma cor diferente para cada um): no primeiro mapa há 5 números de identificação diferentes (com cinco cores); no segundo mapa representa apenas dois tipos (com duas cores).

Quando se deseja exibir somente a geometria usa-se:

plot(sf::st_geometry(geometria2))

Para recuperar a informação do dataframe:

print(as.data.frame(geometria2))
  id                  tipo          geometry
1  1 Unidade de tratamento  POINT Z (3 5 20)
2  2                  Poço POINT Z (10 6 12)
3  3                  Poço  POINT Z (1 6 15)
4  4                  Poço  POINT Z (6 4 11)
5  5                  Poço  POINT Z (4 11 8)

Para criar um objeto que é um subconjunto de outro objeto, por exemplo separando-se os pontos que representam poços:

pocos <- geometria2[tipo=="Poço",]
print(class(pocos))
[1] "sf"         "data.frame"
print(pocos)
Simple feature collection with 4 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 1 ymin: 4 xmax: 10 ymax: 11
z_range:       zmin: 8 zmax: 15
CRS:           NA
  id tipo          geometry
2  2 Poço POINT Z (10 6 12)
3  3 Poço  POINT Z (1 6 15)
4  4 Poço  POINT Z (6 4 11)
5  5 Poço  POINT Z (4 11 8)
plot(sf::st_geometry(pocos))

linhas (simples ou múltiplas)

Uma linha pode ser definida por um ou mais segmentos de reta posicionados pelas coordenadas (x,y) de suas extremidades. Suponha que os poços estejam ligados à unidade de tratamento de água por tubulações, representadas por linhas. Da mesma forma que foi feito com pontos, pode-se criar objetos com coleções de linhas. Por exemplo, executando-se demo_linhas.R (depende de rodar o exemplo anterior para criar a unidade de tratamento e os poços antes):

# recupera as coordenadas dos pontos
# tubulacoes comecam nos pocos
pocos <- geometria2[tipo=="Poço",]
coord_pocos <- sf::st_coordinates(pocos$geometry)
# tubulacoes terminam na unidade de tratamento
tubos <- geometria2[tipo=="Unidade de tratamento",]
coord_trata <- sf::st_coordinates(tubos$geometry)
# maxmin y
ymin <- min(coord_pocos[,2],coord_trata[,2])
ymax <- max(coord_pocos[,2],coord_trata[,2])

# criando uma linha de cada poco ate a un. de tratamento
x <- c(coord_pocos[1,1], 6.3, coord_trata[1,1])
y <- c(coord_pocos[1,2], 6.1, coord_trata[1,2])
xy <- cbind(x,y)
l1 <- sf::st_linestring(xy)
x <- c(coord_pocos[2,1],2.2, 2.7, coord_trata[1,1])
y <- c(coord_pocos[2,2],6.6, 4.8, coord_trata[1,2])
xy <- cbind(x,y)
l2 <- sf::st_linestring(xy)
x <- c(coord_pocos[3,1],5.3, 3.8, coord_trata[1,1])
y <- c(coord_pocos[3,2],4.1, 4.6, coord_trata[1,2])
xy <- cbind(x,y)
l3 <- sf::st_linestring(xy)
x <- c(coord_pocos[4,1],3.3, 3.8, coord_trata[1,1])
y <- c(coord_pocos[4,2],10.2, 8.1, coord_trata[1,2])
xy <- cbind(x,y)
l4 <- sf::st_linestring(xy)
# junta as linhas
adutoras <- sf::st_multilinestring(list(l1,l2,l3,l4))

# mapa
plot(adutoras,lwd=2,col="blue", xlim=c(0,12), ylim=c(0,12))
plot(pocos$geometry,pch=21,bg="lightblue",add=TRUE)
plot(tubos$geometry,pch=21,bg="green",cex=1.5,add=TRUE)

temos:

polígonos (áreas)

Polígonos fechados também podem ser representados em objetos. O cuidado é fornecer as coordenadas \((x,y)\) de forma que o último par de coordenadas coincida com o primeiro para fechar o polígono. Por exemplo, demo_poligonos.R:

library("sf")

## P1 Polígono rio
## Polygon
# Cria uma cadeia de coordenadas em X
x <- c(0,6,10,12,12,11,5,0,0)
y <- c(0,0,1,3,6,5,3,3,0)
c <- cbind(x,y)
# Cria um objeto polygon (forma simples fechada)
rio = sf::st_polygon(list(c))
# Lencol freatico
x <- c(0,5,11,12,12,10,6,3,2,1,0,0)
y <- c(5,3,5,6,7,10,12,12,10,9,7,5)
c <- cbind(x,y)
# cria um buraco, area seca
x <- c(2,3,4,6,5,3,3,2)
y <- c(5,4,4,7,9,9,7,5)
c2 <- cbind(x,y)
lencol = sf::st_polygon(list(c,c2))

# floresta alta
x <- c(0,5,11,12,12,0,0)
y <- c(3,3,5,6,12,12,3)
c <- cbind(x,y)
floresta1 = sf::st_polygon(list(c))
# floresta baixa
x <- c(6,12,12,10,6)
y <- c(0,0,3,1,0)
c <- cbind(x,y)
floresta2 = sf::st_polygon(list(c))

# Junta várias geometrias com sfc (coleção de simple features)
geometria.all = sf::st_sfc(rio, lencol, floresta1, floresta2)

# Tabela de atributos
ID <- c(1,2,3,4)
code <- c("A","L","F","F")
tipo <- c("Agua","Lencol freatico", "Floresta", "Floresta")
tabela <- data.frame(cbind(ID,code,tipo))
# sf object
sfc_all = sf::st_sf(tabela, geometry = geometria.all)
plot(sfc_all$geometry)
plot(st_geometry(sfc_all[tipo=="Agua",]),col="#00ffff",add=TRUE)
plot(st_geometry(sfc_all[tipo=="Floresta",]),col="#00880050",add=TRUE)
plot(st_geometry(sfc_all[tipo=="Lencol freatico",]),col="#ffff00",add=TRUE)

mostra:

Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

agrupando geometrias

Podemos exibir conjuntamente as duas coleções de geometrias que construímos. A primeira coleção contém os poços (em azul claro) e tubulações (em azul escuro) que levam água à unidade de tratamento (verde); a segunda monta um ambiente com um rio (em azul claro), lençol freático (em amarelo) e a mata (em verde, parcialmente transparente, sobrepondo o lençol freático subterrâneo).

plot(sfc_all$geometry, xlim=c(0,12), ylim=c(0,12))
plot(sf::st_geometry(sfc_all[tipo=="Agua",]),col="#00ffff",add=TRUE)
plot(sf::st_geometry(sfc_all[tipo=="Lencol freatico",]),col="#ffff00",add=TRUE)
plot(sf::st_geometry(sfc_all[tipo=="Floresta",]),col="#00880050",add=TRUE)
plot(adutoras,lwd=2,col="blue",add=TRUE)
plot(pocos$geometry,pch=21,bg="lightblue",add=TRUE)
plot(tubos$geometry,pch=21,bg="green",cex=1.5,add=TRUE)

Este é o princípo de vários mapas disponíveis na Web, como por exemplo:

Usando mapas prontos

Os arquivos com este tipo de objeto, contendo geometrias e atributos, são oferecidos em shapefiles, conjuntos formados, no mínimo, por arquivos:

  • .shp (com as geometrias),
  • .shx (um índice) e
  • .dbf (com os atributos).

Este padrão foi estabelecido pela ESRI e amplamente utilizado para a troca de informações geográficas. A ESRI (Environmental Systems Research Institute) tem sede mundial em Redlands na Califórnia, Estados Unidos, produzindo soluções para o uso de informações geográficas. Dependendo do distribuidor, há arquivos com dados adicionais.

Há vários sites que distribuem mapas gratuitos. Busque “shapefiles” e o tipo de mapa que quer na Web.

Mapa do mundo

Por exemplo, o arquivo

  • ne_50m_admin_0_countries.zip

disponível no site Natural Earth tem, ao descompactar, os dados de geometria em arquivo .shp, trazendo os polígonos de cada país. Este site oferece muitos mapas diferentes, com três níveis de detalhamento (o número 50 que aparece no nome do arquivo é o detalhamento moderado - veja Natural Earth, Features{target”blank”}).

Descompactando o arquivo em uma pasta chamada NaturalEarth, executa-se, no estilo sf:

sf_countries <- sf::st_read("NaturalEarth/ne_50m_admin_0_countries.shp")
plot(sf::st_geometry(sf_countries))
Reading layer `ne_50m_admin_0_countries' from data source 
  `D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Geoprocessamento\NaturalEarth\ne_50m_admin_0_countries.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 241 features and 94 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
Geodetic CRS:  WGS 84
Implementado com demo_sf_worldmap1.R.

Equivalentemente, para o modo sp mais tradicional:

sp_countries <- rgdal::readOGR(dsn="NaturalEarth",layer="ne_50m_admin_0_countries")
plot(sp_countries)
Implementado com demo_sp_worldmap1.R.

Embora o resultado visual seja o mesmo, uma diferença importante entre os dois métodos é que rgdal::readOGR separa o caminho (dns=) e o nome de base do objeto espacial sem a extensão .shp (layer=).

Os objetos sf_countries e sp_countries também são diferentes:

cat("\nsf_countries:\n")
print(paste0(class(sf_countries)))
# cat("\nsp_countries:\n")
# print(class(sp_countries))

sf_countries:
[1] "sf"         "data.frame"

Caso os verifique com str() ou print() verá que os dados estão internamente organizados de forma diversa.

Note, ainda, que as funções sf::st_read e rgdal::readOGR ecoam informações sobre o mapa que leram, o que pode não ser conveniente em sua aplicação. Na documentação das funções encontram-se os parâmetros que as silenciam, respectivamente:

  • quiet=TRUE para o pacote sf:
sfc_countries <- sf::st_read("NaturalEarth/ne_50m_admin_0_countries.shp", quiet=TRUE)
plot(sf::st_geometry(sfc_countries))
Implementado com demo_sf_worldmap1d1.R.
  • e verbose=FALSE para o pacote sp:
sp_countries <- rgdal::readOGR(dsn="NaturalEarth",layer="ne_50m_admin_0_countries", verbose=FALSE)
plot(sp_countries)
Implementado com demo_sp_worldmap1d1.R.

Como são mapas feitos com polígonos, podem ser preenchidos com cores, usando os valores prontos em R ou qualquer cor que deseje em RGB, com a importante diferença de que a função sf::st_geometry() é necessária para o primeiro estilo. Por exemplo,

plot(sf::st_geometry(sf_countries),col="lightgreen")
Implementado com demo_sf_worldmap1b.R.

ou

plot(sp_countries,col="lightgreen")
Implementado com demo_sp_worldmap1b.R.

encontrando o data frame

Embora os dois objetos tenham organizações internas diferentes, o acesso aos dados é similar. Além dos polígonos, para descobrir o que existe em seu data frame, use

print(names(as.data.frame(sf_countries)))
 [1] "featurecla" "scalerank"  "LABELRANK"  "SOVEREIGNT" "SOV_A3"    
 [6] "ADM0_DIF"   "LEVEL"      "TYPE"       "ADMIN"      "ADM0_A3"   
[11] "GEOU_DIF"   "GEOUNIT"    "GU_A3"      "SU_DIF"     "SUBUNIT"   
[16] "SU_A3"      "BRK_DIFF"   "NAME"       "NAME_LONG"  "BRK_A3"    
[21] "BRK_NAME"   "BRK_GROUP"  "ABBREV"     "POSTAL"     "FORMAL_EN" 
[26] "FORMAL_FR"  "NAME_CIAWF" "NOTE_ADM0"  "NOTE_BRK"   "NAME_SORT" 
[31] "NAME_ALT"   "MAPCOLOR7"  "MAPCOLOR8"  "MAPCOLOR9"  "MAPCOLOR13"
[36] "POP_EST"    "POP_RANK"   "GDP_MD_EST" "POP_YEAR"   "LASTCENSUS"
[41] "GDP_YEAR"   "ECONOMY"    "INCOME_GRP" "WIKIPEDIA"  "FIPS_10_"  
[46] "ISO_A2"     "ISO_A3"     "ISO_A3_EH"  "ISO_N3"     "UN_A3"     
[51] "WB_A2"      "WB_A3"      "WOE_ID"     "WOE_ID_EH"  "WOE_NOTE"  
[56] "ADM0_A3_IS" "ADM0_A3_US" "ADM0_A3_UN" "ADM0_A3_WB" "CONTINENT" 
[61] "REGION_UN"  "SUBREGION"  "REGION_WB"  "NAME_LEN"   "LONG_LEN"  
[66] "ABBREV_LEN" "TINY"       "HOMEPART"   "MIN_ZOOM"   "MIN_LABEL" 
[71] "MAX_LABEL"  "NE_ID"      "WIKIDATAID" "NAME_AR"    "NAME_BN"   
[76] "NAME_DE"    "NAME_EN"    "NAME_ES"    "NAME_FR"    "NAME_EL"   
[81] "NAME_HI"    "NAME_HU"    "NAME_ID"    "NAME_IT"    "NAME_JA"   
[86] "NAME_KO"    "NAME_NL"    "NAME_PL"    "NAME_PT"    "NAME_RU"   
[91] "NAME_SV"    "NAME_TR"    "NAME_VI"    "NAME_ZH"    "geometry"

ou

# print(names(as.data.frame(sp_countries)))

Repare que as.data.frame() funciona da mesma maneira para ambos.

Como este arquivo foi bem construído, os nomes das colunas são sugestivos de seu conteúdo. Que tal tentar a coluna NAME_PT para verificar?

# as.data.frame(sp_countries) fornece com o mesmo resultado
dt_countries <- as.data.frame(sf_countries) 
cat("\nTotal de ",nrow(dt_countries)," países disponíveis:\n",sep="")
print(dt_countries$NAME_PT)

Total de 241 países disponíveis:
  [1] "Zimbábue"                                 
  [2] "Zâmbia"                                   
  [3] "Iémen"                                    
  [4] "Vietname"                                 
  [5] "Venezuela"                                
  [6] "Vaticano"                                 
  [7] "Vanuatu"                                  
  [8] "Usbequistão"                              
  [9] "Uruguai"                                  
 [10] "Micronésia"                               
 [11] "Ilhas Marshall"                           
 [12] "Marianas Setentrionais"                   
 [13] "Ilhas Virgens Americanas"                 
 [14] "Guam"                                     
 [15] "Samoa Americana"                          
 [16] "Porto Rico"                               
 [17] "Estados Unidos"                           
 [18] "Ilhas Geórgia do Sul e Sandwich do Sul"   
 [19] "Território Britânico do Oceano Índico"    
 [20] "Santa Helena"                             
 [21] "Ilhas Pitcairn"                           
 [22] "Anguilla"                                 
 [23] "Ilhas Malvinas"                           
 [24] "Ilhas Cayman"                             
 [25] "Bermudas"                                 
 [26] "Ilhas Virgens Britânicas"                 
 [27] "Turks e Caicos"                           
 [28] "Montserrat"                               
 [29] "Jersey"                                   
 [30] "Guernsey"                                 
 [31] "Ilha de Man"                              
 [32] "Reino Unido"                              
 [33] "Emirados Árabes Unidos"                   
 [34] "Ucrânia"                                  
 [35] "Uganda"                                   
 [36] "Turquemenistão"                           
 [37] "Turquia"                                  
 [38] "Tunísia"                                  
 [39] "Trinidad e Tobago"                        
 [40] "Tonga"                                    
 [41] "Togo"                                     
 [42] "Timor-Leste"                              
 [43] "Tailândia"                                
 [44] "Tanzânia"                                 
 [45] "Tajiquistão"                              
 [46] "Taiwan"                                   
 [47] "Síria"                                    
 [48] "Suíça"                                    
 [49] "Suécia"                                   
 [50] "eSwatini"                                 
 [51] "Suriname"                                 
 [52] "Sudão do Sul"                             
 [53] "Sudão"                                    
 [54] "Sri Lanka"                                
 [55] "Espanha"                                  
 [56] "Coreia do Sul"                            
 [57] "África do Sul"                            
 [58] "Somália"                                  
 [59] "Somalilândia"                             
 [60] "Ilhas Salomão"                            
 [61] "Eslováquia"                               
 [62] "Eslovénia"                                
 [63] "Singapura"                                
 [64] "Serra Leoa"                               
 [65] "Seychelles"                               
 [66] "Sérvia"                                   
 [67] "Senegal"                                  
 [68] "Arábia Saudita"                           
 [69] "São Tomé e Príncipe"                      
 [70] "San Marino"                               
 [71] "Samoa"                                    
 [72] "São Vicente e Granadinas"                 
 [73] "Santa Lúcia"                              
 [74] "São Cristóvão e Nevis"                    
 [75] "Ruanda"                                   
 [76] "Rússia"                                   
 [77] "Roménia"                                  
 [78] "Catar"                                    
 [79] "Portugal"                                 
 [80] "Polónia"                                  
 [81] "Filipinas"                                
 [82] "Peru"                                     
 [83] "Paraguai"                                 
 [84] "Papua-Nova Guiné"                         
 [85] "Panamá"                                   
 [86] "Palau"                                    
 [87] "Paquistão"                                
 [88] "Omã"                                      
 [89] "Noruega"                                  
 [90] "Coreia do Norte"                          
 [91] "Nigéria"                                  
 [92] "Níger"                                    
 [93] "Nicarágua"                                
 [94] "Nova Zelândia"                            
 [95] "Niue"                                     
 [96] "Ilhas Cook"                               
 [97] "Países Baixos"                            
 [98] "Aruba"                                    
 [99] "Curaçao"                                  
[100] "Nepal"                                    
[101] "Nauru"                                    
[102] "Namíbia"                                  
[103] "Moçambique"                               
[104] "Marrocos"                                 
[105] "Saara Ocidental"                          
[106] "Montenegro"                               
[107] "Mongólia"                                 
[108] "Moldávia"                                 
[109] "Mónaco"                                   
[110] "México"                                   
[111] "Maurícia"                                 
[112] "Mauritânia"                               
[113] "Malta"                                    
[114] "Mali"                                     
[115] "Maldivas"                                 
[116] "Malásia"                                  
[117] "Malawi"                                   
[118] "Madagáscar"                               
[119] "República da Macedónia"                   
[120] "Luxemburgo"                               
[121] "Lituânia"                                 
[122] "Liechtenstein"                            
[123] "Líbia"                                    
[124] "Libéria"                                  
[125] "Lesoto"                                   
[126] "Líbano"                                   
[127] "Letónia"                                  
[128] "Laos"                                     
[129] "Quirguistão"                              
[130] "Kuwait"                                   
[131] "Kosovo"                                   
[132] "Kiribati"                                 
[133] "Quénia"                                   
[134] "Cazaquistão"                              
[135] "Jordânia"                                 
[136] "Japão"                                    
[137] "Jamaica"                                  
[138] "Itália"                                   
[139] "Israel"                                   
[140] "Palestina"                                
[141] "República da Irlanda"                     
[142] "Iraque"                                   
[143] "Irão"                                     
[144] "Indonésia"                                
[145] "Índia"                                    
[146] "Islândia"                                 
[147] "Hungria"                                  
[148] "Honduras"                                 
[149] "Haiti"                                    
[150] "Guiana"                                   
[151] "Guiné-Bissau"                             
[152] "Guiné"                                    
[153] "Guatemala"                                
[154] "Granada"                                  
[155] "Grécia"                                   
[156] "Gana"                                     
[157] "Alemanha"                                 
[158] "Geórgia"                                  
[159] "Gâmbia"                                   
[160] "Gabão"                                    
[161] "França"                                   
[162] "Saint-Pierre e Miquelon"                  
[163] "Wallis e Futuna"                          
[164] "São Martinho"                             
[165] "Coletividade de São Bartolomeu"           
[166] "Polinésia Francesa"                       
[167] "Nova Caledónia"                           
[168] "Terras Austrais e Antárticas Francesas"   
[169] "Åland"                                    
[170] "Finlândia"                                
[171] "Fiji"                                     
[172] "Etiópia"                                  
[173] "Estónia"                                  
[174] "Eritreia"                                 
[175] "Guiné Equatorial"                         
[176] "El Salvador"                              
[177] "Egito"                                    
[178] "Equador"                                  
[179] "República Dominicana"                     
[180] "Dominica"                                 
[181] "Djibouti"                                 
[182] "Gronelândia"                              
[183] "Ilhas Feroe"                              
[184] "Dinamarca"                                
[185] "República Checa"                          
[186] "República Turca de Chipre do Norte"       
[187] "Chipre"                                   
[188] "Cuba"                                     
[189] "Croácia"                                  
[190] "Costa do Marfim"                          
[191] "Costa Rica"                               
[192] "República Democrática do Congo"           
[193] "Congo"                                    
[194] "Comores"                                  
[195] "Colômbia"                                 
[196] "China"                                    
[197] "Macau"                                    
[198] "Hong Kong"                                
[199] "Chile"                                    
[200] "Chade"                                    
[201] "República Centro-Africana"                
[202] "Cabo Verde"                               
[203] "Canadá"                                   
[204] "Camarões"                                 
[205] "Camboja"                                  
[206] "Mianmar"                                  
[207] "Burundi"                                  
[208] "Burkina Faso"                             
[209] "Bulgária"                                 
[210] "Brunei"                                   
[211] "Brasil"                                   
[212] "Botsuana"                                 
[213] "Bósnia e Herzegovina"                     
[214] "Bolívia"                                  
[215] "Butão"                                    
[216] "Benim"                                    
[217] "Belize"                                   
[218] "Bélgica"                                  
[219] "Bielorrússia"                             
[220] "Barbados"                                 
[221] "Bangladesh"                               
[222] "Bahrein"                                  
[223] "Bahamas"                                  
[224] "Azerbaijão"                               
[225] "Áustria"                                  
[226] "Austrália"                                
[227] "Territórios australianos do Oceano Índico"
[228] "Ilha Heard e Ilhas McDonald"              
[229] "Ilha Norfolk"                             
[230] "Ilhas Ashmore e Cartier"                  
[231] "Arménia"                                  
[232] "Argentina"                                
[233] "Antígua e Barbuda"                        
[234] "Angola"                                   
[235] "Andorra"                                  
[236] "Argélia"                                  
[237] "Albânia"                                  
[238] "Afeganistão"                              
[239] "Glaciar de Siachen"                       
[240] "Antártida"                                
[241] "São Martinho"
Implementado com demo_worldmap1b1.R.

Para evitar redundâncias, alguns dos exemplos abaixo são feitos apenas em estilo sf ou sp com resultados visuais idênticos.

É fácil, caso deseje, adaptar o código para o outro estilo.

destacando um dos países (polígono)

Podemos separar o polígono que contém o Brasil e exibi-lo em outra cor (amarelo, neste exemplo):

sf_countries <- sf::st_read("NaturalEarth/ne_50m_admin_0_countries.shp", quiet=TRUE)
plot(st_geometry(sf_countries$geometry),col="lightgray")

dt_countries <- as.data.frame(sf_countries)
linha <- which(dt_countries$NAME_PT=="Brasil")
plot(sf_countries$geometry[linha],col="yellow",add=TRUE)
Implementado com demo_worldmap1c.R.

exibindo somente um país

Obviamente, também é possível exibir apenas um polígono. Embora isto não seja muito interessante por enquanto, mostra que o uso destas geometrias é muito flexível:

sp_countries <- rgdal::readOGR(dsn="NaturalEarth",layer="ne_50m_admin_0_countries", verbose=FALSE)
dt_countries <- as.data.frame(sfc_countries)
linha <- which(dt_countries$NAME_PT=="Brasil")
plot(sp_countries[linha,],col="yellow")
Implementado com demo_worldmap1d.R.

Observe as duas formas, nos dois estilos, para isolar um único polígono.

Em objetos lidos por de sf::st_read(), em vez de usar a função sf::geometry(), acessamos com

sf_countries$geometry[linha]

Para objetos lidos por rgdal::readOGR(), o acesso foi possível com

sp_countries[linha,]

(observe a vírgula após linha)

mapas feitos com linhas

Ainda do mesmo site há arquivos com outras possibilidades:

  • ne_50m_admin_0_boundary_lines_land.zip
  • ne_50m_coastline.zip

Estes dois arquivos fornecem as as fronteiras dos países e as linhas costeiras. Propositalmente exibimos aqui, primeiro, apenas as fronteiras:

sfc_borders <- sf::st_read("NaturalEarth/ne_50m_admin_0_boundary_lines_land.shp", quiet=TRUE)
plot(sf::st_geometry(sfc_borders$geometry))
Implementado com demo_worldmap1e.R.

… exatamente o que promete o arquivo: são linhas, e não polígonos como no caso anterior.

Para ter a aparência do mapa mundial, combine os dois (para clareza, exibimos com duas cores diferentes: preto para os contornos costeiros e cinza para as fronteiras entre os países):

sfc_coastline <- sf::st_read("NaturalEarth/ne_50m_coastline.shp", quiet=TRUE)
sfc_borders <- sf::st_read("NaturalEarth/ne_50m_admin_0_boundary_lines_land.shp", quiet=TRUE)
plot(sf::st_geometry(sfc_coastline$geometry),col="black")
plot(sf::st_geometry(sfc_borders$geometry),col="darkgray",add=TRUE)
Implementado com demo_worldmap1e1.R.

Note o uso de add=TRUE para sobrepor as duas geometrias.

combinando mapas

Interessantemente, mapas podem ser combinados. Neste exemplo, como usaremos o polígono que representa o Brasil, note que o parâmetro col indica a cor do preenchimento e border é a cor da linha; para os exemplos anteriores, que as geometrias eram linhas, col era a cor da linha. Além disto, como estamos usando a função nativa plot(), podemos aproveitar lwd para definir a espessura da linha de contorno do polígono (lwd, do inglês, line width).

Executa-se este exemplo, que destaca o Brasil e a Groenlândia:

sfc_coastline <- sf::st_read("NaturalEarth/ne_50m_coastline.shp", quiet=TRUE)
sfc_borders <- sf::st_read("NaturalEarth/ne_50m_admin_0_boundary_lines_land.shp", quiet=TRUE)
sfc_countries <- sf::st_read("NaturalEarth/ne_50m_admin_0_countries.shp", quiet=TRUE)

plot(sf::st_geometry(sfc_coastline$geometry),col="black")
plot(sf::st_geometry(sfc_borders$geometry),col="darkgray",add=TRUE)

# Brasil em destaque
dt_countries <- as.data.frame(sfc_countries)
linha <- which(dt_countries$NAME_PT=="Brasil")
plot(sf::st_geometry(sfc_countries$geometry[linha]),
     col="yellow",border="darkgreen",
     lwd=2, add=TRUE)

# Groenlandia em destaque
# por erro, aparece como "Gronelândia" no mapa
dt_countries <- as.data.frame(sfc_countries)
linha <- which(dt_countries$NAME_PT=="Gronelândia")
plot(sf::st_geometry(sfc_countries$geometry[linha]),
     col="#a6da9a",border="black",
     add=TRUE)
Implementado com demo_worldmap1e2.R.

Este mapa, pela aparência, usa uma projeção equiretangular (https://en.wikipedia.org/wiki/Equirectangular_projection), um tipo de projeção cilíndrica, das quais a mais famosa é a projeção de Mercator (https://en.wikipedia.org/wiki/Mercator_projection). O problema a ser resolvido pelas projeções é representar em um plano a superfície esférica.

Nas projeções cilíndricas, o principal defeito é distorcer as áreas dos países, ampliando o que está mais próximo dos polos. Por isso a Europa parece maior do que é, a Groenlândia, muito próxima ao polo norte, tem tamanho próximo ao Brasil (quando na verdade é cerca de 4 vezes menor), e a Antártida, no polo sul, é exibida com um tamanho tremendamente exagerado.

Uma animação disponível em https://commons.wikimedia.org/wiki/File:Worlds_animate.gif ilustra bem a distorção das áreas:

Cada projeção tem vantagens e desvantagens. Veja https://en.wikipedia.org/wiki/Map_projection para mais informações.

Existem formas de usar outros tipos de projeção na análise espacial em R. Por exemplo, podemos exibir o mapa mundial com projeção de (Mollweide)[https://en.wikipedia.org/wiki/Mollweide_projection], a qual preserva a proporção entre as áreas, com

sp_countries <- rgdal::readOGR(dsn="NaturalEarth",layer="ne_50m_admin_0_countries", verbose=FALSE)
plot(sp_countries)

sp2 <- sp::spTransform(sp_countries, sp::CRS("+proj=moll"))
plot(sp2)
dt_countries <- as.data.frame(sp_countries)
linha <- which(dt_countries$NAME_PT=="Brasil")
br2 <- sp::spTransform(sp_countries[linha,], sp::CRS("+proj=moll"))
plot(br2, col="yellow", add=TRUE)
Implementado com demo_Mollweide.R.

CRS (Coordinate Reference Systems) é o que governa as projeções. Uma documentação que pode ajudar está em (Overview of Coordinate Reference Systems (CRS) in R)[https://www.nceas.ucsb.edu/sites/default/files/2020-04/OverviewCoordinateReferenceSystems.pdf].

combinando mapas de diferentes fontes

Outra coisa interessante é que, sendo o formato padronizado, mapas de diferentes fontes podem ser combinados. Por exemplo, podemos pegar do website Nereus(USP) o arquivo Brasil.zip, que foi descompactado na pasta BrasilUSP, para exibir o contorno das unidades da federação:

sfc_coastline <- sf::st_read("NaturalEarth/ne_50m_coastline.shp", quiet=TRUE)
sfc_borders <- sf::st_read("NaturalEarth/ne_50m_admin_0_boundary_lines_land.shp", quiet=TRUE)
sfc_countries <- sf::st_read("NaturalEarth/ne_50m_admin_0_countries.shp", quiet=TRUE)

plot(sf::st_geometry(sfc_coastline$geometry),col="black")
plot(sf::st_geometry(sfc_borders$geometry),col="gray",add=TRUE)

# Brasil em destaque
dt_countries <- as.data.frame(sfc_countries)
linha <- which(dt_countries$NAME_PT=="Brasil")
plot(sf::st_geometry(sfc_countries$geometry[linha]),
     col="yellow",border="darkgreen",
     lwd=1, add=TRUE)

# Mapa dos estados
sfc_br <- sf::st_read("BrasilUSP/UFEBRASIL.shp",quiet=TRUE)
plot(sfc_br$geometry,add=TRUE,lwd=0.6)
Implementado com demo_worldmap1e3.R.

… ou, pegando o arquivo South_America.rar (disponível em Efrain Maps) que descompactamos na subpasta EfrainMaps, exibir somente a América do Sul em lugar do mundo inteiro:

# Paises do Mundo
sfc_countries <- sf::st_read("NaturalEarth/ne_50m_admin_0_countries.shp")

# America do Sul
sfc_americasul <- sf::st_read("EfrainMaps/South_America.shp", quiet=TRUE)
plot(sf::st_geometry(sfc_americasul))

# Brasil em destaque
dt_countries <- as.data.frame(sfc_countries)
linha <- which(dt_countries$NAME_PT=="Brasil")
plot(sf::st_geometry(sfc_countries$geometry[linha]),
     col="yellow",border="darkgreen",
     lwd=1, add=TRUE)

# Mapa dos estados
sfc_br <- sf::st_read("BrasilUSP/UFEBRASIL.shp",quiet=TRUE)
plot(sfc_br$geometry,add=TRUE,lwd=0.6)
Reading layer `ne_50m_admin_0_countries' from data source 
  `D:\Usuarios\Jose\Documents\FMUSP\Pos-graduacao\MDR5730\Geoprocessamento\NaturalEarth\ne_50m_admin_0_countries.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 241 features and 94 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
Geodetic CRS:  WGS 84
Implementado com demo_americasul1.R.

Podemos adicionar quantas camadas desejarmos, como por exemplo exibindo as hidrovias brasileiras (https://www.gov.br/infraestrutura/pt-br/assuntos/dados-de-transportes/bit/bitmodosmapas). Neste exemplo alteramos cores e espessura de linhas para dar destaque aos elementos de maior interesse:

library(sf)

# Paises do Mundo
sfc_countries <- sf::st_read("NaturalEarth/ne_50m_admin_0_countries.shp", quiet=TRUE)

# Paises da America do Sul
sfc_americasul <- sf::st_read("EfrainMaps/South_America.shp", quiet=TRUE)
plot(sf::st_geometry(sfc_americasul), border="gray")

# Brasil em destaque
dt_countries <- as.data.frame(sfc_countries)
linha <- which(dt_countries$NAME_PT=="Brasil")
plot(sf::st_geometry(sfc_countries$geometry[linha]),
     col="gray",border="black", lwd=1, add=TRUE)

# Mapa dos estados
sfc_br <- sf::st_read("BrasilUSP/UFEBRASIL.shp",quiet=TRUE)
plot(sfc_br$geometry,border="black",add=TRUE,lwd=0.3)

# Hidrovias
sfc_hidro <- sf::st_read("BrasilGov/Hidrovias.shp",quiet=TRUE)
plot(sfc_hidro$geometry,add=TRUE,col="darkcyan",lwd=1.5)
Implementado com demo_americasul2.R.

escrevendo no mapa

O mesmo mapa pode adicionar informações do data frame associado, por exemplo extraindo-se os nomes dos países da América do Sul e escrevendo-os perto dos centróides dos respectivos polígonos:

# Nomes dos países
paises <- as.data.frame(sfc_americasul)$COUNTRY
for (p in 1:length(paises))
{
  centroide <- sf::st_centroid(sfc_americasul$geometry[[p]])
  text(centroide[1],centroide[2],paises[p],col="brown",cex=0.8)
}
Implementado com demo_americasul2nome.R.

(algum refinamento poderia ser feito, acertando a posição das Guianas e encurtando o nome das ilhas oceânicas ao sul, mas isto vai além da necessidade deste texto)

outras funções do pacote sf

Podemos utilizar a função st_centroid() para definir onde escrever os nomes dos países. Há um grande número de funções neste pacote (e há vários outros pacotes que interagem com ele). Por exemplo:

  • st_distance(): computa a distância entre geometrias
  • st_intersects(): verifica se há intersecção (por exemplo, para saber se determinados pontos estão dentro de um polígono)
  • st_intersection(): que cria um novo objeto com a intersecção de outros.

Por exemplo, para desenhar a parte das hidrovias que estão no Estado do Amazonas, podemos adicionar:

# parte da Hidrovia que está no Estado do Amazonas
linha_UF <- which(sfc_br$NM_ESTADO=="AMAZONAS")
interseccao <- sf::st_intersection(sfc_br$geometry[linha_UF],sfc_hidro)
plot(interseccao,add=TRUE,lwd=2,col="purple")
Implementado com demo_americasul2partes.R.

Os métodos para encontrar os centróides nos objetos lidos por rgdal::readOGR() e funções para encontrar distâncias ou intersecções são diferentes. No entanto, o uso de plot() com objetos de um ou outro tipo misturam-se sem problema no mesmo mapa.

Usando o ggplot

O pacote ggplot2 traz a poderosa função ggplot para a criação de gráficos. Para este exemplo utilizaremos shapefiles dos pacotes rnaturalearth, parte do projeto ROpenSci. Para isto, a seguinte instalação precisa ser feita:

  install.packages("devtools")
  install.packages("rnaturalearth")
  install.packages("rnaturalearthdata")
  devtools::install_github("ropensci/rnaturalearthhires")

Há, nestes pacotes, objetos do tipo sf contendo mapas de várias regiões do mundo. Para separar os shapefiles que deseja, precisará encontrar o nome dos países ou regiões que lhe interessam. Neste exemplo quero exibir o mapa da América do Sul, com o Brasil e, em destaque, o Estado de São Paulo. Para localizá-los, usaremos as funções rnaturalearth::ne_countries e rnaturalearth::ne_states:

df <- rnaturalearth::ne_countries()
continentes <- unique(df$continent)
prmatrix(continentes,collab="",rowlab=rep("",length(continentes)))
                          
 "Oceania"                
 "Africa"                 
 "North America"          
 "Asia"                   
 "South America"          
 "Europe"                 
 "Seven seas (open ocean)"
 "Antarctica"             
prmatrix(df$name_long,collab="",rowlab=rep("",nrow(df)))
                                      
 "Fiji"                               
 "Tanzania"                           
 "Western Sahara"                     
 "Canada"                             
 "United States"                      
 "Kazakhstan"                         
 "Uzbekistan"                         
 "Papua New Guinea"                   
 "Indonesia"                          
 "Argentina"                          
 "Chile"                              
 "Democratic Republic of the Congo"   
 "Somalia"                            
 "Kenya"                              
 "Sudan"                              
 "Chad"                               
 "Haiti"                              
 "Dominican Republic"                 
 "Russian Federation"                 
 "Bahamas"                            
 "Falkland Islands / Malvinas"        
 "Norway"                             
 "Greenland"                          
 "French Southern and Antarctic Lands"
 "Timor-Leste"                        
 "South Africa"                       
 "Lesotho"                            
 "Mexico"                             
 "Uruguay"                            
 "Brazil"                             
 "Bolivia"                            
 "Peru"                               
 "Colombia"                           
 "Panama"                             
 "Costa Rica"                         
 "Nicaragua"                          
 "Honduras"                           
 "El Salvador"                        
 "Guatemala"                          
 "Belize"                             
 "Venezuela"                          
 "Guyana"                             
 "Suriname"                           
 "France"                             
 "Ecuador"                            
 "Puerto Rico"                        
 "Jamaica"                            
 "Cuba"                               
 "Zimbabwe"                           
 "Botswana"                           
 "Namibia"                            
 "Senegal"                            
 "Mali"                               
 "Mauritania"                         
 "Benin"                              
 "Niger"                              
 "Nigeria"                            
 "Cameroon"                           
 "Togo"                               
 "Ghana"                              
 "Côte d'Ivoire"                      
 "Guinea"                             
 "Guinea-Bissau"                      
 "Liberia"                            
 "Sierra Leone"                       
 "Burkina Faso"                       
 "Central African Republic"           
 "Republic of the Congo"              
 "Gabon"                              
 "Equatorial Guinea"                  
 "Zambia"                             
 "Malawi"                             
 "Mozambique"                         
 "Kingdom of eSwatini"                
 "Angola"                             
 "Burundi"                            
 "Israel"                             
 "Lebanon"                            
 "Madagascar"                         
 "Palestine"                          
 "The Gambia"                         
 "Tunisia"                            
 "Algeria"                            
 "Jordan"                             
 "United Arab Emirates"               
 "Qatar"                              
 "Kuwait"                             
 "Iraq"                               
 "Oman"                               
 "Vanuatu"                            
 "Cambodia"                           
 "Thailand"                           
 "Lao PDR"                            
 "Myanmar"                            
 "Vietnam"                            
 "Dem. Rep. Korea"                    
 "Republic of Korea"                  
 "Mongolia"                           
 "India"                              
 "Bangladesh"                         
 "Bhutan"                             
 "Nepal"                              
 "Pakistan"                           
 "Afghanistan"                        
 "Tajikistan"                         
 "Kyrgyzstan"                         
 "Turkmenistan"                       
 "Iran"                               
 "Syria"                              
 "Armenia"                            
 "Sweden"                             
 "Belarus"                            
 "Ukraine"                            
 "Poland"                             
 "Austria"                            
 "Hungary"                            
 "Moldova"                            
 "Romania"                            
 "Lithuania"                          
 "Latvia"                             
 "Estonia"                            
 "Germany"                            
 "Bulgaria"                           
 "Greece"                             
 "Turkey"                             
 "Albania"                            
 "Croatia"                            
 "Switzerland"                        
 "Luxembourg"                         
 "Belgium"                            
 "Netherlands"                        
 "Portugal"                           
 "Spain"                              
 "Ireland"                            
 "New Caledonia"                      
 "Solomon Islands"                    
 "New Zealand"                        
 "Australia"                          
 "Sri Lanka"                          
 "China"                              
 "Taiwan"                             
 "Italy"                              
 "Denmark"                            
 "United Kingdom"                     
 "Iceland"                            
 "Azerbaijan"                         
 "Georgia"                            
 "Philippines"                        
 "Malaysia"                           
 "Brunei Darussalam"                  
 "Slovenia"                           
 "Finland"                            
 "Slovakia"                           
 "Czech Republic"                     
 "Eritrea"                            
 "Japan"                              
 "Paraguay"                           
 "Yemen"                              
 "Saudi Arabia"                       
 "Antarctica"                         
 "Northern Cyprus"                    
 "Cyprus"                             
 "Morocco"                            
 "Egypt"                              
 "Libya"                              
 "Ethiopia"                           
 "Djibouti"                           
 "Somaliland"                         
 "Uganda"                             
 "Rwanda"                             
 "Bosnia and Herzegovina"             
 "North Macedonia"                    
 "Serbia"                             
 "Montenegro"                         
 "Kosovo"                             
 "Trinidad and Tobago"                
 "South Sudan"

Como existe Brazil, seus estados estão em

uf <- rnaturalearth::ne_states(country="Brazil")
prmatrix(uf$name_pt,collab="",rowlab=rep("",nrow(uf)))
                      
 "Rio Grande do Sul"  
 "Roraima"            
 "Pará"               
 "Acre"               
 "Amapá"              
 "Mato Grosso do Sul" 
 "Paraná"             
 "Santa Catarina"     
 "Amazonas"           
 "Rondônia"           
 "Mato Grosso"        
 "Maranhão"           
 "Piauí"              
 "Ceará"              
 "Rio Grande do Norte"
 "Paraíba"            
 "Pernambuco"         
 "Alagoas"            
 "Sergipe"            
 "Bahia"              
 "Espírito Santo"     
 "Rio de Janeiro"     
 "São Paulo"          
 "Goiás"              
 "Federal"            
 "Minas Gerais"       
 "Tocantins"

Mapa do Estado de São Paulo

Começando por São Paulo, podemos ver seu mapa, que é adicionado a ggplot2::ggplot com a função ggplot2::geom_sf (o pacote ggplot2 tem diversas funções especialmente desenvolvidas para mapas):

Brazil <- rnaturalearth::ne_states(country="Brazil", returnclass="sf")
SP <- Brazil[Brazil$name_pt=="São Paulo",]

grf <- ggplot2::ggplot()+
  ggplot2::geom_sf(data=SP, fill="gray", color="black")
print(grf)
Implementado com demo_ggplot_SaoPaulo01.R.

Adicionando municípios

Este mapa é um tanto sem graça. Vamos, então, detalhar o Estado em seus municípios. Como esta informação não existe no rnaturalearth, podemos procurar shapes de outras fontes. O (IBGE)[https://www.ibge.gov.br/geociencias/organizacao-do-territorio/malhas-territoriais/15774-malhas.html], por exemplo, as tem. Baixei o arquivo que nos interessa e coloquei na pasta IBGE (um nível abaixo de onde está o R Script). Executamos, então:

Brazil <- rnaturalearth::ne_states(country="Brazil", returnclass="sf")
SP <- Brazil[Brazil$name_pt=="São Paulo",]

municipios <- sf::st_read("IBGE/SP_Municipios_2021.shp", quiet=TRUE)
grf <- ggplot2::ggplot()+
  ggplot2::geom_sf(data=SP, fill="#CAE0AB", color="black")+
  ggplot2::geom_sf(data=municipios, fill="transparent", color="#507052", size=0.2)
print(grf)
Implementado com demo_ggplot_SaoPaulo02.R.

Observe que o ggplot já plotou latitude e longitude nos eixos. Além disto, como a construção do gráfico é feita em camadas, precisamos encadeá-las com + na ordem em que as queremos sobrepor. A função ggplot2::geom_sf aceita vários parâmetros: neste exemplo destacamos o contorno do estado com uma linha de cor preta e fundo (fill) em um tom de verde claro (RGB “#CAE0AB”); os municípios sem cor, têm apenas seus contornos em verde escuto (RGB “#507052”).

Destacando municípios

Suponha, porém, que eu queira separar apenas os municípios mais populosos. Segundo a Wikipedia há nove municípios com mais que 500.000 habitantes:

  • São Paulo, 12.396.372
  • Guarulhos, 1.404.694
  • Campinas, 1.223.237
  • São Bernardo do Campo, 849.874
  • São José dos Campos, 737.310
  • Santo André, 723.889
  • Ribeirão Preto, 720.116
  • Osasco, 701.428
  • Sorocaba, 695.328

Vamos localizá-los e destacar em nosso mapa. Para isto, criaremos um data frame com os nove municípios. No arquivo do IBGE…

municipios <- sf::st_read("IBGE/SP_Municipios_2021.shp", quiet=TRUE)
print(names(municipios))
[1] "CD_MUN"   "NM_MUN"   "SIGLA"    "AREA_KM2" "geometry"

…os nomes dos municípios estão na coluna NM_MUN. É preciso verificar se estão grafados da mesma forma que estão na Wikipedia (neste caso não houve problema). O R Script que implementa esta abordagem e destaca os municípios em amarelo (RGB “#F7CB45”) é o seguinte:

Brazil <- rnaturalearth::ne_states(country="Brazil", 
                                   returnclass="sf")
SP <- Brazil[Brazil$name_pt=="São Paulo",]

municipios <- sf::st_read(file.path("IBGE",
                                    "SP_Municipios_2021.shp"), 
                          quiet=TRUE)

# data frame com os municipios mais populosos
# fonte: https://pt.wikipedia.org/
mun <- c("São Paulo", 
         "Guarulhos",
         "Campinas",
         "São Bernardo do Campo",
         "São José dos Campos",
         "Santo André",
         "Ribeirão Preto",
         "Osasco",
         "Sorocaba")
pop <- c(12396372,
          1404694,
          1223237,
           849874,
           737310,
           723889,
           720116,
           701428,
           695328)
dtmun <- data.frame(mun,pop)
names(dtmun) <- c("municipio","populacao")

grf <- ggplot2::ggplot()+
  ggplot2::geom_sf(data=SP, fill="#CAE0AB", color="black")+
  ggplot2::geom_sf(data=municipios, fill="transparent", color="#507052", size=0.2)

# adiciona os municipos selecionados ao grafico, em amarelo
for (m in dtmun$municipio)
{
  shapemun <- municipios[municipios$NM_MUN==m,]
  if(nrow(shapemun)>0)
  {
    grf <- grf +
      ggplot2::geom_sf(data=shapemun, fill="#F7CB45", color="#507052", size=0.2)
  } else
  {
    cat("\nAlerta: municipio ",m," não foi localizado no IBGE.\n")
  }
  
}

print(grf)
Implementado com demo_ggplot_SaoPaulo03.R.

Longitude e Latitude

Uma outra alternativa é utilizar longitude e latitude das cidades, em vez de desenhar o contorno dos municípios. Para isto é necessário achar estes dados em algum outro lugar.

O IBGE as tem neste link. Há algumas pastas, incluindo shapefiles, que não explorei porque buscava uma planilha. Não está facilitado: os dados estão no formato .mdb (do Microsoft Access), .kml (do Google Earth) e .dbf (database SQL). Nada em universais arquivos-texto ou pelo menos em um .xls ou .xlsx do popular Excel.

Destes, o mais simples para o R é .dbf, que pode ser lido com

path = file.path("IBGE","BR_Localidades_2010_v1.dbf")
db <- foreign::read.dbf(path)
print(names(db))
 [1] "ID"         "CD_GEOCODI" "TIPO"       "CD_GEOCODB" "NM_BAIRRO" 
 [6] "CD_GEOCODS" "NM_SUBDIST" "CD_GEOCODD" "NM_DISTRIT" "CD_GEOCODM"
[11] "NM_MUNICIP" "NM_MICRO"   "NM_MESO"    "NM_UF"      "CD_NIVEL"  
[16] "CD_CATEGOR" "NM_CATEGOR" "NM_LOCALID" "LONG"       "LAT"       
[21] "ALT"        "GMRotation"
Implementado com demo_read_dbf.R.

Este arquivo, porém, é trabalhoso para arrumar. A coluna que nos interessa é NM_MUNICIP mas seu conteúdo (talvez por usarem sistemas Windows) tem problemas com acentuação (e.g., MISS\(\text{\xc3O}\)VELHA, HIDROL\(\text{\xc2}\)NDIA e GUAI\(\text{\xda}\)BA). Teríamos que encontrar a codificação de todos os possíveis acentos e fazer substituições antes de podermos usar o arquivo.

Neste tipo de situação, é melhor tentar algo mais fácil (deixo esta narreativa para que vejam que nem sempre os dados estão em formato conveniente; sugerimos que somente invista tempo em arrumar este tipo de dado se não houver alternativa).

Encontrei latitudes e longitudes de todos os munícipiso brasileiros (5565 no total) em formato em Excel na página

https://groups.google.com/g/qgisbrasil/c/0kJDo2KH0is

José Herlânio de Lima solicitou planilha dos municípios brasileiros ao IBGE, que lhe enviou esta planilha, a qual ele gentilmente compartilhou. Baixamos esta planilha na mesma pasta onde estão os shapefiles. A coluna de interesse é NOME_MUNICIPIO e seu único inconveniente é que os nomes estão grafados inteiramente em maiúsculas (o que se resolve com a funcão toupper — mas não sei se Windows consegue fazer isto corretamente). O procedimento é:

Warning in readLines("demo_addlonlat.R"): linha final incompleta encontrada em
'demo_addlonlat.R'
path = file.path("IBGE",
                 "anexo_16261_Coordenadas_Sedes_5565_Municípios_2010.xls")
dtibge <- readxl::read_excel(path)
print(names(dtibge))

# data frame com os municipios mais populosos
# fonte: https://pt.wikipedia.org/
mun <- c("São Paulo", 
         "Guarulhos",
         "Campinas",
         "São Bernardo do Campo",
         "São José dos Campos",
         "Santo André",
         "Ribeirão Preto",
         "Osasco",
         "Sorocaba")
pop <- c(12396372,
         1404694,
         1223237,
         849874,
         737310,
         723889,
         720116,
         701428,
         695328)
dtmun <- data.frame(mun,pop)
names(dtmun) <- c("municipio","populacao")

# adicionando latitude e longitude
dtmun$lat <- NA
dtmun$lon <- NA
for (m in 1:nrow(dtmun))
{
  row <- which(dtibge$NOME_MUNICIPIO==toupper(dtmun$municipio[m]))
  if (length(row)>0)
  {
    if(length(row)>1)
    {
      cat("\nAlerta: ",dtmun$municipio[m],"aparece em duas localizacoes: ",row,"\n")
    }
    dtmun$lat[m] <- dtibge$LATITUDE[row]
    dtmun$lon[m] <- dtibge$LONGITUDE[row]
  }  else
  {
    cat("\nAlerta: municipio ",m," não foi localizado na planilha do IBGE.\n")
  }
}
prmatrix(dtmun,quote=FALSE,rowlab=rep("",nrow(dtmun)))
[1] "GEOCODIGO_MUNICIPIO" "NOME_MUNICIPIO"      "LONGITUDE"          
[4] "LATITUDE"           

Alerta:  Santo André aparece em duas localizacoes:  1418 3801 
Warning in dtmun$lat[m] <- dtibge$LATITUDE[row]: número de itens a substituir
não é um múltiplo do comprimento do substituto
Warning in dtmun$lon[m] <- dtibge$LONGITUDE[row]: número de itens a substituir
não é um múltiplo do comprimento do substituto
 municipio             populacao lat        lon      
 São Paulo             12396372  -23.567387 -46.57038
 Guarulhos              1404694  -23.468506 -46.53108
 Campinas               1223237  -22.907343 -47.06016
 São Bernardo do Campo   849874  -23.710305 -46.55026
 São José dos Campos     737310  -23.184062 -45.88418
 Santo André             723889   -7.220489 -36.63098
 Ribeirão Preto          720116  -21.184835 -47.80548
 Osasco                  701428  -23.533612 -46.78881
 Sorocaba                695328  -23.499323 -47.45785
Implementado com demo_addlonlat.R.

Esta planilha tem um problema grave: existem os munícipios, mas não consta o Estado. Existem duas cidades chamadas “Santo André”. Podemos ver que a primeira Santo André encontrada foi tranferida para dtmun e esta não é a do Estado de São Paulo.

Para saber qual é a de São Paulo, podemos verificar se suas coordenadas estão entre as latitudes e longitudes máximas do shapefile do Estado de São Paulo (armazenado em SP.

Uma forma de resolver este problema é utilizar a função sf::st_bbox, aceitando apenas as localidades cujas coordenadas estejam dentro dos limites, como fizemos nesta segunda versão:

Warning in readLines("demo_addlonlat02.R"): linha final incompleta encontrada
em 'demo_addlonlat02.R'
# shapefile de São Paulo
Brazil <- rnaturalearth::ne_states(country="Brazil", returnclass="sf")
SP <- Brazil[Brazil$name_pt=="São Paulo",]

# data frame com os municipios mais populosos
# fonte: https://pt.wikipedia.org/
mun <- c("São Paulo", 
         "Guarulhos",
         "Campinas",
         "São Bernardo do Campo",
         "São José dos Campos",
         "Santo André",
         "Ribeirão Preto",
         "Osasco",
         "Sorocaba")
pop <- c(12396372,
         1404694,
         1223237,
         849874,
         737310,
         723889,
         720116,
         701428,
         695328)
dtmun <- data.frame(mun,pop)
names(dtmun) <- c("municipio","populacao")

# adicionando latitude e longitude
# coordenadas do IBGE
path = file.path("IBGE",
                 "anexo_16261_Coordenadas_Sedes_5565_Municípios_2010.xls")
dtibge <- readxl::read_excel(path)
print(names(dtibge))
# pega o boundary box 
limites <- sf::st_bbox(SP)
# cria as colunas lat e lon
dtmun$lat <- NA
dtmun$lon <- NA
warnings <- ""
for (m in 1:nrow(dtmun))
{
  if(!is.na(dtmun$lat[m])) {next}
  row <- which(dtibge$NOME_MUNICIPIO==toupper(dtmun$municipio[m]))
  if (length(row)>0)
  {
    if(length(row)>1)
    {
      warnings <- paste0(warnings,"\nAlerta: ",dtmun$municipio[m]," aparece em duas localizacoes: ")
      for (r in 1:length(row))
      {
        warnings <- paste0(warnings,row[r]," ")
      }
      prmatrix(dtibge[row,])
      row.org <- row
      for (r in 1:length(row.org))
      {
        # latitude fora da faixa
        if (dtibge$LONGITUDE[row.org[r]] < as.numeric(limites[1]) |
            dtibge$LONGITUDE[row.org[r]] > as.numeric(limites[3]))
        {
          row[r] <- NA
        }
        # longitude fora da faixa
        if (dtibge$LATITUDE[row.org[r]] < as.numeric(limites[2]) |
            dtibge$LATITUDE[row.org[r]] > as.numeric(limites[4]))
        {
          row[r] <- NA
        }
      }
      row <- row[!is.na(row)]
      if (length(row)==1)
      {
        warnings <- paste0(warnings,"\n\tcorrigido, escolhida a linha ",row," \n")
      } else
      {
        warnings <- paste0(warnings,"\n\tnao corrigido, ha duas localidades com o mesmo nome em São Paulo: ",row,"\n")
      }
    }
    dtmun$lat[m] <- dtibge$LATITUDE[row]
    dtmun$lon[m] <- dtibge$LONGITUDE[row]
  }  else
  {
    warnings <- paste0(warnings,"\nAlerta: municipio ",m," não foi localizado na planilha do IBGE.\n")
  }
}
cat(warnings)
prmatrix(dtmun,quote=FALSE,rowlab=rep("",nrow(dtmun)))
[1] "GEOCODIGO_MUNICIPIO" "NOME_MUNICIPIO"      "LONGITUDE"          
[4] "LATITUDE"           
     GEOCODIGO_MUNICIPIO NOME_MUNICIPIO LONGITUDE   LATITUDE    
[1,] "2513851"           "SANTO ANDRÉ"  "-36.63098" " -7.220489"
[2,] "3547809"           "SANTO ANDRÉ"  "-46.53087" "-23.657510"

Alerta: Santo André aparece em duas localizacoes: 1418 3801 
    corrigido, escolhida a linha 3801 
 municipio             populacao lat       lon      
 São Paulo             12396372  -23.56739 -46.57038
 Guarulhos              1404694  -23.46851 -46.53108
 Campinas               1223237  -22.90734 -47.06016
 São Bernardo do Campo   849874  -23.71030 -46.55026
 São José dos Campos     737310  -23.18406 -45.88418
 Santo André             723889  -23.65751 -46.53087
 Ribeirão Preto          720116  -21.18483 -47.80548
 Osasco                  701428  -23.53361 -46.78881
 Sorocaba                695328  -23.49932 -47.45785
Implementado com demo_addlonlat02.R.

Agora obtivemos as coordenadas corretas, de forma que podemos colocá-las no mapa (vamos, também, colocar as cidades em ordem decrescente de população, transformando o município em fator, com a ordem adequada, para que o ggplot2::ggplot a utilize), implementando:

Brazil <- rnaturalearth::ne_states(country="Brazil", 
                                   returnclass="sf")
SP <- Brazil[Brazil$name_pt=="São Paulo",]

# data frame com os municipios mais populosos
# fonte: https://pt.wikipedia.org/
mun <- c("São Paulo", 
         "Guarulhos",
         "Campinas",
         "São Bernardo do Campo",
         "São José dos Campos",
         "Santo André",
         "Ribeirão Preto",
         "Osasco",
         "Sorocaba")
pop <- c(12396372,
          1404694,
          1223237,
           849874,
           737310,
           723889,
           720116,
           701428,
           695328)
dtmun <- data.frame(mun,pop)
names(dtmun) <- c("municipio","populacao")
#m ordem decrescente do tamanho da populacao
dtmun <- dtmun[order(dtmun$populacao, decreasing=TRUE),]
# forcando a mesma ordem para municipio ser fator
# (serve para o ggplot utilizar esta ordem)
dtmun$municipio <- factor(dtmun$municipio,
                          levels=dtmun$municipio)

# adicionando latitude e longitude
# coordenadas do IBGE
path = file.path("IBGE",
                 "anexo_16261_Coordenadas_Sedes_5565_Municípios_2010.xls")
dtibge <- readxl::read_excel(path)
print(names(dtibge))
# pega o boundary box 
limites <- sf::st_bbox(SP)
# cria as colunas lat e lon
dtmun$lat <- NA
dtmun$lon <- NA
for (m in 1:nrow(dtmun))
{
  row <- which(dtibge$NOME_MUNICIPIO==toupper(dtmun$municipio[m]))
  if (length(row)>0)
  {
    if(length(row)>1)
    {
      cat("\nAlerta: ",as.character(dtmun$municipio[m]),"aparece em duas localizacoes: ",row,"\n")
      prmatrix(dtibge[row,])
      row.org <- row
      for (r in 1:length(row.org))
      {
        # latitude fora da faixa
        if (dtibge$LONGITUDE[row.org[r]] < as.numeric(limites[1]) |
            dtibge$LONGITUDE[row.org[r]] > as.numeric(limites[3]))
        {
          row[r] <- NA
        }
        # longitude fora da faixa
        if (dtibge$LATITUDE[row.org[r]] < as.numeric(limites[2]) |
            dtibge$LATITUDE[row.org[r]] > as.numeric(limites[4]))
        {
          row[r] <- NA
        }
      }
      row <- row[!is.na(row)]
      if (length(row)==1)
      {
        cat("\tcorrigido, escolhida a linha ",row,"\n\n")
      } else
      {
        cat("\tnao corrigido, ha duas localidades com o mesmo nome em São Paulo: ",row,"\n\n")
      }
    }
    dtmun$lat[m] <- dtibge$LATITUDE[row]
    dtmun$lon[m] <- dtibge$LONGITUDE[row]
  }  else
  {
    cat("\nAlerta: municipio ",m," não foi localizado na planilha do IBGE.\n")
  }
}

grf <- ggplot2::ggplot()+
  ggplot2::geom_sf(data=SP, fill="#CAE0AB", color="black")+
  ggplot2::geom_point(data=dtmun, ggplot2::aes(x=lon,y=lat,color=municipio))+
  ggplot2::labs(x="Longitude",y="Latitude",color="Cidade")
print(grf)
[1] "GEOCODIGO_MUNICIPIO" "NOME_MUNICIPIO"      "LONGITUDE"          
[4] "LATITUDE"           

Alerta:  Santo André aparece em duas localizacoes:  1418 3801 
     GEOCODIGO_MUNICIPIO NOME_MUNICIPIO LONGITUDE   LATITUDE    
[1,] "2513851"           "SANTO ANDRÉ"  "-36.63098" " -7.220489"
[2,] "3547809"           "SANTO ANDRÉ"  "-46.53087" "-23.657510"
    corrigido, escolhida a linha  3801
Implementado com demo_ggplot_SaoPaulo04.R.

Neste caso as coordenadas foram adicionadas com ggplot2::geom_point que utiliza o dataframe que criamos, dtmun, informando quais colunas coordenadas para longitude e latitude, além mudar a cor de acordo com a coluna dtmun$municipio. Usamos, também, a função ggplot2::labs para nomear os eixos e a legenda.

São Paulo em contexto

Podemos localizar o Estado de São Paulo dentro do Brasil, e o Brasil dentro do continente com [demo_ggplot_Brazil.R]:

AmericaSul <- rnaturalearth::ne_countries(continent = "South America", returnclass="sf")
Brazil <- rnaturalearth::ne_states(country="Brazil", returnclass="sf")
SP <- Brazil[Brazil$name_pt=="São Paulo",]

grf <- ggplot2::ggplot()+
  ggplot2::geom_sf(data=AmericaSul, fill="white", color="gray")+
  ggplot2::geom_sf(data=Brazil, fill="lightgray", color="gray")+
  ggplot2::geom_sf(data=SP, fill="#CAE0AB", color="black")+
  ggplot2::geom_point(data=dtmun, ggplot2::aes(x=lon,y=lat,color=municipio))+
  ggplot2::labs(x="Longitude",y="Latitude",color="Cidade")
print(grf)
Implementado com demo_ggplot_Brazil.R.

Este mapa tem o defeito de fazer a América do Sul, sem a conexão com a América Central, parecer-se com uma ilha. Há um problema como os dados estão organizados nos shapefiles de rnaturalearth. Observe:

rnaturalearth::ne_countries(returnclass="sf")
dtpaises <- rnaturalearth::ne_countries(returnclass="sf")
cat("\nCampos:\n")
print(names(dtpaises))
cat("\nContinentes:\n")
print(unique(dtpaises$continent))

Campos:
  [1] "featurecla" "scalerank"  "labelrank"  "sovereignt" "sov_a3"    
  [6] "adm0_dif"   "level"      "type"       "tlc"        "admin"     
 [11] "adm0_a3"    "geou_dif"   "geounit"    "gu_a3"      "su_dif"    
 [16] "subunit"    "su_a3"      "brk_diff"   "name"       "name_long" 
 [21] "brk_a3"     "brk_name"   "brk_group"  "abbrev"     "postal"    
 [26] "formal_en"  "formal_fr"  "name_ciawf" "note_adm0"  "note_brk"  
 [31] "name_sort"  "name_alt"   "mapcolor7"  "mapcolor8"  "mapcolor9" 
 [36] "mapcolor13" "pop_est"    "pop_rank"   "pop_year"   "gdp_md"    
 [41] "gdp_year"   "economy"    "income_grp" "fips_10"    "iso_a2"    
 [46] "iso_a2_eh"  "iso_a3"     "iso_a3_eh"  "iso_n3"     "iso_n3_eh" 
 [51] "un_a3"      "wb_a2"      "wb_a3"      "woe_id"     "woe_id_eh" 
 [56] "woe_note"   "adm0_iso"   "adm0_diff"  "adm0_tlc"   "adm0_a3_us"
 [61] "adm0_a3_fr" "adm0_a3_ru" "adm0_a3_es" "adm0_a3_cn" "adm0_a3_tw"
 [66] "adm0_a3_in" "adm0_a3_np" "adm0_a3_pk" "adm0_a3_de" "adm0_a3_gb"
 [71] "adm0_a3_br" "adm0_a3_il" "adm0_a3_ps" "adm0_a3_sa" "adm0_a3_eg"
 [76] "adm0_a3_ma" "adm0_a3_pt" "adm0_a3_ar" "adm0_a3_jp" "adm0_a3_ko"
 [81] "adm0_a3_vn" "adm0_a3_tr" "adm0_a3_id" "adm0_a3_pl" "adm0_a3_gr"
 [86] "adm0_a3_it" "adm0_a3_nl" "adm0_a3_se" "adm0_a3_bd" "adm0_a3_ua"
 [91] "adm0_a3_un" "adm0_a3_wb" "continent"  "region_un"  "subregion" 
 [96] "region_wb"  "name_len"   "long_len"   "abbrev_len" "tiny"      
[101] "homepart"   "min_zoom"   "min_label"  "max_label"  "label_x"   
[106] "label_y"    "ne_id"      "wikidataid" "name_ar"    "name_bn"   
[111] "name_de"    "name_en"    "name_es"    "name_fa"    "name_fr"   
[116] "name_el"    "name_he"    "name_hi"    "name_hu"    "name_id"   
[121] "name_it"    "name_ja"    "name_ko"    "name_nl"    "name_pl"   
[126] "name_pt"    "name_ru"    "name_sv"    "name_tr"    "name_uk"   
[131] "name_ur"    "name_vi"    "name_zh"    "name_zht"   "fclass_iso"
[136] "tlc_diff"   "fclass_tlc" "fclass_us"  "fclass_fr"  "fclass_ru" 
[141] "fclass_es"  "fclass_cn"  "fclass_tw"  "fclass_in"  "fclass_np" 
[146] "fclass_pk"  "fclass_de"  "fclass_gb"  "fclass_br"  "fclass_il" 
[151] "fclass_ps"  "fclass_sa"  "fclass_eg"  "fclass_ma"  "fclass_pt" 
[156] "fclass_ar"  "fclass_jp"  "fclass_ko"  "fclass_vn"  "fclass_tr" 
[161] "fclass_id"  "fclass_pl"  "fclass_gr"  "fclass_it"  "fclass_nl" 
[166] "fclass_se"  "fclass_bd"  "fclass_ua"  "geometry"  

Continentes:
[1] "Oceania"                 "Africa"                 
[3] "North America"           "Asia"                   
[5] "South America"           "Europe"                 
[7] "Seven seas (open ocean)" "Antarctica"
Implementado com demo_naotemCentral.R.

Não existe “Central America”. Se usarmos a América do Norte fica assim:

AmericaNorte <- rnaturalearth::ne_countries(continent = "North America", returnclass="sf")
AmericaSul <- rnaturalearth::ne_countries(continent = "South America", returnclass="sf")
Brazil <- rnaturalearth::ne_states(country="Brazil", returnclass="sf")
SP <- Brazil[Brazil$name_pt=="São Paulo",]

grf <- ggplot2::ggplot()+
  ggplot2::geom_sf(data=AmericaNorte, fill="white", color="gray")+
  ggplot2::geom_sf(data=AmericaSul, fill="white", color="gray")+
  ggplot2::geom_sf(data=Brazil, fill="lightgray", color="gray")+
  ggplot2::geom_sf(data=SP, fill="#CAE0AB", color="black")+
  ggplot2::geom_point(data=dtmun, ggplot2::aes(x=lon,y=lat,color=municipio))+
  ggplot2::labs(x="Longitude",y="Latitude",color="Cidade")
print(grf)
Implementado com demo_ggplot_Brazil02.R.

Depois de alguma investigação percebe-se que a América do Norte está subdivida em três áreas, sendo “Central America” uma delas, além da América do Norte propriamente dita e o Caribe:

rnaturalearth::ne_countries(returnclass="sf")
dtpaises <- rnaturalearth::ne_countries(returnclass="sf")
cat("\nCampos:\n")
print(names(dtpaises))
cat("\nContinentes:\n")
print(unique(dtpaises$continent))
cat("\nSub regioes:\n")
print(unique(dtpaises$subregion[dtpaises$continent=="North America"]))

Campos:
  [1] "featurecla" "scalerank"  "labelrank"  "sovereignt" "sov_a3"    
  [6] "adm0_dif"   "level"      "type"       "tlc"        "admin"     
 [11] "adm0_a3"    "geou_dif"   "geounit"    "gu_a3"      "su_dif"    
 [16] "subunit"    "su_a3"      "brk_diff"   "name"       "name_long" 
 [21] "brk_a3"     "brk_name"   "brk_group"  "abbrev"     "postal"    
 [26] "formal_en"  "formal_fr"  "name_ciawf" "note_adm0"  "note_brk"  
 [31] "name_sort"  "name_alt"   "mapcolor7"  "mapcolor8"  "mapcolor9" 
 [36] "mapcolor13" "pop_est"    "pop_rank"   "pop_year"   "gdp_md"    
 [41] "gdp_year"   "economy"    "income_grp" "fips_10"    "iso_a2"    
 [46] "iso_a2_eh"  "iso_a3"     "iso_a3_eh"  "iso_n3"     "iso_n3_eh" 
 [51] "un_a3"      "wb_a2"      "wb_a3"      "woe_id"     "woe_id_eh" 
 [56] "woe_note"   "adm0_iso"   "adm0_diff"  "adm0_tlc"   "adm0_a3_us"
 [61] "adm0_a3_fr" "adm0_a3_ru" "adm0_a3_es" "adm0_a3_cn" "adm0_a3_tw"
 [66] "adm0_a3_in" "adm0_a3_np" "adm0_a3_pk" "adm0_a3_de" "adm0_a3_gb"
 [71] "adm0_a3_br" "adm0_a3_il" "adm0_a3_ps" "adm0_a3_sa" "adm0_a3_eg"
 [76] "adm0_a3_ma" "adm0_a3_pt" "adm0_a3_ar" "adm0_a3_jp" "adm0_a3_ko"
 [81] "adm0_a3_vn" "adm0_a3_tr" "adm0_a3_id" "adm0_a3_pl" "adm0_a3_gr"
 [86] "adm0_a3_it" "adm0_a3_nl" "adm0_a3_se" "adm0_a3_bd" "adm0_a3_ua"
 [91] "adm0_a3_un" "adm0_a3_wb" "continent"  "region_un"  "subregion" 
 [96] "region_wb"  "name_len"   "long_len"   "abbrev_len" "tiny"      
[101] "homepart"   "min_zoom"   "min_label"  "max_label"  "label_x"   
[106] "label_y"    "ne_id"      "wikidataid" "name_ar"    "name_bn"   
[111] "name_de"    "name_en"    "name_es"    "name_fa"    "name_fr"   
[116] "name_el"    "name_he"    "name_hi"    "name_hu"    "name_id"   
[121] "name_it"    "name_ja"    "name_ko"    "name_nl"    "name_pl"   
[126] "name_pt"    "name_ru"    "name_sv"    "name_tr"    "name_uk"   
[131] "name_ur"    "name_vi"    "name_zh"    "name_zht"   "fclass_iso"
[136] "tlc_diff"   "fclass_tlc" "fclass_us"  "fclass_fr"  "fclass_ru" 
[141] "fclass_es"  "fclass_cn"  "fclass_tw"  "fclass_in"  "fclass_np" 
[146] "fclass_pk"  "fclass_de"  "fclass_gb"  "fclass_br"  "fclass_il" 
[151] "fclass_ps"  "fclass_sa"  "fclass_eg"  "fclass_ma"  "fclass_pt" 
[156] "fclass_ar"  "fclass_jp"  "fclass_ko"  "fclass_vn"  "fclass_tr" 
[161] "fclass_id"  "fclass_pl"  "fclass_gr"  "fclass_it"  "fclass_nl" 
[166] "fclass_se"  "fclass_bd"  "fclass_ua"  "geometry"  

Continentes:
[1] "Oceania"                 "Africa"                 
[3] "North America"           "Asia"                   
[5] "South America"           "Europe"                 
[7] "Seven seas (open ocean)" "Antarctica"             

Sub regioes:
[1] "Northern America" "Caribbean"        "Central America"
Implementado com demo_naotemCentral02.R.

Basta, então, filtrá-la para colocar a imagem:

AmericaNorte <- rnaturalearth::ne_countries(continent = "North America", returnclass="sf")
AmericaCentral <- AmericaNorte[AmericaNorte$subregion=="Central America",]
Caribe <- AmericaNorte[AmericaNorte$subregion=="Caribbean",]
AmericaSul <- rnaturalearth::ne_countries(continent = "South America", returnclass="sf")
Brazil <- rnaturalearth::ne_states(country="Brazil", returnclass="sf")
SP <- Brazil[Brazil$name_pt=="São Paulo",]

grf <- ggplot2::ggplot()+
  ggplot2::geom_sf(data=AmericaCentral, fill="white", color="gray")+
  ggplot2::geom_sf(data=Caribe, fill="white", color="gray")+
  ggplot2::geom_sf(data=AmericaSul, fill="white", color="gray")+
  ggplot2::geom_sf(data=Brazil, fill="lightgray", color="gray")+
  ggplot2::geom_sf(data=SP, fill="#CAE0AB", color="black")+
  ggplot2::geom_point(data=dtmun, ggplot2::aes(x=lon,y=lat,color=municipio))+
  ggplot2::labs(x="Longitude",y="Latitude",color="Cidade")
print(grf)
Implementado com demo_ggplot_Brazil03.R.

Outra possibilidade é cortar o gráfico, para aparecer somente o que nos interessa (o Estado de São Paulo e sua vizinhança):

Brazil <- rnaturalearth::ne_states(country="Brazil", returnclass="sf")
SP <- Brazil[Brazil$name_pt=="São Paulo",]

grf <- ggplot2::ggplot()+
  ggplot2::geom_sf(data=Brazil, fill="lightgray", color="gray")+
  ggplot2::geom_sf(data=SP, fill="#CAE0AB", color="black")+
  ggplot2::geom_point(data=dtmun, ggplot2::aes(x=lon,y=lat,color=municipio))+
  ggplot2::labs(x="Longitude",y="Latitude",color="Cidade")+
  ggplot2::ylim(-26,-19)+
  ggplot2::xlim(-55,-43)
print(grf)
Implementado com demo_ggplot_Brazil04.R.

Não satisfeito com o quadriculado em cinza aparecendo ao fundo? Há temas e configurações que podem ser usadas:

Brazil <- rnaturalearth::ne_states(country="Brazil", returnclass="sf")
SP <- Brazil[Brazil$name_pt=="São Paulo",]

grf <- ggplot2::ggplot()+
  ggplot2::geom_sf(data=Brazil, fill="lightgray", color="gray")+
  ggplot2::geom_sf(data=SP, fill="#CAE0AB", color="black")+
  ggplot2::geom_point(data=dtmun, ggplot2::aes(x=lon,y=lat,color=municipio))+
  ggplot2::labs(x="Longitude",y="Latitude",color="Cidade")+
  ggplot2::ylim(-26,-19)+
  ggplot2::xlim(-55,-43)+
  ggplot2::theme_light()+
  ggplot2::theme(panel.grid.major = ggplot2::element_line(color = gray(.5),
                                                          linetype = 'dashed',
                                                          size = 0.2),
        panel.grid.minor = ggplot2::element_blank(),
        panel.background = ggplot2::element_rect(fill=NA,color = 'black'),
        panel.ontop = TRUE)
print(grf)
Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
ℹ Please use the `linewidth` argument instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.
Implementado com demo_ggplot_Brazil05.R.

Falta o Oceano. Um pedaço que parece adequado está em MarineRegions.org. O shapefile foi colocado nem pasta chamada Ocean e tem o nome eez_iho.shp. Obtém-se:

Ocean <- sf::st_read(file.path("Ocean","eez_iho.shp"), quiet=TRUE)
Brazil <- rnaturalearth::ne_states(country="Brazil", returnclass="sf")
SP <- Brazil[Brazil$name_pt=="São Paulo",]

grf <- ggplot2::ggplot()+
  ggplot2::geom_sf(data=Ocean, fill="#A0D8D8", color="#A0D8D8")+
  ggplot2::geom_sf(data=Brazil, fill="lightgray", color="gray")+
  ggplot2::geom_sf(data=SP, fill="#CAE0AB", color="black")+
  ggplot2::geom_point(data=dtmun, ggplot2::aes(x=lon,y=lat,color=municipio))+
  ggplot2::labs(x="Longitude",y="Latitude",color="Cidade")+
  ggplot2::ylim(-26,-19)+
  ggplot2::xlim(-55,-43)+
  ggplot2::theme_light()+
  ggplot2::theme(panel.grid.major = ggplot2::element_line(color = gray(.5),
                                                          linetype = 'dashed',
                                                          size = 0.2),
        panel.grid.minor = ggplot2::element_blank(),
        panel.background = ggplot2::element_rect(fill=NA,color = 'black'),
        panel.ontop = TRUE)
print(grf)
Implementado com demo_ggplot_Brazil06.R.

O área dos círculos pode ser proporcional à população (o município de São Paulo concentra quase 64% dos habitantes destas nove regiões). Além disto, podemos controlar as cores das cidades selecionadas (usando ggplot2::scale_color_manual), usando uma paleta que facilita o discernimento para daltônicos. Por exemplo:

source("eiras.friendlycolor.R")

Ocean <- sf::st_read(file.path("Ocean","eez_iho.shp"), quiet=TRUE)
Brazil <- rnaturalearth::ne_states(country="Brazil", returnclass="sf")
SP <- Brazil[Brazil$name_pt=="São Paulo",]

# circle diameters proportional to population
dtmun$pop.perc <- (dtmun$populacao/sum(dtmun$populacao)*100)
# area do circulo é area=pi*(radius^2) -> radius <- (area/pi)^0.5
# diameter = 2 radius
# k se precisar para ajeitar a escala
k <- 0.7
dtmun$diameter <- 2*((dtmun$pop.perc/pi)^0.5)*k
prmatrix(dtmun,quote=FALSE,rowlab=rep("",nrow(dtmun)))

grf <- ggplot2::ggplot()+
  ggplot2::geom_sf(data=Ocean, fill="#A0D8D8", color="#A0D8D8")+
  ggplot2::geom_sf(data=Brazil, fill="lightgray", color="gray")+
  ggplot2::geom_sf(data=SP, fill="#CAE0AB", color="black")+
  ggplot2::geom_point(data=dtmun, size=dtmun$diameter[1:nrow(dtmun)], 
                      ggplot2::aes(x=lon,y=lat,color=municipio))+
  ggplot2::scale_color_manual(labels = dtmun$municipio, 
                              values = c(eiras.friendlycolor(seq(1,25,by=6)),
                                         eiras.friendlycolor(seq(4,28,by=6)))) +
  ggplot2::labs(x="Longitude",y="Latitude",color="Cidade")+
  ggplot2::ylim(-26,-19)+
  ggplot2::xlim(-55,-43)+
  ggplot2::theme_light()+
  ggplot2::theme(panel.grid.major = ggplot2::element_line(color = gray(.5),
                                                          linetype = 'dashed',
                                                          size = 0.2),
        panel.grid.minor = ggplot2::element_blank(),
        panel.background = ggplot2::element_rect(fill=NA,color = 'black'),
        panel.ontop = TRUE)
print(grf)
 municipio             populacao lat       lon       pop.perc  diameter
 São Paulo             12396372  -23.56739 -46.57038 63.727195 6.305441
 Guarulhos              1404694  -23.46851 -46.53108  7.221243 2.122556
 Campinas               1223237  -22.90734 -47.06016  6.288409 1.980722
 São Bernardo do Campo   849874  -23.71030 -46.55026  4.369027 1.650994
 São José dos Campos     737310  -23.18406 -45.88418  3.790359 1.537777
 Santo André             723889  -23.65751 -46.53087  3.721364 1.523717
 Ribeirão Preto          720116  -21.18483 -47.80548  3.701968 1.519740
 Osasco                  701428  -23.53361 -46.78881  3.605897 1.499891
 Sorocaba                695328  -23.49932 -47.45785  3.574538 1.493355
Implementado com demo_ggplot_Brazil07.R.

Há vários outros recursos, ainda. Imagine a representação de seus dados como quiser, e encontrará alguém que fez algo similar para adaptar ao seu caso.

Dados de COVID

Entre outras fontes, o Ministério da Saúde disponibiliza dados de casos e óbitos por COVID-19 no Brasil, estratificados por Unidade da Federação e por Município em https://covid.saude.gov.br/.

Já modificaram o formato do arquivo algumas vezes. Desde o ano passado é atualizado um arquivo compactado no formato .zip, acessível a partir de um botão com o texto Arquivo CSV no canto superior-direito da página.

Parte da tela de https://covid.saude.gov.br/ (acesso em 21 de junho de 2022).

O arquivo compactado .rar baixado em 21 de junho de 2022 continha cinco arquivos .csv:

  • HIST_PAINEL_COVIDBR_2020_Parte1_20jun2022.csv
  • HIST_PAINEL_COVIDBR_2020_Parte2_20jun2022.csv
  • HIST_PAINEL_COVIDBR_2021_Parte1_20jun2022.csv
  • HIST_PAINEL_COVIDBR_2021_Parte2_20jun2022.csv
  • HIST_PAINEL_COVIDBR_2022_Parte1_20jun2022.csv

O painel muda, atualizando a situação. Em 31 de maio de 2023 e agosto de 2024 encontramos o seguinte:

Parte da tela de https://covid.saude.gov.br/ (acessos em 31 de maio de 2023 e 29 de agosto de 2024).

O arquivo baixado tem formato .zip e traz subdivisões com o formato:

  • HIST_PAINEL_COVIDBR_20XX_Parte1_DDmesYYYY.csv
  • HIST_PAINEL_COVIDBR_20XX_Parte2_DDmesYYYY.csv

onde 20XX vai de 2020 a 2024, Partes 1 e 2 correspondem ao primeiro e segundo semestres, e DDmesYYYY é a data da totalização dos dados. Atualmente os dados são apresentados com a soma de ocorrências de cada semana (no momento da pandemia o dado era diário).

Não encontrei documentação alguma sobre a organização dos dados nestes arquivos (e já fizeram mudanças ao longo do tempo) mas, através de um exercício rápido de exploração, é possível entendê-lo. Todos os arquivos têm a mesma estrutura, de forma que podem ser lidos e integrados em um único data frame. Descompactando o arquivo .zip em uma pasta subordinada MS_COVID, a leitura dos dados pode ser feita com:

COVID <- NULL
folder.in <- "MS_COVID"
folder.out <- "COVID2024"
sufixo <- "23ago2024"
year <- 2020
ok <- TRUE

try(dir.create(folder.out),silent=TRUE)

cat("\nReading COVID data:\n")
while(ok)
{
  for(p.aux in 1:2)
  {
    filename <- file.path(folder.in,
                          paste0("HIST_PAINEL_COVIDBR_",year,
                                 "_Parte", p.aux,"_",sufixo,".csv"))
    if(file.exists(filename))
    {
      dt_tmp <- read.csv2(filename)
      if(is.null(COVID))
      {
        COVID <- dt_tmp
      } else
      {
        COVID <- rbind(COVID,dt_tmp)
      }
      cat("\n",filename)
    } else
    {
      ok <- FALSE
    }
  }
  year <- year+1
}
rm(dt_tmp)
cat("\nFinished\n")

COVID$datanum <- as.integer(as.Date(COVID$data,origin = "1899-12-30"))+25569
saveRDS(COVID,file.path(folder.out,"COVID.rds"))

Este arquivo tem …

# exibe o nome das colunas
folder.out <- "COVID2024"
COVID <- readRDS(file.path(folder.out,"COVID.rds"))
print(names(COVID))
 [1] "regiao"                 "estado"                 "municipio"             
 [4] "coduf"                  "codmun"                 "codRegiaoSaude"        
 [7] "nomeRegiaoSaude"        "data"                   "semanaEpi"             
[10] "populacaoTCU2019"       "casosAcumulado"         "casosNovos"            
[13] "obitosAcumulado"        "obitosNovos"            "Recuperadosnovos"      
[16] "emAcompanhamentoNovos"  "interior.metropolitana" "datanum"               
total_linhas <- nrow(COVID)
cat(total_linhas," linhas",sep="")
9053077 linhas

Este número de linhas é explicado porque nem todas as colunas estão preenchidas ao longo do data frame. Por exemplo, observe que região contém “Brasil” em …

cat(sum(COVID$regiao=="Brasil")," linhas",sep="")
1642 linhas

e

datas <- COVID$data[COVID$regiao=="Brasil"]
cat(paste0("Total de datas diferentes: ",length(unique(datas)),
           "\nde ",min(datas)," a ",max(datas)))
Total de datas diferentes: 1642
de 2020-02-25 a 2024-08-23

São os 1642 dias de registro da pandemia, iniciando-se em “2020-02-25” e terminando em “2024-08-23”.

O número acumulado de casos e de óbitos podem ser exibidos:

datas <- COVID$data[COVID$regiao=="Brasil"]
datas.num <- as.integer(as.Date(datas,origin = "1899-12-30"))+25569

# strings dos meses e numero correspondente ao dia 1 de cada mes
ano.inicial <- as.numeric(substr(datas[1],1,4))
mes.inicial <- as.numeric(substr(datas[1],6,7))
ano.final <-as.numeric(substr(datas[length(datas)],1,4))
mes.final <- as.numeric(substr(datas[length(datas)],6,7))
label.datas <- c()
label.num <- c()
y.aux <- ano.inicial
m.aux <- mes.inicial-1; if(m.aux<=0){m.aux<-1}
mes.final <- mes.final+1; if(mes.final>12){mes.final<-12}
while (!(y.aux == ano.final & m.aux == mes.final))
{
  m.aux <- m.aux+1
  if (m.aux>12)
  {
    m.aux <- 1
    y.aux <- y.aux+1
  }
  label.datas <- c(label.datas,paste0(m.aux,"/",y.aux))
  label.num <- c(label.num,  as.integer(as.Date(paste0(y.aux,"-",m.aux,"-01"),origin = "1899-12-30"))+25569)
}

# graficos
casos.acm <- COVID$casosAcumulado[COVID$regiao=="Brasil"]
obitos.acm <- COVID$obitosAcumulado[COVID$regiao=="Brasil"]
plot(datas.num, casos.acm, 
     ylim = c(0,max(casos.acm,na.rm=TRUE)),
     xlab="", ylab="Número de casos (acumulado)", 
     type="l", lwd=2, col="darkblue", axes=FALSE)
axis(1, labels = label.datas, at = label.num, las=2)
axis(2)

plot(datas.num, obitos.acm, 
     ylim = c(0,max(obitos.acm,na.rm=TRUE)),
     xlab="", ylab="Número de óbitos (acumulado)", 
     type="l", lwd=2, col="darkred", axes=FALSE)
axis(1, labels = label.datas, at = label.num, las=2)
axis(2)

A partir dos dados conta-se 38867008 casos e 712957 mortes no Brasil, os mesmos valores que o painel registra.

Há mais conteúdo nestes dados, pois vimos que há 9053077 linhas no arquivo. Após um processo de tentativa-e-erro vasculhando o banco encontramos…

regioes <- unique(COVID$regiao)
cat("Regiões encontradas:",regioes,"\n")
Regiões encontradas: Brasil Norte Nordeste Sudeste Sul Centro-Oeste 
estados <- unique(COVID$estado)
estados <- estados[!is.na(estados)]
estados <- estados[nchar(estados)>0]
cat("Estados encontrados: ",estados," (total de ",length(estados)," estados)\n")
Estados encontrados:  RO AC AM RR PA AP TO MA PI CE RN PB PE AL SE BA MG ES RJ SP PR SC RS MS MT GO DF  (total de  27  estados)
municipios <- c()
for (e in estados)
{
  m <- unique(COVID$municipio[COVID$estado==e])
  m <- m[!is.na(m)]
  m <- m[nchar(m)>0]
  cat("\t",e,": ",length(m)," municipios\n")
  municipios <- c(municipios,m)
}
cat("\nTotal de municípios encontrados:",length(municipios),"\n")
     RO :  52  municipios
     AC :  22  municipios
     AM :  62  municipios
     RR :  15  municipios
     PA :  144  municipios
     AP :  16  municipios
     TO :  139  municipios
     MA :  217  municipios
     PI :  224  municipios
     CE :  184  municipios
     RN :  167  municipios
     PB :  223  municipios
     PE :  185  municipios
     AL :  102  municipios
     SE :  75  municipios
     BA :  417  municipios
     MG :  853  municipios
     ES :  78  municipios
     RJ :  92  municipios
     SP :  645  municipios
     PR :  399  municipios
     SC :  295  municipios
     RS :  497  municipios
     MS :  79  municipios
     MT :  141  municipios
     GO :  246  municipios
     DF :  1  municipios

Total de municípios encontrados: 5570

Este banco está organizado com base em regiao, estado e municipio. O segmento de interesse pode ser separado combinando-se esses três elementos:

  • quando só regiao está preenchida com “Brasil” (estado e municipio em branco), os dados são a totalização do país.
  • quando só regiao está preenchida (estado e municipio em branco), os dados são a totalização da região.
  • quando regiao e estado estão preenchidos (municipio em branco), os dados são a totalização do Estado.
  • quando estado e municipio estão preenchidos, os dados são a totalização do município.

Assim:

total.linhas <- 0

df_brasil <- COVID[COVID$regiao=="Brasil",]
cat("df_brasil: data frame com ",nrow(df_brasil)," linhas\n", sep="")
total.linhas <- total.linhas+nrow(df_brasil)

df_estados <- COVID[COVID$regiao!="Brasil" & nchar(COVID$estado)>0 & nchar(COVID$municipio)==0,]
cat("df_estados: data frame com ",nrow(df_estados)," linhas\n", sep="")
total.linhas <- total.linhas+nrow(df_estados)

estados <- unique(df_estados$estado)
for (idx in 1:length(estados))
{
  regiao <- unique(COVID$regiao[COVID$estado==estados[idx]])
  df_municipios <- COVID[COVID$regiao==regiao & COVID$estado==estados[idx] & nchar(COVID$municipio)>0,]
  cat(estados[idx], " (regiao ",regiao,"): data frame com ",nrow(df_municipios)," linhas, incluindo municipios\n", sep="")
  
  total.linhas <- total.linhas+nrow(df_municipios)
}

cat("\nTotal de linhas aproveitadas do arquivo: ",total.linhas," de ",nrow(COVID),"\n",sep="")

if (total.linhas==nrow(COVID)) {cat("\nTodas as linhas foram encontradas\n")}
df_brasil: data frame com 1642 linhas
df_estados: data frame com 78165 linhas
RO (regiao Norte): data frame com 83772 linhas, incluindo municipios
AC (regiao Norte): data frame com 35442 linhas, incluindo municipios
AM (regiao Norte): data frame com 99882 linhas, incluindo municipios
RR (regiao Norte): data frame com 24165 linhas, incluindo municipios
PA (regiao Norte): data frame com 231984 linhas, incluindo municipios
AP (regiao Norte): data frame com 25776 linhas, incluindo municipios
TO (regiao Norte): data frame com 223929 linhas, incluindo municipios
MA (regiao Nordeste): data frame com 349587 linhas, incluindo municipios
PI (regiao Nordeste): data frame com 360864 linhas, incluindo municipios
CE (regiao Nordeste): data frame com 296424 linhas, incluindo municipios
RN (regiao Nordeste): data frame com 269037 linhas, incluindo municipios
PB (regiao Nordeste): data frame com 359253 linhas, incluindo municipios
PE (regiao Nordeste): data frame com 298035 linhas, incluindo municipios
AL (regiao Nordeste): data frame com 164322 linhas, incluindo municipios
SE (regiao Nordeste): data frame com 120825 linhas, incluindo municipios
BA (regiao Nordeste): data frame com 671787 linhas, incluindo municipios
MG (regiao Sudeste): data frame com 1374183 linhas, incluindo municipios
ES (regiao Sudeste): data frame com 125658 linhas, incluindo municipios
RJ (regiao Sudeste): data frame com 148212 linhas, incluindo municipios
SP (regiao Sudeste): data frame com 1039095 linhas, incluindo municipios
PR (regiao Sul): data frame com 642789 linhas, incluindo municipios
SC (regiao Sul): data frame com 475245 linhas, incluindo municipios
RS (regiao Sul): data frame com 800667 linhas, incluindo municipios
MS (regiao Centro-Oeste): data frame com 127269 linhas, incluindo municipios
MT (regiao Centro-Oeste): data frame com 227151 linhas, incluindo municipios
GO (regiao Centro-Oeste): data frame com 396306 linhas, incluindo municipios
DF (regiao Centro-Oeste): data frame com 1611 linhas, incluindo municipios

Total de linhas aproveitadas do arquivo: 9053077 de 9053077

Todas as linhas foram encontradas

Desta forma pudemos produzir os gráficos da evolução diária de casos e óbitos e armazenamos os resultados em um PDF, indo até junho de 2023 (os dados deixaram de ser anotados diariamente para serem atualizados semanamente). Nesses gráficos, as linhas mais finas correspondem aos dados e a linha mais grossa é uma curva ajustada:

Criamos uma sequência de imagens para visualizar no mapa a progressão da pandemia no Brasil, colorindo os polígonos que representam as Unidades da Federação. O resultado está em COVID_UF.mp4 (esta é a versão de junho/2022).

Referências