Dados de origem

Alguns desses dados de origem podem ser diretamente obtidos desde R

Pacotes

library(maps) #mapas simples, eixos, escala, cidades 
library(mapdata) #base de dados WorldHires e rios
library(rworldmap) #outra base de dados de mapas do mundo
library(maptools) #Ler ESRI shapefiles 
library(mapproj) #Projeções e grids
library(ggmap) #Gmaps, OSM + mapas baseados em ggplot2
library(rgdal)

Mapas desde os pacotes de R

library(maps)
par(mar=c(1,1,1,1))
map("worldHires","Brazil")

Mapas desde os pacotes de R

par(mar=c(1,1,1,1))
map("world","Brazil")
map.axes()

Mapas desde os pacotes de R

map("world","Brazil")
map.axes()
map.scale(ratio = F, cex = 0.7) #tentem ratio = T

Mapas desde os pacotes de R

Dá para usar qualquer elemento de plot do base R (points, abline, text, legend) e os parâmetros de par() como pch (o símbolo), cex (o tamanho do símbolo), lty, lwd (tipo e largura de linha), font (1= normal, 2= itálica, 3= bold).

map("world","Brazil")
map.axes()
map.scale(ratio = F, cex = 0.7)
abline(h = 0, lty = 2)

Mapas desde os pacotes de R

map("world","Brazil", fill=T, col="grey90")
map(,,add=T)
map.axes()
map.scale(ratio=F, cex=0.7)
abline(h=0, lty = 2)
map.cities(country = "Brazil",minpop = 2000000,pch=19, cex=1.2)# pacote maps

Mapas desde os pacotes de R

m <- map("world","Brazil", fill=T, col="grey95")
map(,,add=T)
map.axes()
map.scale(ratio=F, cex=0.7)
abline(h=0, lty = 2)
map.grid(m, nx = 5, ny = 5, col="grey50", font=1, cex=0.7 , pretty = T)#library(mapproj)

Projeções

R suporta diferentes tipos de projeção. Pacote rgdal e base de dados epsg, que é um padrão para sistemas de referência espacial. http://spatialreference.org, https://github.com/OSGeo/proj.4/wiki

m <- map(,"brazil", plot=FALSE)
map(,"brazil", project="albers", par=c(-2, -22),lwd=2)
map.grid(m,lwd=2,font=1,cex=0.7, col="grey70") 
legend(-0.3, 4.37, legend= c("Albers", "WGS84"), col=c("black","grey70"), lty=c(1,2),lwd=2,horiz = T,bty = "o",title = "Projeção",box.col = "white")

Ler shapefiles

#library(rgdal)
# maptools::readShapePoly
br <- readShapePoly("./shapefiles/brasil/BRA_adm_shp/BRA_adm0.shp") # 0: país, 1: estados, 2: municípios
estados1 <- readShapePoly("./shapefiles/brasil/BRA_adm_shp/BRA_adm1.shp") # 0: país, 1: estados, 2: municípios
cidades <- readShapePoly("./shapefiles/brasil/BRA_adm_shp/BRA_adm2")
biomas <- readShapePoly("./shapefiles/brasil/BR_BIOMAS_IBGE.shp")

# rgdal::readOGR
estados <- readOGR("./shapefiles/brasil/BRA_adm_shp/", "BRA_adm1")
## OGR data source with driver: ESRI Shapefile 
## Source: "./shapefiles/brasil/BRA_adm_shp/", layer: "BRA_adm1"
## with 27 features
## It has 12 fields

readOGR inclui uma projeção padrão no shapefile - não será necessário incluí-la manualmente depois, caso seja necessário

proj4string(estados1)
## [1] NA
proj4string(estados)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

Subsetar shapefiles

names(cidades)
##  [1] "ID_0"      "ISO"       "NAME_0"    "ID_1"      "NAME_1"   
##  [6] "ID_2"      "NAME_2"    "HASC_2"    "CCN_2"     "CCA_2"    
## [11] "TYPE_2"    "ENGTYPE_2" "NL_NAME_2" "VARNAME_2"
head(cidades$NAME_1)
## [1] Acre Acre Acre Acre Acre Acre
## 27 Levels: Acre Alagoas Amapá Amazonas Bahia Ceará ... Tocantins
RJ <- cidades[cidades@data$NAME_1=="Rio de Janeiro",] ##seleciona e cria um novo shapefile

