Calculate selected area from stacked raster object

问题:

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) and x-y (Lat-Long) position. 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)


 # 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)

Then, my question would be: what is the area (within total area) with both depth (layer.1) and slope (layer.2) meeting given conditions ? In this case, with elevation between -1000 and -4000 m and with a slope angle of >10 degrees.

My initial thought was to do:

> area(selectAtt)

giving the answer:

class       : RasterLayer 
dimensions  : 623, 815, 507745  (nrow, ncol, ncell)
resolution  : 0.008323108, 0.008319957  (x, y)
extent      : -38.50417, -31.72083, -33.8875, -28.70417  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 0.7083593, 0.7481782  (min, max)

Which is the basic information about the Raster object...it gave me the strange sensation that I was not getting the answer to the question I posed. Maybe I was not asking the correct question ? In any case, it did not tell me anything about the size of the area meeting my conditions.

Then I did:

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, and it appears to be a mismatch between S3 and S4 functionalities, even though I don't exactly know what those are.

Anyhow, I thought I was posing a relatively simple query to the spatial data, namely what is the area corresponding to a selection based on information from several layers from a raster stack ? What am I missing here ? Any help is greatly appreciated !

回答1:

There are several methods how to access the layers and the values of a raster stack. I prefer the following:

#select first layer of raster stack 
stk[[1]]

#get values of second layer
stk[[2]][]

Now back to your question:

When I want to calculate the area of pixels that meet certain criteria, I do the following (when working with small rasters):

numberOfPixels <- sum(stk[[1]][] <= -1000 & stk[[1]][] >= -4000 & stk[[2]][] >= 10, na.rm=T)

This gives you the number of pixels that meet the defined criteria. If you would be working in a projected coordinate system (you are working in WGS 84 and therefore you can't calculate accurate areas from this) you would simply multiply the numberOfPixels by the resolution of your raster:

area <- numberOfPixels * (res(stk)[1] * res(stk)[2])

If you want to get the area in squared meters, reproject your raster to a projected coordinate system. For example into UTM. In your case this one might be a good fit (note that your extent is stretching across multiple zones): http://www.spatialreference.org/ref/epsg/wgs-84-utm-zone-24s/

stk <- projectRaster(from=stk, crs="+proj=utm +zone=24 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

Then again:

numberOfPixels <- sum(stk[[1]][] <= -1000 & stk[[1]][] >= -4000 & stk[[2]][] >= 10, na.rm=T)
area <- numberOfPixels * (res(stk)[1] * res(stk)[2])
area
[1] 258898897920 

回答2:

Ok, here is a possible different approach, based on converting the raster cells to polygons and then computing the area of the polygons exploiting package geosphere:

require(geosphere)
polys = rasterToPolygons(stk)
polys_sub = subset(polys, (layer.1 <= -1000 & layer.1 >= -4000 & layer.2 >= 10))
> sum(areaPolygon(polys_sub))/1E6  #in km2

[1] 157035.2

This is quite similar with what the OP got from the cellStats method (about 2% off), and still completely different from the reprojection method. From what I understand from the help of areaPolygon this should be a numerically correct approach. I still don't understand the different results with respect to @maRtin answer, though.

By the way, for the OP: why are you setting-up an example with a so "strange" reference system (lat/long, different x/y resolution) ?

回答3:

So a couple of sanity checks first...

What is the actual area we're dealing with here? We'll base that on your stack and we'll use res() to pull the resolutions. It would be easy enough to grab them manually, but I prefer to call them.

 length(stk)*res(stk)[1]*res(stk)[2]

 70.32058

So I'm getting ~70.3 as the actual area of the stack. I'm only getting this from what you provided for data, so if the CRS or projections are wonky, thats on you. If we do the work to find the area for our subset and find it to be wildly different than this, we went astray.

Now I'm going to start at your 4th block of code where you did this:

 a <- stk[stk$layer.1 <= -1000 & stk$layer.1 >= -4000 & stk$layer.2 >= 10]

So now we have a layer, a where we've subsetted stk such that it is now a list of values which meet the criteria. We don't need that list, we really just want to know how many elements it contains.

Replace with:

 a<-length(stk[stk$layer.1 <= -1000 & stk$layer.1 >= -4000 & stk$layer.2 >= 10])

since we only want to know how many elements are in the list...

now, using res() to grab the x and y resolutions:

 a*res(stk)[1]*res(stk)[2]

Results in:

 29.80777

Based on the area of the raster we calculated previously, this seems reasonable, ~42% of the total area. This seems reasonable to me.

Regarding the issues with projections and coordinate systems in the other answers comments, YES that's a problem. But its on you to make sure your projections aren't wonky and make sense for what you are trying to do. Also remember that simply reassigning a projection system is not the same as re-projection. Probably one of the most common mistakes I see working with rasters in R is simply reassigning a CRS rather than reprojecting a raster. Its on you to know what form your data is coming in as and what you need to do with it to do the kind of math you want to do on it.