A Annexo 1
Exemplos de código R usado para - baixar Mapbiomas, - recortar para uma região de interesse - plotar com legenda
A.1 Download mapbiomas cover
Site: https://brasil.mapbiomas.org/colecoes-mapbiomas/. Any zeros are NA (NODATA). http://forum.mapbiomas.ecostage.com.br/t/pixels-com-valor-zero/170/5 . And “27” “Não observado”. No exemplo os arquivos vai ficar no working directory.
#Get mapbiomas data
#
url_text <- "https://storage.googleapis.com/mapbiomas-public/initiatives/brasil/collection_8/lclu/coverage/"
layer_names <- paste("brasil_coverage_",1985:2022, sep="")
file_type <- ".tif"
filenames_mapbioma_8 <- paste(layer_names, file_type, sep="")
urlnames_mapbioma_8 <- paste(url_text, layer_names, file_type, sep="")
#set timeout to 20 minutes as big files and teeny tiny internet connection
options(timeout = max((60*20), getOption("timeout")))
# test with one year - 1985
res <- utils::download.file(url=urlnames_mapbioma_8[1],
destfile=filenames_mapbioma_8[1],
method="auto", quiet = FALSE, mode = "wb",
cacheOK = TRUE)
#make list
dflist <- data.frame(layer_names = layer_names,
urls = urlnames_mapbioma_8,
filenames = filenames_mapbioma_8)
alist <- split(dflist, f = dflist$layer_names)
#make function to automatically process list elements
get_mapbiomas <- function(x){
res <- utils::download.file(url=x$urls,
destfile=x$filenames,
method="auto", quiet = FALSE, mode = "wb",
cacheOK = TRUE)
res
}
#run function in blocks
lapply(alist[2], get_mapbiomas)
lapply(alist[3:5], get_mapbiomas)
lapply(alist[6:8], get_mapbiomas)
lapply(alist[9:20], get_mapbiomas)
lapply(alist[21:30], get_mapbiomas)
lapply(alist[31:37], get_mapbiomas)A.2 Crop mapbiomas cover
Identify all tif files in a folder.
#works one folder at a time
get_files <- function(folder_name = NA) {
library(tidyverse)
folder_location <- folder_name
in_files <- list.files(folder_location,
pattern = ".tif", full.names = TRUE)
data.frame(folder_id = folder_location, file_id = in_files) |>
group_by(folder_id, file_id) |>
dplyr::summarise(file_count = dplyr::n()) |>
ungroup() -> df_muni_tif
return(df_muni_tif)
}
infolder <- "E:/mapbiomas"
df_muni_tif <- get_files(folder_name = infolder)
#update
df_muni_tif[1:38, ] |>
dplyr::mutate(state_code = str_sub(folder_id, -2, -1),
) -> df_muni_tif
df_muni_tif$state_code <- "AP"Area of interest - format extent for use by terra.
#extent from 50 km arouund river upstream of cachoeira caldeirao
meuSIG <- file.choose()
# "C:\\Users\\user\\Documents\\CA\\CA_2023\\EP\\vectors\\AP_rodoviario\\br156.gpkg"
#"C:\\Users\\Darren\\Documents\\gisdata\\vector\\rivers.GPKG"
# "C:\\Users\\user\\Documents\\Articles\\gis_layers\\gisdata\\inst\\vector\\rivers.gpkg"
rsl <- st_read(meuSIG, layer = "centerline")
rsl_50km <- st_union(st_buffer(rsl, dist=50000))
br156 <- st_read(meuSIG, layer = "extent_br156")
myexent <- ext(vect(rsl_50km))
myexent_jari <- vect(br156[1, ])
myexent_portogrande <- vect(br156[2, ])Now to crop and project to desired coordinate system (epsg:31976).
# This function reads in file, crops and projects to new coordinate system.
# Then exports resutl as Geotiff (.tif).
crop_proj <- function(x, myextent = NA,
folder_name = NA, file_name = NA,
state_id = NA, sf_state = NA) {
library(plyr)
library(tidyverse)
library(terra)
library(sf)
library(stringi)
state_sigla <- x$state_code
rin <- x$file_id
#rin <- df_muni_tif$file_id[1]
rbig <- terra::rast(rin)
layer_name <- names(rbig)
layer_year <- stri_sub(layer_name,-4,-1)
out_name <- paste("mapbiomas_cover", layer_year, sep="_")
e2 <- ext(project(myextent, rbig))
#Crop
rtest_mask <- crop(rbig, e2, snap="out")
rm("rbig")
#Project
new_crs_utm <- "epsg:31976"
rtest_mask <- project(rtest_mask, new_crs_utm, method = "near")
names(rtest_mask) <- out_name
#Export
outfile_name <- paste(file_name,layer_year, ".tif", sep = "")
f <- file.path(folder_name, outfile_name)
writeRaster(rtest_mask, f, datatype = "INT1U", NAflag=27,
gdal=c("COMPRESS=DEFLATE"), overwrite = TRUE)
#clear temporary files
tmpFiles(current =TRUE, remove = TRUE)
endtime <- Sys.time()
textout <- paste(outfile_name, ": ", endtime, sep="")
print(textout)
}
#run br156 jari
plyr::a_ply(df_muni_tif[2:38, ], .margins = 1,
.fun = crop_proj, myextent = myexent_jari,
folder_name = "E:/mapbiomas/br156_jari",
file_name = "br156_jari_")
#run br156 Porto Grande
plyr::a_ply(df_muni_tif[2:38, ], .margins = 1,
.fun = crop_proj, myextent = myexent_portogrande,
folder_name = "E:/mapbiomas/br156_portogrande",
file_name = "br156_portogrande_")A.3 Plot with legend
Combine all years into multi layer raster.
# check
rtest <- rast("E:/mapbiomas/br156_jari/br156_jari_1985.tif")
rtest <- rast("E:/mapbiomas/br156_portogrande/br156_portogrande_1985.tif")
plot(rtest)
infolder_rio <- "E:/mapbiomas/AP_utm_rio"
df_tif_rio <- get_files(folder_name = infolder_rio)
r85a22 <- rast(df_tif_rio$file_id)
writeRaster(r85a22, "E:/mapbiomas/utm_cover_AP_rio_85a22.tif",
datatype = "INT1U", NAflag=27,
gdal=c("COMPRESS=DEFLATE"), overwrite = TRUE)
# br156 jari
infolder_jari <- "E:/mapbiomas/br156_jari"
df_tif_jari <- get_files(folder_name = infolder_jari)
r85a22 <- rast(df_tif_jari$file_id)
writeRaster(r85a22, "E:/mapbiomas/br156_jari_85a22.tif",
datatype = "INT1U", NAflag=27,
gdal=c("COMPRESS=DEFLATE"), overwrite = TRUE)
# br156 Porto Grande
infolder_pg <- "E:/mapbiomas/br156_portogrande"
df_tif_pg <- get_files(folder_name = infolder_pg)
r85a22 <- rast(df_tif_pg$file_id)
writeRaster(r85a22, "br156_portogrande_85a22.tif",
datatype = "INT1U", NAflag=27,
gdal=c("COMPRESS=DEFLATE"), overwrite = TRUE)
# Garimpo Lorenco
infolder_lor <- "E:/mapbiomas/garimpo_lorenco"
df_tif_lor <- get_files(folder_name = infolder_lor)
r85a22 <- rast(df_tif_lor$file_id)
writeRaster(r85a22, "E:/mapbiomas/lorenco_85a22.tif",
datatype = "INT1U", NAflag=27,
gdal=c("COMPRESS=DEFLATE"), overwrite = TRUE)
# Garimpo Vila Nova
infolder_vn <- "E:/mapbiomas/garimpo_vilanova"
df_tif_vn <- get_files(folder_name = infolder_vn)
r85a22 <- rast(df_tif_vn$file_id)
writeRaster(r85a22, "E:/mapbiomas/vilanova_85a22.tif",
datatype = "INT1U", NAflag=27,
gdal=c("COMPRESS=DEFLATE"), overwrite = TRUE)Sample for BR156. Processed in QGIS to get same extent. Downloaded from IBGE, Bases cartográficas contínuas Brasil - 1:250 000, versão 2023: https://www.ibge.gov.br/geociencias/downloads-geociencias.html?caminho=cartas_e_mapas/bases_cartograficas_continuas/bc250/versao2023/
Need buffers and sample points.
vin <- "C:\\Users\\user\\Documents\\CA\\CA_2023\\EP\\vectors\\AP_rodoviario\\br156.gpkg"
st_layers(vin)## Driver: GPKG
## Available layers:
## layer_name geometry_type features fields crs_name
## 1 extent_br156 Multi Polygon 2 7 WGS 84
## 2 BR156_jari_portogrande Multi Line String 35 21 SIRGAS 2000## Reading layer `extent_br156' from data source
## `C:\Users\user\Documents\CA\CA_2023\EP\vectors\AP_rodoviario\br156.gpkg'
## using driver `GPKG'
## Simple feature collection with 2 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -52.50101 ymin: -0.8184735 xmax: -50.67069 ymax: 1.540734
## Geodetic CRS: WGS 84## Reading layer `BR156_jari_portogrande' from data source
## `C:\Users\user\Documents\CA\CA_2023\EP\vectors\AP_rodoviario\br156.gpkg'
## using driver `GPKG'
## Simple feature collection with 35 features and 21 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -52.50101 ymin: -0.8184735 xmax: -50.88881 ymax: 1.539533
## Geodetic CRS: SIRGAS 2000br156_j <- br156 |> filter(layer == "BR156_mazagao_jari")
br156_p <- br156 |> filter(layer != "BR156_mazagao_jari")
ext_j <- extents |> filter(trecho == "Jari_Mazagao")
ext_p <- extents |> filter(trecho != "Jari_Mazagao")
# buffers roads
j_b250 <- st_buffer(st_union(br156_j), dist = 250 ) |>
st_as_sf() |> mutate(buff_dist = 250)
j_b500 <- st_buffer(st_union(br156_j), dist = 500 ) |>
st_as_sf() |> mutate(buff_dist = 500)
j_b1000 <- st_buffer(st_union(br156_j), dist = 1000 ) |>
st_as_sf() |> mutate(buff_dist = 1000)
j_b2000 <- st_buffer(st_union(br156_j), dist = 2000 ) |>
st_as_sf() |> mutate(buff_dist = 2000)
j_b4000 <- st_buffer(st_union(br156_j), dist = 4000 ) |>
st_as_sf() |> mutate(buff_dist = 4000)
j_b8000 <- st_buffer(st_union(br156_j), dist = 8000 ) |>
st_as_sf() |> mutate(buff_dist = 8000)
p_b250 <- st_buffer(st_union(br156_p), dist = 250 ) |>
st_as_sf() |> mutate(buff_dist = 250)
p_b500 <- st_buffer(st_union(br156_p), dist = 500 ) |>
st_as_sf() |> mutate(buff_dist = 500)
p_b1000 <- st_buffer(st_union(br156_p), dist = 1000 ) |>
st_as_sf() |> mutate(buff_dist = 1000)
p_b2000 <- st_buffer(st_union(br156_p), dist = 2000 ) |>
st_as_sf() |> mutate(buff_dist = 2000)
p_b4000 <- st_buffer(st_union(br156_p), dist = 4000 ) |>
st_as_sf() |> mutate(buff_dist = 4000)
p_b8000 <- st_buffer(st_union(br156_p), dist = 8000 ) |>
st_as_sf() |> mutate(buff_dist = 8000)
# bind together
jari_buff <- bind_rows(j_b250, j_b500, j_b1000, j_b2000, j_b4000, j_b8000) |>
mutate(trecho = "Jari") |> st_crop(ext_j)## Warning: attribute variables are assumed to be spatially constant throughout all geometriesporto_buff <- bind_rows(p_b250, p_b500, p_b1000, p_b2000, p_b4000, p_b8000) |>
mutate(trecho = "Porto Grande") |> st_crop(ext_p)## Warning: attribute variables are assumed to be spatially constant throughout all geometriesroad_buff <- bind_rows(jari_buff, porto_buff)
# sample points density of one every 3 km
j_line <- br156_j |>
st_sf() |>
st_combine() |>
st_sf() |>
st_line_merge() |>
st_cast("LINESTRING")
# 49
j_points <- st_line_sample(j_line,
density = units::set_units(0.3, 1/km)) |>
st_cast("POINT") |> data.frame() |> st_as_sf() |>
mutate(trecho = "Jari")
#Porto Grande
p_line <- br156_p |>
st_sf() |>
st_combine() |>
st_sf() |>
st_line_merge() |>
st_cast("LINESTRING")
# 37
p_points <- st_line_sample(p_line,
density = units::set_units(0.3, 1/km)) |>
st_cast("POINT") |> data.frame() |> st_as_sf() |>
mutate(trecho = "Porto Grande")
# Point buffers
br156_points <- bind_rows(j_points, p_points)
br156_points$aid <- rank(-st_coordinates(br156_points)[,2])
# length(unique(br156_points$aid)) # 86
points_b250 <- st_buffer(br156_points, dist = 250 ) |>
st_as_sf() |> mutate(buff_dist = 250)
points_b500 <- st_buffer(br156_points, dist = 500 ) |>
st_as_sf() |> mutate(buff_dist = 500)
points_b1000 <- st_buffer(br156_points, dist = 1000 ) |>
st_as_sf() |> mutate(buff_dist = 1000)
points_b2000 <- st_buffer(br156_points, dist = 2000 ) |>
st_as_sf() |> mutate(buff_dist = 2000)
points_b4000 <- st_buffer(br156_points, dist = 4000 ) |>
st_as_sf() |> mutate(buff_dist = 4000)
points_b8000 <- st_buffer(br156_points, dist = 8000 ) |>
st_as_sf() |> mutate(buff_dist = 8000)
# bind
point_buff_all <- bind_rows(points_b250, points_b500, points_b1000,
points_b2000, points_b4000, points_b8000) |>
mutate(aid = paste(trecho, aid, buff_dist, sep="_")) |>
st_crop(extents)## Warning: attribute variables are assumed to be spatially constant throughout all geometries## Warning: attribute variables are assumed to be spatially constant throughout all geometries## Warning: attribute variables are assumed to be spatially constant throughout all geometriespoint_buff_all <- bind_rows(point_buff_j, point_buff_p)
# Export
outfile <- "C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg"
# roads
st_write(extents, dsn = outfile,
layer = "extents_br156", delete_layer = TRUE, append = TRUE)## Deleting layer `extents_br156' using driver `GPKG'
## Updating layer `extents_br156' to data source `C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg' using driver `GPKG'
## Writing 2 features with 7 fields and geometry type Multi Polygon.## Deleting layer `extent_jari' using driver `GPKG'
## Updating layer `extent_jari' to data source `C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg' using driver `GPKG'
## Writing 1 features with 7 fields and geometry type Multi Polygon.## Deleting layer `extent_porto' using driver `GPKG'
## Updating layer `extent_porto' to data source `C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg' using driver `GPKG'
## Writing 1 features with 7 fields and geometry type Multi Polygon.st_write(br156, dsn = outfile,
layer = "br156_jari_portogrande", delete_layer = TRUE, append = TRUE)## Deleting layer `br156_jari_portogrande' using driver `GPKG'
## Updating layer `br156_jari_portogrande' to data source `C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg' using driver `GPKG'
## Writing 35 features with 21 fields and geometry type Multi Line String.## Deleting layer `br156_jari_mazagao' using driver `GPKG'
## Updating layer `br156_jari_mazagao' to data source `C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg' using driver `GPKG'
## Writing 1 features with 0 fields and geometry type Line String.st_write(p_line, dsn = outfile,
layer = "br156_portogrande_itaubal", delete_layer = TRUE, append = TRUE)## Deleting layer `br156_portogrande_itaubal' using driver `GPKG'
## Updating layer `br156_portogrande_itaubal' to data source `C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg' using driver `GPKG'
## Writing 1 features with 0 fields and geometry type Line String.## Deleting layer `br156_buffers' using driver `GPKG'
## Updating layer `br156_buffers' to data source `C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg' using driver `GPKG'
## Writing 12 features with 2 fields and geometry type Polygon.st_write(jari_buff, dsn = outfile,
layer = "br156_buffers_jari", delete_layer = TRUE, append = TRUE)## Deleting layer `br156_buffers_jari' using driver `GPKG'
## Updating layer `br156_buffers_jari' to data source `C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg' using driver `GPKG'
## Writing 6 features with 2 fields and geometry type Polygon.st_write(porto_buff, dsn = outfile,
layer = "br156_buffers_portogrande", delete_layer = TRUE, append = TRUE)## Deleting layer `br156_buffers_portogrande' using driver `GPKG'
## Updating layer `br156_buffers_portogrande' to data source `C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg' using driver `GPKG'
## Writing 6 features with 2 fields and geometry type Polygon.## Deleting layer `br156_points' using driver `GPKG'
## Updating layer `br156_points' to data source `C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg' using driver `GPKG'
## Writing 86 features with 2 fields and geometry type Point.## Deleting layer `br156_points_jari' using driver `GPKG'
## Updating layer `br156_points_jari' to data source `C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg' using driver `GPKG'
## Writing 49 features with 1 fields and geometry type Point.st_write(p_points, dsn = outfile,
layer = "br156_points_portogrande", delete_layer = TRUE, append = TRUE)## Deleting layer `br156_points_portogrande' using driver `GPKG'
## Updating layer `br156_points_portogrande' to data source `C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg' using driver `GPKG'
## Writing 37 features with 1 fields and geometry type Point.st_write(point_buff_all, dsn = outfile,
layer = "br156_points_buffers", delete_layer = TRUE, append = TRUE)## Deleting layer `br156_points_buffers' using driver `GPKG'
## Updating layer `br156_points_buffers' to data source `C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg' using driver `GPKG'
## Writing 516 features with 3 fields and geometry type Polygon.st_write(point_buff_j, dsn = outfile,
layer = "br156_points_buffers_jari", delete_layer = TRUE, append = TRUE)## Deleting layer `br156_points_buffers_jari' using driver `GPKG'
## Updating layer `br156_points_buffers_jari' to data source `C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg' using driver `GPKG'
## Writing 294 features with 3 fields and geometry type Polygon.st_write(point_buff_p, dsn = outfile,
layer = "br156_points_buffers_porto", delete_layer = TRUE, append = TRUE)## Deleting layer `br156_points_buffers_porto' using driver `GPKG'
## Updating layer `br156_points_buffers_porto' to data source `C:/Users/user/Documents/Articles/gis_layers/gisdata/inst/vector/br156.gpkg' using driver `GPKG'
## Writing 222 features with 3 fields and geometry type Polygon.## Driver: GPKG
## Available layers:
## layer_name geometry_type features fields crs_name
## 1 extents_br156 Multi Polygon 2 7 SIRGAS 2000 / UTM zone 22N
## 2 extent_jari Multi Polygon 1 7 SIRGAS 2000 / UTM zone 22N
## 3 extent_porto Multi Polygon 1 7 SIRGAS 2000 / UTM zone 22N
## 4 br156_jari_portogrande Multi Line String 35 21 SIRGAS 2000 / UTM zone 22N
## 5 br156_jari_mazagao Line String 1 0 SIRGAS 2000 / UTM zone 22N
## 6 br156_portogrande_itaubal Line String 1 0 SIRGAS 2000 / UTM zone 22N
## 7 br156_buffers Polygon 12 2 SIRGAS 2000 / UTM zone 22N
## 8 br156_buffers_jari Polygon 6 2 SIRGAS 2000 / UTM zone 22N
## 9 br156_buffers_portogrande Polygon 6 2 SIRGAS 2000 / UTM zone 22N
## 10 br156_points Point 86 2 SIRGAS 2000 / UTM zone 22N
## 11 br156_points_jari Point 49 1 SIRGAS 2000 / UTM zone 22N
## 12 br156_points_portogrande Point 37 1 SIRGAS 2000 / UTM zone 22N
## 13 br156_points_buffers Polygon 516 3 SIRGAS 2000 / UTM zone 22N
## 14 br156_points_buffers_jari Polygon 294 3 SIRGAS 2000 / UTM zone 22N
## 15 br156_points_buffers_porto Polygon 222 3 SIRGAS 2000 / UTM zone 22NLegend.
lf <- "data/raster/mapbiomas/legend/Codigos-da-legenda-colecao-8.csv"
mapbiomas8_legenda <- read.csv2(lf, as.is = TRUE)Plot.
class_vals <- mapbiomas8_legenda$Class_ID
class_color <- mapbiomas8_legenda$Color
names(class_color) <- mapbiomas8_legenda$Class_ID
#labels
my_label <- mapbiomas8_legenda$Descricao
names(my_label) <- mapbiomas8_legenda$Class_ID
# Plot one year to check
tm_shape(subset(r85a22n_ag, c(1, 2, 37,38))) +
tm_raster(style = "cat", palette = class_color, labels = my_label, title="") +
tm_scale_bar(breaks = c(0, 10), text.size = 1,
position=c("right", "bottom"))
# Now all years together.
tm_shape(r85a22n_ag) +
tm_raster(style = "cat", palette = class_color, labels = my_label, title="") +
tm_scale_bar(breaks = c(0, 10), text.size = 1,
position=c("right", "bottom")) +
tm_facets(ncol=5) +
tm_layout(legend.bg.color="white", legend.outside = TRUE,
legend.outside.position = "right",
panel.labels = c('1985', '1986', '1987', '1988','1989', '1990',
'1991', '1992', '1993', '1994', '1995','1996',
'1997', '1998',
'1999', '2000', '2001', '2002', '2003',
'2004', '2005', '2006', '2007',
'2008', '2009', '2010', '2011', '2012',
'2013', '2014', '2015',
'2016', '2017', '2018', '2019',
'2020', '2021', '2022')) -> figmap
png("figmap.png", width=6, height=9,
units="in", res = 600)
figmap
invisible(dev.off())