rrasterterra

R calculate area within contourlines from raster


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


Solution

  • 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