##### SHAPE #####
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(rgdal)
## Carregando pacotes exigidos: sp
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-27, (SVN revision 1148)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/Mario/Documents/R/win-library/4.1/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/Mario/Documents/R/win-library/4.1/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was C:/Users/Mario/Documents/R/win-library/4.1/rgdal/proj
# carregar shape BAHIA em SIRGAS 2000, UTM 24S -------------------------------------------------------------------
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1
shpBA <- sf::read_sf("C:/Territorio_do_Baixo_Sul_da_Bahia/Shapes/BA-UTM24S.shp")
plot(shpBA)
#### filtro Linhas da tabela, neste caso por os nomes dos municipios
shp <- shpBA %>%
filter(NM_MUNICIP %in% c("IBICOARA"))
plot(shp)
library(raster)
##
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
##
## select
anos <- list.files("C:/Users/Mario/Desktop/Articulos/IBICOARA/Imagens2/imgR/", pattern = ".tif", full.names = T)
# Stacking vectors concatenates multiple vectors into a single vector ------------------
stackanos <- stack(anos)
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Sistema de Referencia Geocentrico para las
## AmericaS 2000 in Proj4 definition
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Sistema de Referencia Geocentrico para las
## AmericaS 2000 in Proj4 definition
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Sistema de Referencia Geocentrico para las
## AmericaS 2000 in Proj4 definition
plot(stackanos)
# CROP THEN MASK imagens mapbiomas com shape contorno do baixo sul BA---------------
anos.crop <- crop(x = stackanos,y = shp)
anos.mask <- mask(anos.crop,shp)
##### Unstacking vectors concatenates -----------------------------------
unstack(anos.mask)
## [[1]]
## class : RasterLayer
## dimensions : 1114, 1613, 1796882 (nrow, ncol, ncell)
## resolution : 29.80822, 29.80822 (x, y)
## extent : 223920.3, 272000.9, 8503589, 8536795 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=24 +south +ellps=GRS80 +units=m +no_defs
## source : memory
## names : RE24Smapbiomas.60.ibicoara.ba.1990
## values : 2, 128 (min, max)
##
##
## [[2]]
## class : RasterLayer
## dimensions : 1114, 1613, 1796882 (nrow, ncol, ncell)
## resolution : 29.80822, 29.80822 (x, y)
## extent : 223920.3, 272000.9, 8503589, 8536795 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=24 +south +ellps=GRS80 +units=m +no_defs
## source : memory
## names : RE24Smapbiomas.60.ibicoara.ba.2020
## values : 3, 128 (min, max)
plot(anos.mask)
ano1990 = stackanos[[1]]
ano2020 = stackanos[[2]]
# crea um tibble coom os valores do pixel de cada imagen
landcover_tibble = tibble(ano1990 = values(ano1990), ano2020 = values(ano2020))
landcover_tibble
## # A tibble: 1,825,740 x 2
## ano1990 ano2020
## <int> <int>
## 1 128 128
## 2 128 128
## 3 128 128
## 4 128 128
## 5 128 128
## 6 128 128
## 7 128 128
## 8 128 128
## 9 128 128
## 10 128 128
## # ... with 1,825,730 more rows
landcover = filter(landcover_tibble, ano1990 != "128" & ano2020 != "128")#apagar classe 4 "CDA"
landcover
## # A tibble: 923,635 x 2
## ano1990 ano2020
## <int> <int>
## 1 2 33
## 2 2 33
## 3 2 33
## 4 4 33
## 5 4 33
## 6 4 33
## 7 4 33
## 8 4 33
## 9 4 33
## 10 4 33
## # ... with 923,625 more rows
# crea matrix cruzada dos anos agrupando os valores dos pixels
landcover_change = table(landcover)
landcover_change
## ano2020
## ano1990 3 4 9 12 24 29 33
## 2 260998 5124 16579 111 18 28 365
## 4 3636 185814 48144 1399 580 83 770
## 9 8294 52967 138913 22769 3877 389 95
## 12 379 14400 57749 93407 592 448 52
## 24 11 311 876 813 2332 91 5
## 29 12 112 318 228 5 116 21
## 33 49 63 6 11 0 21 254
# realiza la conversion de pixl para km, multiplica por tamanho de pixel e divide por 1 km
#### landcover_change_matrix = round(addmargins(landcover_change)*100*100*1/1000000, digits = 2)
landcover_change_matrix = round(landcover_change*300*300*1/1000000, digits = 2)
landcover_change_matrix
## ano2020
## ano1990 3 4 9 12 24 29 33
## 2 23489.82 461.16 1492.11 9.99 1.62 2.52 32.85
## 4 327.24 16723.26 4332.96 125.91 52.20 7.47 69.30
## 9 746.46 4767.03 12502.17 2049.21 348.93 35.01 8.55
## 12 34.11 1296.00 5197.41 8406.63 53.28 40.32 4.68
## 24 0.99 27.99 78.84 73.17 209.88 8.19 0.45
## 29 1.08 10.08 28.62 20.52 0.45 10.44 1.89
## 33 4.41 5.67 0.54 0.99 0.00 1.89 22.86
rownames(landcover_change_matrix) = c("FLO1 1990", "FLO2 1990", "AGR 1990", "FLO3 1990", "URB 1990", "ROC 1990", "CDA 1990")
colnames(landcover_change_matrix) = c("FLO1 2020", "FLO2 2020", "AGR 2020", "FLO3 2020", "URB 2020", "ROC 2020", "CDA 2020")
landcover_change_matrix
## ano2020
## ano1990 FLO1 2020 FLO2 2020 AGR 2020 FLO3 2020 URB 2020 ROC 2020 CDA 2020
## FLO1 1990 23489.82 461.16 1492.11 9.99 1.62 2.52 32.85
## FLO2 1990 327.24 16723.26 4332.96 125.91 52.20 7.47 69.30
## AGR 1990 746.46 4767.03 12502.17 2049.21 348.93 35.01 8.55
## FLO3 1990 34.11 1296.00 5197.41 8406.63 53.28 40.32 4.68
## URB 1990 0.99 27.99 78.84 73.17 209.88 8.19 0.45
## ROC 1990 1.08 10.08 28.62 20.52 0.45 10.44 1.89
## CDA 1990 4.41 5.67 0.54 0.99 0.00 1.89 22.86
# Sankey
# Libraries
library(tidyverse)
#library(viridis)
library(patchwork)
library(hrbrthemes)
#library(circlize)
# Package
data = as.data.frame.matrix(landcover_change_matrix)
# I need a long format
data_long <- data %>%
rownames_to_column %>%
gather(key = 'key', value = 'value', -rowname) %>%
filter(value > 0)
colnames(data_long) <- c("source", "target", "value")
data_long$target <- paste(data_long$target, " ", sep="")
# From these flows we need to create a node data frame: it lists every entities involved in the flow
nodes <- data.frame(name=c(as.character(data_long$source), as.character(data_long$target)) %>% unique())
# With networkD3, connection must be provided using id, not using real name like in the links dataframe.. So we need to reformat it.
data_long$IDsource = match(data_long$source, nodes$name)-1
data_long$IDtarget = match(data_long$target, nodes$name)-1
# prepare color scale: I give one specific color for each node.
colores <- 'd3.scaleOrdinal() .domain(["FLO1", "FLO2", "AGR", "FLO3", "URB", "ROC", "CDA"]) .range(["129912", "00ff00", "f1c232", "B8AF4F", "red", "gray", "0000FF"])'
library(htmlwidgets)
library(networkD3)
# Make the Network
sn = sankeyNetwork(Links = data_long, Nodes = nodes,
Source = 'IDsource', Target = 'IDtarget', Value = 'value',
NodeID = 'name',
sinksRight = FALSE,
nodeWidth = 40,
fontSize = 14,
fontFamily = "Ariel", units = "%",
nodePadding = 20,
colourScale = colores)
sn