[R-sig-Geo] calculate selected area from stacked raster object

[Loïc Dutrieux]

Thanks for providing an example, see the suggestion inline

On 22/02/2017 10:45, Francesc Montserrat wrote:
> Dear List,
> 
> As I have only recently started using R for spatial analysis, and I am
> not a geographer or spatial data specialist by any means, I have a -what
> I presume to be- relatively simple question. I am trying to calculate
> the area of part of a stacked raster object that meets certain
> conditions. More specifically, from a dataset from the deep sea in the
> south Atlantic, I have stacked two raster objects (depth and slope) that
> are further identical in coordinate system (WGS84), x-y (Lat-Long)
> position and projection. From the stacked raster object, I would like to
> extract the part that sits between (say) 1000 and 4000 m depth, with a
> slope of more than 10 degrees. I would like to know what the areal
> extent is in square km and I would like to add it to a previously
> plotted map. Below is a reproducible example:
> 
> # Raster object containing depth values
> dpt <-  raster(ncol=623, nrow=815, xmx=-31.72083, xmn=-38.50417,
> ymn=-33.8875, ymx=-28.70417)
> values(dpt) <- sample(-200:-5000, size=(nrow(dpt)*ncol(dpt)), replace=T)
> 
> # Raster object containing slope values
> slp <- raster(ncol=623, nrow=815, xmx=-31.72083, xmn=-38.50417,
> ymn=-33.8875, ymx=-28.70417)
> values(slp) <- sample(0:30, size=(nrow(slp)*ncol(slp)), replace=T)
> 
> # Stack raster objects
> stk <- stack(dpt,slp)
> 
# Define function to be passed to calc
fun <- function(x) {
  # consider x to be a vector of length 2
  ifelse(x[1] <= -1000 & x[1] >= -4000 & x[2] >= 10, 1, NA)
}

condition_raster <- calc(x = stk, fun = fun)
area_raster <- area(condition_raster, na.rm = TRUE)
cellStats(area_raster, 'sum')

Cheers,
Loïc

> 
> # Colour palette
> colrs <-
> colorRampPalette(c("navyblue","dodgerblue3","cyan2","green2","darkgoldenrod1"))
> 
> # Plot raster map; does not look like ocean floor because of "sample"
> plot(dpt, xlab="Longitude", ylab="Latitude", col=colrs(100), font.lab=2,
> cex.lab=1.5, las=1)
> 
> 
> # Create a blank copy of previous raster plot
> selectAtt <- raster(dpt)
> # Fill in cells where Attribute(s) meet(s) conditions
> selectAtt[stk$layer.1 <= -1000 & stk$layer.1 >= -4000 & stk$layer.2 >=
> 10] <- 90
> # Set object projection
> projection(selectAtt) <- CRS("+proj=longlat +ellps=WGS84")
> # Plot selection in previous raster
> plot(selectAtt, col="red", add=TRUE, legend=F, proj4string=crswgs84)
> 
> #####>>> What is the area (within total area) with both depth (layer.1)
> and slope (layer.2) meeting given conditions ?
> 
> a <- stk[stk$layer.1 <= -1000 & stk$layer.1 >= -4000 & stk$layer.2 >= 10]
> area(a, na.rm=T)
> 
> # This gives me the error message:
> Error in (function (classes, fdef, mtable)  :
>             unable to find an inherited method for function ‘area’ for
> signature ‘"matrix"’
> 
> I have tried to find what this actually means, but I cannot find a way
> to solve this. What am I missing here ? Any help is greatly appreciated !
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo