I am rather new to GIS and I don't know how to perform the following operation: [ R programming language] I am developing a project evaluating forest cover change over 10 years (2010-2020) for a specific area in the state of São Paulo/ Brazil. I have this area as a shapefile, and I have managed to clip this area from the raster of forest cover change in a single year (I can sure do it for the other years). I then proceeded to extract the value of the resulting raster to inspect the output. However, I don't know where to go from there to compile the forest cover change data over the years and compute statistics. The output from the extraction is a list I cannot interpret. Does anyone know how to do this? So far, I have used the sf, sp and raster packages.
1 Answer
After extracting you just need to get the area of each class, which is resolution_x*resolution_y (or resolution^2 if the same) by the number of pixels of that class.
In the following code data is grouped by year and class, after the extract part you already did.
library(raster)
library(tidyr)
library(ggplot2)
rs = stack() # build a random stack just to use for the example
for(i in 1:8) {
m = matrix(sample(1:8, 100, TRUE), nrow = 10)
r = raster(m, xmn = 700000, xmx = 700100, ymn = 2100000, ymx = 2100100, crs = "+proj=utm +zone=14 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
r = as.factor(r)
names(r) = paste0("year", i)
rs = stack(rs, r)
}
sample_pol = st_as_sfc("POLYGON ((700013 2100005, 700050 2100005, 700050 2100070, 700013 2100070, 700013 2100005))")
sample_pol = as(sample_pol, "Spatial")
rs_extract = raster::extract(rs, sample_pol) # you already did this part
as.data.frame(rs_extract) %>%
gather() %>% # change the shape of the data frame
group_by(key, value) %>% summarise(cover_area = n()*100) %>% # 10*10 m resolution
ggplot(aes(key, cover_area, fill = as.factor(value))) + geom_col() +
scale_fill_discrete(name = "Class")
-
Thank you for the reply! I don't understand all of it though: where do I retrieve the values in this line of code: r = raster(m, xmn = 700000, xmx = 700100, ymn = 2100000, ymx = 2100100, crs = "+proj=utm +zone=14 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")– Julia R.Mar 27, 2021 at 8:16
-
That's just a sample raster I built for demonstrating the code, I wrote them by hand from XY values Mar 27, 2021 at 14:00
-
1Kudos, that is a nice graph, lots of information in a small space. Feb 11, 2022 at 16:08