从堆叠光栅对象计算选定区域 [英] Calculate selected area from stacked raster object

问题描述

因为我最近才开始使用 R 进行空间分析，所以我无论如何都不是地理学家或空间数据专家，我有一个 - 什么我认为是 - 相对简单的问题.我正在尝试计算满足特定条件的堆叠光栅对象的一部分区域使适应.更具体地说，来自深海的数据集南大西洋，我堆叠了两个光栅对象(深度和坡度)在坐标系 (WGS84) 和 x-y (Lat-Long) 中进一步相同位置.从堆叠的光栅对象，我想提取位于(例如)1000 到 4000 m 深度之间的部分，使用坡度超过10度.我想知道什么地区范围以平方公里为单位，我想将其添加到以前的绘制的地图.下面是一个可重现的示例:

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)

然后，我的问题是:深度(layer.1)和坡度(layer.2)都满足给定条件的面积(在总面积内)是多少?在这种情况下，海拔在 -1000 到 -4000 米之间，坡度大于 10 度.

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.

我最初的想法是:

> area(selectAtt)

给出答案:

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)

哪些是关于 Raster 对象的基本信息...它给了我一种奇怪的感觉，我没有得到我提出的问题的答案.也许我没有问正确的问题?无论如何，它没有告诉我任何符合我条件的区域大小.

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.

然后我做了:

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"’

我试图找出这实际上意味着什么，这似乎是 S3 和 S4 功能之间的不匹配，即使我不完全知道它们是什么.

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 !

推荐答案

有几种方法可以访问图层和栅格堆栈的值.我更喜欢以下内容:

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]][]

现在回到你的问题:

当我想计算满足特定条件的像素面积时，我会执行以下操作(使用小栅格时):

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)

这会为您提供符合定义标准的像素数.如果您在投影坐标系中工作(您在 WGS 84 中工作，因此您无法从中计算出准确的面积)，您只需将 numberOfPixels 乘以您的光栅分辨率:

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

如果您想以平方米为单位获得面积，请将您的栅格重新投影到投影坐标系.例如进入 UTM.在您的情况下，这可能很合适(请注意，您的范围跨越多个区域):http://www.spatialreference.org/ref/epsg/wgs-84-utm-zone-24s/

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

再说一遍:

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 

这篇关于从堆叠光栅对象计算选定区域的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！