My objective is to calculate the area within a closed isoline from a raster.
I know the topic has been already covered in this forum (link) but I'm here wondering if there is a less convoluted way to proceed compared to what I've managed so far (which I admit it's pretty hawful and I'm not so sure it's even correct, indeed!).
Given this reproducible example
library(terra)
r <- rast(system.file("ex/elev.tif", package="terra"))
r <- project(r, "EPSG:3035")
e <- ext(4020500, 4029000, 2975000, 2982000)
rc <- crop(r, e)
plot(rc)
contour(rc, levels = 500, add = TRUE)
rc_df <- flip(rc) |> as.data.frame(xy = TRUE)
# check plot
contour(x = sort(unique(rc_df$x)),
y = sort(unique(rc_df$y)),
z = matrix(rc_df$elevation, nrow = 11, ncol = 9, byrow = F),
levels = 500
)
cl_500 <- contourLines(x = sort(unique(rc_df$x)),
y = sort(unique(rc_df$y)),
z = matrix(rc_df$elevation, nrow = 11, ncol = 9, byrow = F),
levels = 500
)
library(splancs)
sum(sapply(cl_500, function(cl){areapl(cbind(cl$x, cl$y))}))
I would like to know if someone could point me to a different direction for accomplishing the task. It seems to me very weak the part of the code transforming the raster to a dataframe which is then fed into the contourLines function, after some preliminary manipulations (the part I've been struggling a lot!).
thanks
I would use terra::as.countour()
and then a bit of massaging in {sf} - we have to convert (MULTI)LINESTRINGs to POLYGONs:
p <- terra::as.contour(rc, levels = 500)
plot(rc)
plot(p, add = TRUE)
p |>
sf::st_as_sf() |>
sf::st_cast(to = "LINESTRING") |>
sf::st_cast(to = "POLYGON") |>
dplyr::summarise() |>
sf::st_area()
#> 4011447 [m^2]
Or just in {terra}:
p <- terra::as.contour(rc, levels = 500)
p |>
terra::as.polygons() |>
terra::expanse() |>
sum()
#> Warning: [as.polygons] dropped attributes
#> [1] 4011447
Created on 2024-12-09 with reprex v2.1.1.9000