plot(RJ)
map.axes()
map.scale(ratio=F)
plot(RJ[RJ$NAME_2=="Silva Jardim",], add=T, col="grey70") #seleciona mas não cria um novo shapefile, apenas na hora de plotar
plot(br, add=T, border=alpha("grey50",0.3))

Projeções

Exemplo com a Reserva Biológica de Poço das Antas: dados obtidos do Plano de Manejo estão em UTM e fragmentos Florestais do SOS Mata Atläntica www.sosma.org.br

PDA <- readOGR("./shapefiles/rbpda/","LimitedaRB")
## OGR data source with driver: ESRI Shapefile 
## Source: "./shapefiles/rbpda/", layer: "LimitedaRB"
## with 1 features
## It has 2 fields
br101 <- readOGR("./shapefiles/rbpda/","BR_101")
## OGR data source with driver: ESRI Shapefile 
## Source: "./shapefiles/rbpda/", layer: "BR_101"
## with 1 features
## It has 6 fields
saojoao <- readOGR("./shapefiles/rbpda/","riosaojoao")
## OGR data source with driver: ESRI Shapefile 
## Source: "./shapefiles/rbpda/", layer: "riosaojoao"
## with 1 features
## It has 1 fields
SOSma <- readOGR(dsn = "./shapefiles/brasil/SOSMA/","rema_fd4d801731725513a4d77aa9bb35534b5305") 
## OGR data source with driver: ESRI Shapefile 
## Source: "./shapefiles/brasil/SOSMA/", layer: "rema_fd4d801731725513a4d77aa9bb35534b5305"
## with 200 features
## It has 2 fields
plot(PDA)
plot(br101, add=T, col="grey")
plot(saojoao, col="dodgerblue", add=T)
map.axes()
plot(SOSma, add=T) #### NAO FUNCIONA PORQUE ESTA EM OUTRA PROJECAO!

#map.scale() também não fica bom. 

proj4string(SOSma) 
## [1] "+proj=poly +lat_0=0 +lon_0=-45 +x_0=0 +y_0=0 +ellps=aust_SA +units=m +no_defs"
plot(SOSma)
map.axes()#SOS está tudo em outra projeção

Vamos reprojetar tudo para WGS84.

O padrão brasileiro atual é sirgas2000 que é teoricamente igual a WGS84

proj4string(PDA)
## [1] NA
proj4string(PDA) <- CRS("+init=epsg:29193")## isto é UTM em SAD69 brasileiro
proj4string(PDA)
## [1] "+init=epsg:29193 +proj=utm +zone=23 +south +ellps=aust_SA +towgs84=-66.87,4.37,-38.52,0,0,0,0 +units=m +no_defs"
PDA.wsg84 <- spTransform(PDA, CRS("+proj=longlat +datum=WGS84"))
SOS.wgs84 <- spTransform(SOSma, CRS("+proj=longlat +datum=WGS84"))

plot(PDA.wsg84, add=F, border="red", lwd=2)
plot(SOS.wgs84, add=T, col=alpha("green", 0.5), border=alpha("white",0))

#plot(RJ, add=T, names=T) ##nao fica legal


plot(SOS.wgs84, add=F, col=alpha("green", 0.5), border=alpha("white",0))
plot(PDA.wsg84, add=T, border="red", lwd=2)
plot(RJ, add=T) 
class(RJ)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
pts.RJ <- getSpatialPolygonsLabelPoints(RJ)
text(pts.RJ@coords, labels=RJ$NAME_2)
map.axes()
map.scale(y=-22.78,ratio=F,cex=0.7, font=2)

Exercicío: a partir dos shapefiles de Estados e Biomas, tentar reproduzir um mapa similar a este:

ex

Outro framework para fazer mapas em R: ggplot2

source("./fct/create_scale_bar.R")
brasil.mapa <- borders("worldHires", regions = "Brazil", fill = "grey90", colour = "black")
bra.map <- map_data("worldHires","Brazil")

mapa.1 <- ggplot() + brasil.mapa + coord_equal()  ##borders inclui o mapa como uma camada no fundo 
mapa.1

mapa.2 <- ggplot() + 
    geom_polygon(data = bra.map, aes(x = long, y = lat, group = group)) +
    coord_equal() # map_data cria um dataframe que deve ser acrescentado como geom_polygon()
mapa.2

Os objetos espaciais devem passar por fortify() para poder ser lidos por ggplot2. A lógica de plotar funciona igual que os gráficos de ggplot, somando camadas sequencialmente.

