0

I want to calculate the area under every possible combination of two classified rasters. I am using the following code

library(terra)

#First create two rasters
r1 <- r2 <- rast(nrow=100, ncol=100)

#Assign random cell values
set.seed(123)
values(r1) <- runif(ncell(r1), min=0, max=1)
values(r2) <- runif(ncell(r2), min=0, max=1)

# classify the values into two groups
m_r1 <- c(min(global(r1, "min", na.rm=TRUE)), 0.2, 1,
               0.2, max(global(r1, "max", na.rm=TRUE)), 2)

m_r2 <- c(min(global(r2, "min", na.rm=TRUE)), 0.2, 1,
              0.2, max(global(r2, "max", na.rm=TRUE)), 2)

#Reclassify the rasters
rclmat_r1 <- matrix(m_r1, ncol=3, byrow=TRUE)
rc_r1 <- classify(r1, rclmat_r1, include.lowest=TRUE)

rclmat_r2 <- matrix(m_r2, ncol=3, byrow=TRUE)
rc_r2 <- classify(r2, rclmat_r2, include.lowest=TRUE)

plot(rc_r1)
plot(rc_r2)

#Convert to polygons
r1_poly <- as.polygons(rc_r1, dissolve=TRUE)
r2_poly <- as.polygons(rc_r2, dissolve=TRUE)
plot(r1_poly)
plot(r2_poly)

#Perform intersections
x <- intersect(r1_poly, r2_poly)
x
#> class       : SpatVector 
#> geometry    : polygons 
#> dimensions  : 2747, 2  (geometries, attributes)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> names       : lyr.1 lyr.1
#> type        : <int> <int>
#> values      :     1     1
#>                   1     2
#>                   2     1

As you can see from the output, one combination i.e. 2-2 is missing. Why is this happening? When I am trying to calculate the area for each combination using expanse(x), it returns a long result. How can I get the area in km2 for the following combinations?

Combination Area (km2)
1-1
1-2
2-1
2-2
UseR10085
  • 7,120
  • 3
  • 24
  • 54
  • I suspect 2-2 missing is a difference between applying `from-to-becomes` vs `cuts`, or suggest it as a line to explore. – Chris Apr 23 '22 at 15:07

1 Answers1

2

With this example it would be better to stay with raster data.

x = 10 * rc_r1 + rc_r2
a = cellSize(x, unit="km")
zonal(a, x, sum)
#  lyr.1      area
#1    11  19886611
#2    12  81946082
#3    21  84763905
#4    22 323469024

By multiplying with 10, the values in the first layer become 10 (if they were 1) or 20 (if they were 2). If you then add the second layer, you get 10 + 1 or 2 and 20 + 1 or 2, so you end up with four classes: 11, 12, 21, and 22. These show the value in the first raster (first digit) and in the second raster (second digit).

When you show a SpatVector only the first three records are printed, and there is a 2-2 record. Nevertheless, intersect did not work properly and I have now fixed this.

UseR10085
  • 7,120
  • 3
  • 24
  • 54
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63