从堆叠的栅格对象计算选定区域

Question

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

# 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 m 之间，坡度角 >10 度。

我最初的想法是：

> 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 对象的基本信息是什么……它给了我一种奇怪的感觉，即我没有得到我提出的问题的答案。也许我没有问正确的问题？无论如何，它并没有告诉我任何符合我条件的区域的大小。

然后我做了：

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 功能之间的不匹配，即使我不完全知道它们是什么。

无论如何，我以为我在对空间数据提出一个相对简单的查询，即基于栅格堆栈中多个图层的信息的选择对应的区域是多少？我在这里想念什么？任何帮助是极大的赞赏 ！

Answer 1

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

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

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

现在回到你的问题：

当我想计算满足特定标准的像素面积时，我会执行以下操作（使用小栅格时）：

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

这将为您提供满足定义标准的像素数。如果您将在投影坐标系中工作（您在 WGS 84 中工作，因此您无法从中计算准确的区域），您只需将 numberOfPixels 乘以您的栅格分辨率即可：

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

如果您想获得以平方米为单位的面积，请将栅格重新投影到投影坐标系。例如进入 UTM。在您的情况下，这可能很合适（请注意，您的范围跨越多个区域）：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 

Answer 2

好的，这是一种可能的不同方法，基于将栅格单元转换为多边形，然后利用包 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

这与 OP 从 cellStats 方法获得的非常相似（大约 2% 的折扣），但与 reprojection 方法仍然完全不同。根据我在areaPolygon 的帮助下的理解，这应该是一种数字正确的方法。不过，我仍然不明白@maRtin answer 的不同结果。

顺便说一句，对于 OP：你为什么要使用如此“奇怪”的参考系统（纬度/经度，不同的 x/y 分辨率）设置一个示例？

Answer 3

所以首先进行一些健全性检查...

我们在这里处理的实际区域是什么？ 我们将基于您的堆栈，我们将使用res() 来提取分辨率。手动抓取它们很容易，但我更喜欢调用它们。

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

 70.32058

所以我得到 ~70.3 作为堆栈的实际区域。我只是从你提供的数据中得到这个，所以如果 CRS 或预测有问题，那就是你了。如果我们努力为我们的子集找到区域并发现它与此大不相同，我们就误入歧途了。

现在我将从您执行此操作的第 4 个代码块开始：

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

所以现在我们有了一个层，我们在其中对stk 进行了子集化，现在它是一个符合条件的值列表。我们不需要那个列表，我们只是想知道它包含多少元素。

替换为：

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

因为我们只想知道列表中有多少元素...

现在，使用res() 获取 x 和 y 分辨率：

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

结果：

 29.80777

根据我们之前计算的栅格面积，这似乎是合理的，约占总面积的 42%。这在我看来是合理的。

关于其他答案 cmets 中的投影和坐标系问题，是的这是一个问题。但是你要确保你的预测不会有问题，并且对你想要做的事情有意义。还要记住，简单地重新分配投影系统与重新投影不同。我看到在 R 中使用栅格时最常见的错误之一可能是简单地重新分配 CRS 而不是重新投影栅格。您需要知道您的数据是以什么形式进入的，以及您需要用它做什么才能对它进行您想要做的那种数学运算。