library(raster)
library(rgdal)
library(dplyr)
library(sf)
#library(stars)
anos <- list.files('C:/Users/Mario/Desktop/Articulos/Artigo2/Resultados/R/imagens/', pattern = ".tif", full.names = T)
anos
## [1] "C:/Users/Mario/Desktop/Articulos/Artigo2/Resultados/R/imagens/24Sclass_Op_Ra2016.tif"
## [2] "C:/Users/Mario/Desktop/Articulos/Artigo2/Resultados/R/imagens/24Sclass_Op_Ra2017.tif"
## [3] "C:/Users/Mario/Desktop/Articulos/Artigo2/Resultados/R/imagens/24Sclass_Op_Ra2018.tif"
## [4] "C:/Users/Mario/Desktop/Articulos/Artigo2/Resultados/R/imagens/24Sclass_Op_Ra2019.tif"
## [5] "C:/Users/Mario/Desktop/Articulos/Artigo2/Resultados/R/imagens/24Sclass_Op_Ra2020.tif"
stackanos <- stack(anos)
unstack(stackanos)
## [[1]]
## class : RasterLayer
## dimensions : 6310, 2982, 18816420 (nrow, ncol, ncell)
## resolution : 9.832879, 9.832879 (x, y)
## extent : 484332.4, 513654, 8472877, 8534923 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=24 +south +ellps=GRS80 +units=m +no_defs
## source : 24Sclass_Op_Ra2016.tif
## names : X24Sclass_Op_Ra2016
## values : 0, 10 (min, max)
##
##
## [[2]]
## class : RasterLayer
## dimensions : 6310, 2982, 18816420 (nrow, ncol, ncell)
## resolution : 9.832879, 9.832879 (x, y)
## extent : 484332.4, 513654, 8472877, 8534923 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=24 +south +ellps=GRS80 +units=m +no_defs
## source : 24Sclass_Op_Ra2017.tif
## names : X24Sclass_Op_Ra2017
## values : 0, 10 (min, max)
##
##
## [[3]]
## class : RasterLayer
## dimensions : 6310, 2982, 18816420 (nrow, ncol, ncell)
## resolution : 9.832879, 9.832879 (x, y)
## extent : 484332.4, 513654, 8472877, 8534923 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=24 +south +ellps=GRS80 +units=m +no_defs
## source : 24Sclass_Op_Ra2018.tif
## names : X24Sclass_Op_Ra2018
## values : 0, 10 (min, max)
##
##
## [[4]]
## class : RasterLayer
## dimensions : 6310, 2982, 18816420 (nrow, ncol, ncell)
## resolution : 9.832879, 9.832879 (x, y)
## extent : 484332.4, 513654, 8472877, 8534923 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=24 +south +ellps=GRS80 +units=m +no_defs
## source : 24Sclass_Op_Ra2019.tif
## names : X24Sclass_Op_Ra2019
## values : 0, 10 (min, max)
##
##
## [[5]]
## class : RasterLayer
## dimensions : 6310, 2982, 18816420 (nrow, ncol, ncell)
## resolution : 9.832879, 9.832879 (x, y)
## extent : 484332.4, 513654, 8472877, 8534923 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=24 +south +ellps=GRS80 +units=m +no_defs
## source : 24Sclass_Op_Ra2020.tif
## names : X24Sclass_Op_Ra2020
## values : 0, 10 (min, max)
ano16 = stackanos[[1]]
ano20 = stackanos[[5]]
plot(ano16)
# crea um tibble coom os valores do pixel de cada imagen
landcover_tibble = tibble(ano16 = values(ano16), ano20 = values(ano20))
landcover_tibble
## # A tibble: 18,816,420 x 2
## ano16 ano20
## <int> <int>
## 1 9 9
## 2 9 9
## 3 9 9
## 4 9 9
## 5 9 9
## 6 9 9
## 7 9 9
## 8 9 1
## 9 9 9
## 10 9 9
## # ... with 18,816,410 more rows
landcover = filter(landcover_tibble, ano16 != '4' & ano20 != '4')#apagar classe 4 'CDA'
landcover
## # A tibble: 12,515,753 x 2
## ano16 ano20
## <int> <int>
## 1 9 9
## 2 9 9
## 3 9 9
## 4 9 9
## 5 9 9
## 6 9 9
## 7 9 9
## 8 9 1
## 9 9 9
## 10 9 9
## # ... with 12,515,743 more rows
# crea matrix cruzada dos anos agrupando os valores dos pixels
landcover_change = table(landcover)
landcover_change
## ano20
## ano16 0 1 2 3 5 6 7 8 9
## 0 66267 2815 833 11253 131 268 32602 8158 39718
## 1 520 5782873 14 99467 83 67298 9684 7625 453576
## 2 288 2 59417 296 397 0 171 733 7
## 3 4392 340376 1136 346442 314 5704 9868 32379 293970
## 5 74 1076 199 3651 22925 624 9508 16434 5124
## 6 29 163186 44 7349 710 1630755 988 4122 19489
## 7 11194 1204 1708 9282 1237 277 402760 49082 37319
## 8 3473 9950 3082 11031 1358 8630 71504 89348 25263
## 9 14309 162643 186 62202 53 1248 64682 25346 1023648
## 10 1947 47440 710 12831 4397 51325 238725 130214 66092
## ano20
## ano16 10
## 0 479
## 1 2777
## 2 18
## 3 2599
## 5 969
## 6 44259
## 7 25224
## 8 12105
## 9 1104
## 10 215155
# 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*100*100*1/1000000, digits = 2)
landcover_change_matrix
## ano20
## ano16 0 1 2 3 5 6 7 8
## 0 662.67 28.15 8.33 112.53 1.31 2.68 326.02 81.58
## 1 5.20 57828.73 0.14 994.67 0.83 672.98 96.84 76.25
## 2 2.88 0.02 594.17 2.96 3.97 0.00 1.71 7.33
## 3 43.92 3403.76 11.36 3464.42 3.14 57.04 98.68 323.79
## 5 0.74 10.76 1.99 36.51 229.25 6.24 95.08 164.34
## 6 0.29 1631.86 0.44 73.49 7.10 16307.55 9.88 41.22
## 7 111.94 12.04 17.08 92.82 12.37 2.77 4027.60 490.82
## 8 34.73 99.50 30.82 110.31 13.58 86.30 715.04 893.48
## 9 143.09 1626.43 1.86 622.02 0.53 12.48 646.82 253.46
## 10 19.47 474.40 7.10 128.31 43.97 513.25 2387.25 1302.14
## ano20
## ano16 9 10
## 0 397.18 4.79
## 1 4535.76 27.77
## 2 0.07 0.18
## 3 2939.70 25.99
## 5 51.24 9.69
## 6 194.89 442.59
## 7 373.19 252.24
## 8 252.63 121.05
## 9 10236.48 11.04
## 10 660.92 2151.55
rownames(landcover_change_matrix) = c("URB", "FLO", "PED", "AGR", "CAR", "MGE", "SOE", "PHM", "SAV", "AHU")
colnames(landcover_change_matrix) = c("URB", "FLO", "PED", "AGR", "CAR", "MGE", "SOE", "PHM", "SAV", "AHU")
landcover_change_matrix
## ano20
## ano16 URB FLO PED AGR CAR MGE SOE PHM
## URB 662.67 28.15 8.33 112.53 1.31 2.68 326.02 81.58
## FLO 5.20 57828.73 0.14 994.67 0.83 672.98 96.84 76.25
## PED 2.88 0.02 594.17 2.96 3.97 0.00 1.71 7.33
## AGR 43.92 3403.76 11.36 3464.42 3.14 57.04 98.68 323.79
## CAR 0.74 10.76 1.99 36.51 229.25 6.24 95.08 164.34
## MGE 0.29 1631.86 0.44 73.49 7.10 16307.55 9.88 41.22
## SOE 111.94 12.04 17.08 92.82 12.37 2.77 4027.60 490.82
## PHM 34.73 99.50 30.82 110.31 13.58 86.30 715.04 893.48
## SAV 143.09 1626.43 1.86 622.02 0.53 12.48 646.82 253.46
## AHU 19.47 474.40 7.10 128.31 43.97 513.25 2387.25 1302.14
## ano20
## ano16 SAV AHU
## URB 397.18 4.79
## FLO 4535.76 27.77
## PED 0.07 0.18
## AGR 2939.70 25.99
## CAR 51.24 9.69
## MGE 194.89 442.59
## SOE 373.19 252.24
## PHM 252.63 121.05
## SAV 10236.48 11.04
## AHU 660.92 2151.55
# 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(["URB", "FLO", "PED", "AGR", "CAR", "MGE", "SOE", "PHM", "SAV", "AHU"]) .range(["#cd0003", "#299339", "#ffb6f0", "#c9e549", "black", "#1fcc00", "#ce8b43", "#DDA0DD", "#859329", "#518887"])'
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
# you save it as an html
saveNetwork(sn, "sn.html")
saveWidget(sn, file="C:/Mario/Sankey1.html")
saveNetwork(sn, "C:/Mario/sn.html", selfcontained = TRUE)