Elaborando Mapas no R
Fellipe Gomes
24 de novembro de 2017
Introdução
Neste documento será visto diversas maneiras de se criar mapas de padrões de pontos e de dados de áres utilizando arquivos shapefile ou latitudes e longitudes com o R.
Leitura dos dados
Leitura do shape do mapa do brasil separado por municipios
#Leitura do shapefile
#mapa=readShapePoly("BRMUE250GC_SIR.shp") #Outra opcao de funcao para a leitura dos dados
mapa = readOGR("BRMUE250GC_SIR.shp")## OGR data source with driver: ESRI Shapefile
## Source: "BRMUE250GC_SIR.shp", layer: "BRMUE250GC_SIR"
## with 5572 features
## It has 2 fields#Alterando a codificação da base do arquivo .shp (shape):
mapa@data[,1]=str_conv(mapa@data[,1],encoding = "UTF8")
#Removendo acentos:
mapa@data[,1]=rm_accent(mapa@data[,1]) #Obs: esta funcao se encontra no arquivo Funcao-rm_accent no mesmo diretorio
#Removendo os tracos "-" e "'" :
mapa@data[,1]=str_replace_all(mapa@data[,1], "[']", " ")
mapa@data[,1]=str_replace_all(mapa@data[,1], "[-]", " ")
#Removendo espaços em branco desnecssarios:
mapa@data[,1]=str_trim(mapa@data[,1])
#Alterando o cabecalho:
names(mapa@data)[1]="NOME"
#A base organizada:
head(mapa@data)## NOME CD_GEOCMU
## 0 ALTA FLORESTA D OESTE 1100015
## 1 ARIQUEMES 1100023
## 2 CABIXI 1100031
## 3 CACOAL 1100049
## 4 CEREJEIRAS 1100056
## 5 COLORADO DO OESTE 1100064Leitura da base de dados que será utilizada
#Fazendo a leitura da base de dados que contem os dados:
library(readxl) #Pacote exigido para leitura de arquivo excell
dadjov <- as.data.frame(read_excel("Prrenchimento_sexo_idade_ign_por_mun.xlsx"))
#Removendo os tracos "-" e "'" :
dadjov$NOME=str_replace_all(dadjov$NOME, "[-]", " ")
head(dadjov)## MUNICIPIO_PFA FREQ BASE PFA PERC_TOTAL_REGISTROS PERC_IGNORADO FREQ_IGN
## 1 8314 3292 4.727502e-05 0.005413599 5
## 2 11001 1964 2.820551e-05 0.0021654396 2
## 3 8001 3042 4.368817e-05 0.0021654396 2
## 4 11002 8148 1.170239e-04 0.0173235167 16
## 5 12001 10273 1.475611e-04 0.0303161542 28
## 6 5001 891 1.280019e-05 0.005413599 5
## Percentual Ign Mun % Ign Mun Ranking A Perc ign Ranking Freq Ign
## 1 0.0015188335 0.15188335 3229 1325
## 2 0.0010183299 0.10183299 3800 2101
## 3 0.0006574622 0.06574622 4316 2135
## 4 0.0019636721 0.19636721 2831 1
## 5 0.0027255914 0.27255914 2313 243
## 6 0.0056116723 0.56116723 1277 1253
## NOME CODUF
## 1 ABADIA DE GOIAS 8
## 2 ABADIA DOS DOURADOS 11
## 3 ABADIANIA 8
## 4 ABAETE 11
## 5 ABAETETUBA 12
## 6 ABAIARA 5Unindo os dois arquivos
#Organizando os dados em ordem alfabetica de acordo com o nome dos municipios:
mapa@data=mapa@data%>%
arrange(NOME)
#Acrescentando a proporcao calculada dentro do shapefile ### IMPORTANTE: A VARIAVEL PRECISA ESTAR NA ORDEM DAS REGIOES
mapa@data[,3]=dadjov%>%
arrange(NOME)%>%
select(PERC_TOTAL_REGISTROS)
mapa@data[,4]=dadjov%>%
arrange(NOME)%>%
select(PERC_IGNORADO)
mapa@data[,5]=dadjov%>%
arrange(NOME)%>%
select(`FREQ BASE PFA`)
mapa@data[,6]=dadjov%>%
arrange(NOME)%>%
select(FREQ_IGN)
#Retirando NAs
for(i in 1:length(mapa@data[,3])){if(is.na(mapa@data[i,3])){mapa@data[i,3]=0}}
for(i in 1:length(mapa@data[,4])){if(is.na(mapa@data[i,4])){mapa@data[i,4]=0}}
for(i in 1:length(mapa@data[,5])){if(is.na(mapa@data[i,5])){mapa@data[i,5]=0}}
for(i in 1:length(mapa@data[,6])){if(is.na(mapa@data[i,6])){mapa@data[i,6]=0}}
#Por alguma razao esta variavel foi interpretada como caracter, entao, convertendo para numerico:
mapa@data$PERC_IGNORADO=as.numeric(mapa@data$PERC_IGNORADO)
#Retirando percentuais superiores a 1 (corrigir, alguns dados talvez precisem ser novamente gerados)
for(i in 1:length(mapa@data$PERC_IGNORADO)){
if(mapa@data$PERC_IGNORADO[i]>1){
mapa@data$PERC_IGNORADO[i]=0
}
}
summary(mapa@data)## NOME CD_GEOCMU PERC_TOTAL_REGISTROS PERC_IGNORADO
## Length:5572 1100015: 1 Min. :0.0000008 Min. :0.000000
## Class :character 1100023: 1 1st Qu.:0.0000186 1st Qu.:0.002165
## Mode :character 1100031: 1 Median :0.0000528 Median :0.004495
## 1100049: 1 Mean :0.0013114 Mean :0.013371
## 1100056: 1 3rd Qu.:0.0002006 3rd Qu.:0.010827
## 1100064: 1 Max. :1.0000000 Max. :0.989606
## (Other):5566
## FREQ BASE PFA FREQ_IGN
## Min. : 1 Min. : 0.00
## 1st Qu.: 669 1st Qu.: 2.00
## Median : 1581 Median : 4.00
## Mean : 12497 Mean : 16.58
## 3rd Qu.: 4948 3rd Qu.: 10.00
## Max. :7327342 Max. :5936.00
## Criando mapas utilizando o shapefile + banco de dados
1 método: plot()
A primeira maneira consiste em utilizar o comando plot() (comando nativo do R), veja:
#Numero de categorias
n.cat = 8
#Cores utilizadas para cada categoria
cores = c("#03E06EFF", "#38FA2AFF", "#9FFF15FF", "#DEFE0BFF", "#F3FB06FF", "#FFCD00FF", "#FF8500FF", "#FF3300FF")
#Criando os intervalos de classes para atribuir as cores escolhidas
intervalos=quantile(mapa@data$PERC_IGNORADO, probs = seq(0,1,0.125))
#Criando as classes:
classes <- classIntervals(mapa@data$PERC_IGNORADO, n.cat, style= "fixed",
fixedBreaks = round(intervalos,2))
#Associando uma cor a uma categoria
corcat <- findColours(classes, cores)
#Plot do mapa com as proporcoes
plot(mapa, col=corcat)
# #Acrescentando uma caixa a figura
# box()
# #Acrescentando uma grade
# grid(nx = 12)
#Acrescentando um titulo
title("Proporcao de dados faltantes")
#Acresentando uma legenda
legend(-70.86522, -18.59567, legend=names(attr(corcat, "table")),
fill=attr(corcat, "palette"), cex=0.9, bty="n", title = "Legenda:", title.adj = 0)2 método: spplot()
Criando mapas com a função spplot(), muito útil para criar mapas a partir de um arquivo shapefile e para plots de dados de área, veja:
#Carregando os pacotes RColorBrewer e sp
library(RColorBrewer)
library(sp)
# Mapas coropléticos
#Criando os intervalos de cores
intervalos=quantile(mapa$PERC_IGNORADO, probs = seq(0,1,0.125)) + c(-0.001,rep(0,7),0.001)
#Plotando o mapa com tons avermelhados
spplot(mapa,c("PERC_IGNORADO"),
at=intervalos,
col.regions =brewer.pal(8, "Reds")) #Outras opções de cores: Greens, BrBG, Accent
Para criar mapas com intervalos do mesmo tamanho, os seguitnes códigos podem ser utilizados:
#Criando os intervalos de mesmo tamanho
ampli_dados = max(mapa$PERC_IGNORADO) - min(mapa$PERC_IGNORADO)
k = 6 #número de classes + 1
ampli_classe = ampli_dados/k
intervalos2 = NULL
intervalos2[1] = min(mapa$PERC_IGNORADO)
for(i in 2:k){
intervalos2[i] = intervalos2[i-1] + ampli_classe
}Portanto, o mapa:
#Plotando o mapa com tons azulados e instervalos de mesmo tamanho:
spplot(mapa,c("PERC_IGNORADO"),at=intervalos2,col.regions =brewer.pal(5, "Blues"))
Dois mapas podem ser comparados da eguinte forma:
# Comparando dois mapas
#Criando os intervalos com base em quantis
intervalos3=quantile(c(mapa$PERC_IGNORADO,mapa$PERC_TOTAL_REGISTROS), probs = seq(0,1,0.125))
#Plotando o mapa com tons alaranjados
spplot(mapa,c("PERC_IGNORADO","PERC_TOTAL_REGISTROS"),at=intervalos3,col.regions =brewer.pal(8,
"Oranges"),names.attr = c("Percentual de ignorados","Percentual total"))
Sintaxe do ggplot
3 método: ggplot()
O comando map_data() faz a importanção automática do shapefile do Brasil, porém ele não conta com as linhas das separações dos estados, veja:
library(ggplot2)
library(maps) # For map data## Warning: package 'maps' was built under R version 3.4.2##
## Attaching package: 'maps'## The following object is masked from 'package:plyr':
##
## ozone# Obtendo o mapa do brasil
states_map <- map_data("world", region="Brazil")
#Veja como podem ser plotados os mapas via ggplot
#Plot
g1=ggplot(states_map, aes(x=long, y=lat, group=group)) +
geom_polygon(fill="white", colour="black")
# geom_path (no fill) e Mercator projection
g2=ggplot(states_map, aes(x=long, y=lat, group=group)) +
geom_path() + coord_map("mercator")
gridExtra::grid.arrange(g1,g2,ncol=2)
Com a sintaxe para o plot com o ggplot é baseada nas coordenadas de latitude e longitude, foi necessário fazer uma busca das coordenadas dos municípios do Brasil para em seguida realizar o plot do padrão de pontos, veja:
#Fazendo a leitura das bases de dados:
dadjov <- as.data.frame(readxl::read_excel("Prrenchimento_sexo_idade_ign_por_mun.xlsx"))
#Dados de latitude e longitude
latlong <- read_excel("latlong.xls")
#Removendo os tracos "-" e "'" :
dadjov$NOME=str_replace_all(dadjov$NOME, "[-]", " ")
base1=dadjov[,c("NOME","FREQ BASE PFA", "FREQ_IGN")]
base2=latlong[,c("LATITUDE", "LONGITUDE", "MUNICIPIO", "UF")]
names(base2)[3]="NOME"
#Criando a base de dados que sera utilizada:
base=left_join(base1, base2, by="NOME")
#Por fim, a base de dados utilizada sera:
head(base)## NOME FREQ BASE PFA FREQ_IGN LATITUDE LONGITUDE UF
## 1 ABADIA DE GOIAS 3292 5 -16.75 -49.43 GO
## 2 ABADIA DOS DOURADOS 1964 2 -18.48 -47.40 MG
## 3 ABADIANIA 3042 2 -16.20 -48.70 GO
## 4 ABAETE 8148 16 -19.16 -45.44 MG
## 5 ABAETETUBA 10273 28 -1.71 -48.88 PA
## 6 ABAIARA 891 5 -7.35 -39.04 CEElaborando os diferentes tipos de mapa com ggplot
#Plotando o mapa:
g1=ggplot(states_map, aes(x=long, y=lat, group=group)) +
geom_polygon(fill="white", colour="black")+
geom_point(data=base, aes(y=base$LATITUDE, x=base$LONGITUDE, colour=base$`FREQ BASE PFA`, group=base$UF))+
labs(colour="Frequencia total")+
geom_path() + coord_map("polyconic")
g2=ggplot(states_map, aes(x=long, y=lat, group=group)) +
geom_polygon(fill="white", colour="black")+
geom_point(data=base, aes(y=base$LATITUDE, x=base$LONGITUDE, colour=base$FREQ_IGN, group=base$UF))+
labs(colour="Frequencia Ignorados")+
geom_path() + coord_map("mercator")
gridExtra::grid.arrange(g1,g2,ncol=1)## Warning: Removed 171 rows containing missing values (geom_point).
## Warning: Removed 171 rows containing missing values (geom_point).
Problema
A duvida surgiu em uma tentativa de elaboraçao de mapas sobre dados dos municipios brasileiros. A maior parte dos tutoriais do ggplot2 já partem de um data.frame contendo colunas “lat” e “long”, enquanto eu tinha em mãos um .shp de polígonos.
De acordo com esta solução será possível elaborar um mapa do Brasil com o ggplot
Iniciando
Fazendo o carregamento dos pacotes necessários:
Em seguida, será feita a leitura dos dados, o shapefile pode ser baixado em http://downloads.ibge.gov.br/downloads_geociencias.html
pelo caminho:
organizacao_do_territorio/malhas_territoriais/…
mapa = readOGR("Brasil5565.shp")## OGR data source with driver: ESRI Shapefile
## Source: "Brasil5565.shp", layer: "Brasil5565"
## with 5565 features
## It has 20 fields
## Integer64 fields read as strings: gid pop2010 homens2010 mulheres20 urbana2010 rural2010 domicilios V07 V08 V09 V10 V11 V12 V13 V14 codibge6 codibge7Os objetos da classe SpatialPolygonsDataFrame são formados por partes chamadas de slots, sendo que os principais são:
- slot data: contém os atributos dos polígonos. Os mesmos que estão no arquivo .dbf
- slot polygons: contém os polígonos originalmente armazenados no arquivo .shp
Os nomes dos slots podem ser conferidos pelo comando:
slotNames(mapa)## [1] "data" "polygons" "plotOrder" "bbox" "proj4string"Vejamos a seguir, como fazer o tratamendo no arquivo shapefile para utilizar com o ggplot
Criando base com UF e REGIAO:
#Leitura do arquivo contendo o nome UF e REGIOES:
regioes <- read_excel("regioes.xlsx")
#Selecionando apenas a coluna do arquivo .shp com as UF
mapa@data=data.frame(mapa@data$UF)
names(mapa@data)=c("UF")
#Unindo o data.frame das regioes com a UF:
mapa@data=join(mapa@data,regioes,by = "UF")
#Criando a Regiao para o DF:
for(i in 1:dim(mapa@data[1])[1]){
if(mapa@data$UF[i]=="DF"){mapa@data$REGIAO[i]="DF"}
}Preparando o arquivo:
Finalmente, o arquivo SpatialPolygonsDataFrame já esta preparado para sofrer as adaptações para ser utilizado com o ggplot, veja:
Primeiro, é necessário declarar um id dos polígonos no slot $@data$.
mapa@data$id = rownames(mapa@data)Outro ajuste é necessário pelo fato de que o pacote ggplot2 é otimizado para data.frame. É necessário realizar esta conversão sem no entanto perder a localização geográfica do polígono que está sendo utilizado.
A função ggplot2::fortify realiza essa tarefa e transforma um SpatialPolygonsDataFrame em arquivo de pontos, veja:
mapa.points <- fortify(mapa,region = "id")Por último, juntaremos os atributos do slot data anterior no novo objeto mapa.points:
mapa.df <- join(mapa.points,mapa@data,by = "id")O mapa do Brasil
Finalmente, o codigo abaixo cria um mapa tematico a partir de um shapefile original com geometria em poligonos:
pal <- colorRampPalette(c("#559999", "grey80", "#BB650B"))(6)
ggplot(mapa.df) +
aes(long,lat,group=group,fill=REGIAO) +
geom_polygon() +
coord_equal()+
scale_fill_manual(values=pal)
Camadas
A partir deste mapa, podemos incluir a quantidade de camadas que quisermos, veja:
ggplot(data=mapa.df) +
geom_polygon( aes(long,lat,group=group,fill=REGIAO)) +
scale_fill_manual(values=pal)+
coord_equal()+
geom_point(data=base, aes(y=base$LATITUDE, x=base$LONGITUDE, colour=base$`FREQ BASE PFA`, group=base$UF),size=0.5)## Warning: Removed 171 rows containing missing values (geom_point).
ggplot(data=mapa.df) +
geom_polygon( aes(long,lat,group=group,fill=REGIAO)) +
coord_equal()+
scale_fill_manual(values=pal)+
geom_point(data=base, aes(y=base$LATITUDE, x=base$LONGITUDE, colour=base$FREQ_IGN, group=base$UF),size=0.5)## Warning: Removed 171 rows containing missing values (geom_point).
Criando intervalos de valores:
# Find the quantile bounds
qa <- quantile(base$`FREQ BASE PFA`, c(0, 0.2, 0.4, 0.6, 0.8, 1.0))
# Add a column of the quantile category
base$porcentagem_total <- cut(base$`FREQ BASE PFA`, qa,
labels=c("0-20%", "20-40%", "40-60%", "60-80%", "80-100%"),
include.lowest=TRUE)
base$porcentagem_ignorados <- cut(base$FREQ_IGN, qa,
labels=c("0-20%", "20-40%", "40-60%", "60-80%", "80-100%"),
include.lowest=TRUE)Portanto:
ggplot(data=mapa.df) +
geom_polygon( aes(long,lat,group=group,fill=REGIAO)) +
coord_equal()+
scale_fill_manual(values=pal)+
geom_point(data=base, aes(y=base$LATITUDE, x=base$LONGITUDE, colour=base$porcentagem_total, group=base$UF),size=0.8)## Warning: Removed 171 rows containing missing values (geom_point).
Outra forma de apresentar os dados:
ggplot(data=mapa.df) +
geom_polygon( aes(long,lat,group=group,fill=REGIAO)) +
coord_equal()+
scale_fill_manual(values=pal)+
geom_point(data=base, aes(y=base$LATITUDE, x=base$LONGITUDE, colour=base$`FREQ BASE PFA`,size=base$FREQ_IGN, group=base$UF))## Warning: Removed 745 rows containing missing values (geom_point).
4 método Gviz (googleViz)
Este método funciona interagindo como google, ao rodar o código somos direcionados para uma página html iterativa de padrão de pontos, veja:
suppressPackageStartupMessages(library(googleVis))
#Fazendo a leitura das bases de dados:
dadjov <- as.data.frame(readxl::read_excel("Prrenchimento_sexo_idade_ign_por_mun.xlsx"))
#Dados de latitude e longitude
latlong <- read_excel("latlong.xls")
#Removendo os tracos "-" e "'" :
dadjov$NOME=str_replace_all(dadjov$NOME, "[-]", " ")
base1=dadjov[,c("NOME","FREQ BASE PFA", "FREQ_IGN")]
base2=latlong[,c("LATITUDE", "LONGITUDE", "MUNICIPIO", "UF")]
names(base2)[3]="NOME"
#Criando a base de dados que sera utilizada:
base=left_join(base1, base2, by="NOME")
#Criando coluna no formato lat:long
base$latlong=str_c(base$LATITUDE,":",base$LONGITUDE)
#Organizando o dataframe para utilizar a funcao:
base=as.data.frame(cbind(latlong=base$latlong,FREQ=base$`FREQ BASE PFA`))
#Omitindo dados faltantes:
base=na.omit(base)
G <- gvisGeoChart(na.omit(base), "latlong", colorvar='FREQ',
options=list(region="BR",
width=800, height=600))
#plot(G)Resultado deste tipo de mapa:
resultado do comando plot(G)
Para mais informações, consultar o link do RDocumentacion ou ainda os exemplos do manual do pacote no CRAN
5 método RgoogleMaps()
Criando mapa estático de padrão de pontos:
Este método permite utilizar um mapa estático do googple para fazer o plot dos padrões de pontos, varias opções podem ser utilizadas, veja:
#Carregando pacotes
library(rworldmap)## Warning: package 'rworldmap' was built under R version 3.4.2## ### Welcome to rworldmap ##### For a short introduction type : vignette('rworldmap')library(RgoogleMaps)## Warning: package 'RgoogleMaps' was built under R version 3.4.2#Fazendo a leitura das bases de dados:
dadjov <- as.data.frame(readxl::read_excel("Prrenchimento_sexo_idade_ign_por_mun.xlsx"))
#Dados de latitude e longitude
latlong <- read_excel("latlong.xls")
#Removendo os tracos "-" e "'" :
dadjov$NOME=str_replace_all(dadjov$NOME, "[-]", " ")
base1=dadjov[,c("NOME","FREQ BASE PFA", "FREQ_IGN")]
base2=latlong[,c("LATITUDE", "LONGITUDE", "MUNICIPIO", "UF")]
names(base2)[3]="NOME"
#Criando a base de dados que sera utilizada:
base=left_join(base1, base2, by="NOME")
#Com a base de dados criada:
base=na.omit(base)
coordinates(base) <- c("LONGITUDE", "LATITUDE")
# Definir projeção espacial
crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84") # geographical, datum WGS84
proj4string(base) <- crs.geo
# Mapa estático usando RgoogleMaps
coords <- as.data.frame(coordinates(base))
#Definindo o centro do mapa que será requerido como o centróide das coordenadas geográficas:
mapa_plot <- GetMap(center = c(mean(coords$LATITUDE), mean(coords$LONGITUDE)), zoom = 4, maptype = "mobile") # "roadmap", "sat ellite", "terrain", "hybrid", "mapmaker-roadmap", "mapmaker-hybrid"
#PLot:
PlotOnStaticMap(mapa_plot, lat = coords$LATITUDE, lon = coords$LONGITUDE, cex
= 0.5, pch = 19, col = "blue", FUN = points)
Criando mapa interativo com plotGoogleMaps():
Por fim, algumas tentativas (ainda sem sucesso) de fazer o plot dessas informaçõe sdiretamente no site do googlemaps
library(leaflet)
#Fazendo a leitura das bases de dados:
dadjov <- as.data.frame(readxl::read_excel("Prrenchimento_sexo_idade_ign_por_mun.xlsx"))
#Dados de latitude e longitude
latlong <- read_excel("latlong.xls")
#Removendo os tracos "-" e "'" :
dadjov$NOME=str_replace_all(dadjov$NOME, "[-]", " ")
base1=dadjov[,c("NOME","FREQ BASE PFA", "FREQ_IGN")]
base2=latlong[,c("LATITUDE", "LONGITUDE", "MUNICIPIO", "UF")]
names(base2)[3]="NOME"
#Criando a base de dados que sera utilizada:
base=left_join(base1, base2, by="NOME")
#Com a base de dados criada:
base=na.omit(base)
coordinates(base) <- ~LONGITUDE+LATITUDE
proj4string(base) <- CRS('+init=epsg:28992')
#Mapa 1
plot(base)#Mapa 2
#Carregabndo o pacote plotGoogleMaps
library(plotGoogleMaps)
library(RColorBrewer)
#Fazendo a leitura das bases de dados:
dadjov <- as.data.frame(readxl::read_excel("Prrenchimento_sexo_idade_ign_por_mun.xlsx"))
#Dados de latitude e longitude
latlong <- read_excel("latlong.xls")
#Removendo os tracos "-" e "'" :
dadjov$NOME=str_replace_all(dadjov$NOME, "[-]", " ")
base1=dadjov[,c("NOME","FREQ BASE PFA", "FREQ_IGN")]
base2=latlong[,c("LATITUDE", "LONGITUDE", "MUNICIPIO", "UF")]
names(base2)[3]="NOME"
#Criando a base de dados que sera utilizada:
base=left_join(base1, base2, by="NOME")
base=na.omit(base)
coordinates(base)<-c("LATITUDE","LONGITUDE")
proj4string(base) <- CRS("+init=epsg:28992")
mapa2 <- plotGoogleMaps(base,
filename='mapa2.htm',
zcol = "FREQ_IGN",
colPalette= brewer.pal(5, "Reds"),
mapTypeId = "TERRAIN")
map <- plotGoogleMaps(base,
zcol = "FREQ_IGN",
zoom = 8,
filename = "Map_GoogleMaps.html",
layerName = "Proporcao de dados faltantes",
colPalette = brewer.pal(5, "Oranges"),
mapTypeId = "ROADMAP") #Outras opções de mapTypeId: HYBRID,SATELLITE, TERRAIN,ROADMAP
bubbleGoogleMaps(base,zcol='FREQ_IGN', max.radius = 80,
filename='myMap3.htm')
meuse_proj = base[,c(base$)]
proj4string(meuse_proj) <- CRS("+init=epsg:28992")
mapa2 <- plotGoogleMaps(meuse_proj, filename='zinco_mapa2.htm', zcol = "zinc", colPalette
= brewer.pal(5, "Reds"))