MA <- biomas[biomas$CD_LEGENDA=="MATA ATLANTICA",]

MAf <- fortify(MA)
PDAf <- fortify (PDA.wsg84)

mapa.2 <- mapa.1 +
    scaleBar(lon = -40, lat = -30, distanceLon = 500, distanceLat = 70, distanceLegend = 150, dist.unit = "km", orientation = FALSE) +
    geom_polygon(data = MAf, aes(x = long, y = lat, group = group), fill = "darkgreen")
mapa.2

Incluir camadas de google e open street maps

library(ggmap)
 pda <- get_map("Reserva Biológica de Poço das Antas",zoom=12)
 pda.map <-   ggmap(pda, extent= "panel")
pda.map

#pda <- get_map("Reserva Biológica de Poço das Antas",zoom=12, source= "stamen",maptype="toner")
# pda.toner <-   ggmap(pda, extent= "panel")
#pda.toner

 #pda <- get_map("Reserva Biológica de Poço das Antas",zoom=12, source = "osm")
 #pda.osm <-   ggmap(pda, extent= "panel")
#pda.osm

mapa.test <- pda.map +
      geom_polygon(data = PDAf, aes(x = long, y = lat, group = group),
        fill = alpha("white",0.2), color = "black", size = 0.5) 

mapa.test          

#ggsave(filename = "map1.png", device="png")

Ler arquivos raster

O pacote raster faz a maior parte das operações necessárias para ler, cortar e escrever rasters em vários formatos.

library(raster)
alt <- raster("./rasters/h_dem/hdr.adf")
plot(alt)

Ler arquivos raster

Vários rasters podem ser juntos num “stack”.

arquivos <- list.files("./rasters", full.names = T)
head(arquivos)
## [1] "./rasters/aet_yr" "./rasters/alpha"  "./rasters/bio_1" 
## [4] "./rasters/bio_10" "./rasters/bio_11" "./rasters/bio_12"
env <- stack(arquivos)

Ler arquivos raster

Cortar arquivos raster

MA <- biomas[biomas$CD_LEGENDA=="MATA ATLANTICA",]
plot(MA)

#### CROP THEN MASK
alt.crop <- crop(x = alt,y = MA)
plot(alt.crop)

alt.crop <- mask(alt.crop,MA)
plot(alt.crop)
plot(MA,add=T)

Arquivos binários

plot(alt.crop)

plot(alt.crop>0, legend= F)

Escrever arquivos raster

dir.create("./MataAtlantica")
writeRaster(alt.crop,filename = "./MataAtlantica/alt", format= "GTiff",overwrite=T)
### depois dá para ler diretamente estes arquivos ----
alt.mat <- raster("./MataAtlantica/alt.tif")
plot(alt.mat)

Acrescentando dados própios

caryocar <- read.delim("./occurrence data/caryocar.txt")
head(caryocar)
##                     sp       lon       lat
## 1 Caryocar brasiliense -49.24670 -20.55282
## 2 Caryocar brasiliense -47.51667 -14.11667
## 3 Caryocar brasiliense -46.25000 -20.91667
## 4 Caryocar brasiliense -48.91667 -10.08333
## 5 Caryocar brasiliense -51.76667 -12.81667
## 6 Caryocar brasiliense -48.96043 -22.32330
plot(alt>0,legend=F, col="grey")
points(caryocar$lon, caryocar$lat, pch=19,cex=0.5)

library(spatstat)
ch <- convexhull.xy(cbind(caryocar$lon, caryocar$lat))
plot(ch, add=T, col=alpha("green",0.5))

Registros de ocorrência a partir de GBIF

Global Biodiversity Information Facility, www.gbif.org

library(rgbif)
library(sp)
euterpe <- occ_data(scientificName = "Euterpe edulis",hasCoordinate = T)
euterpe <- SpatialPoints(cbind(euterpe$data$decimalLongitude,euterpe$data$decimalLatitude))
if(!dir.exists("occs")){
dir.create("./occs")
write.table(euterpe,"./occs/euterpe.txt",  row.names = F)
}

Registros de ocorrência a partir de GBIF

euterpe2 <- read.delim("./occs/euterpe.txt",sep="")
plot(MA)
points(euterpe2[,c(1,2)],pch=19,col="red",cex=0.8)
map.axes()
map.scale(ratio=F,relwidth = 0.3, x=-45)

Informação adicional