在 R 中折叠栅格数据答案

【问题标题】：Collapsing raster data in R在 R 中折叠栅格数据
【发布时间】：2021-04-02 18:10:25
【问题描述】：

我目前正在处理来自https://www.worldclim.org/data/worldclim21.html 的气候数据，并且正在尝试获取秘鲁 196 个省份的最高和最低温度。到目前为止，我已经使用了这段代码：

my_shape <- st_read("multipolygon shapfile with 196 provinces inside Peru") 
r <- raster("raster file of temp with 1km2 grid")
r2 <- crop(r, extent(my_shape))

到目前为止，这给了我秘鲁的形状。由此，我想获得 196 个省份中每个省份的平均（最高和最低）温度，并最终得到一个包含 196 个观测值的数据框。有没有有效的方法？

【问题讨论】：

嗨阿方索！请提供一个可重现的示例：以便试图帮助您的人可以复制粘贴您的代码，运行它并能够重现您的步骤。在您的解释中附上图片也很有帮助，让您更清楚地了解您的问题

标签： r raster shapefile temperature

【解决方案1】：

您可以为此使用zonal。类似的东西

library(raster)
my_shape <- shapefile("multipolygon shapfile with 196 provinces inside Peru") 
r <- raster("raster file of temp with 1km2 grid")
r2 <- crop(r, my_shape)
zones <- rasterize(my_shape, r2)
zonal(r2, zones)

【讨论】：

【解决方案2】：

您可以使用 raster::extract() 根据您的多边形汇总栅格数据。
我经常发现有必要指定要提取的包以确保您获取正确的版本。

library(sp) # I'm assuming that's the source for st_read
library(raster)  

my_shape <- st_read("multipolygon shapfile with 196 provinces inside Peru") 
r <- raster("raster file of temp with 1km2 grid")
## not necessary to crop source raster

my_shape$mean <- raster::extract(r, my_shape, fun = mean, na.rm = TRUE)

这应该会在您的省份多边形图层中添加一个属性列。你可以用 fun = min 和 fun = max 做同样的事情。

如果您的栅格图层很大，或者您想多次执行此操作，则使用 velox 等替代包可能会获得更好的性能。

library(velox)

vs <- velox("raster file of temp with 1km2 grid")
my_shape$mean2 <- vs$extract(my_shape, fun = function(x){mean(x, na.rm= TRUE)})

velox 方法要快得多（见下文）

另一种选择是使用 stars 包，并使用aggregate() 来汇总最小和最大数据。见How to extract values from a raster using polygons with the R stars package?

library(stars)
library(microbenchmark)

rstar = read_stars(raster_file)
st_crs(rstar) <- st_crs(my_shape) 

f1 = function(){
      my_shape$mean = raster::extract(r, my_shape, fun = mean, na.rm = TRUE)
}
f2 = function(){
      vs <- velox(raster_file)
      my_shape$mean2 <- vs$extract(my_shape, fun = function(x){mean(x, na.rm= TRUE)})
}

f3 = function(){
      zones = rasterize(my_shape, r)
      raster::zonal(r, zones)
}

f4 = function(){
x = aggregate(rstar, my_shape, mean, na.rm = T)
st_as_sf(x)
}

res = microbenchmark(f1(), f2(), f3(), f4(), times = 3)
res
# Unit: seconds
# expr        min         lq       mean     median       uq        max neval
# f1() 172.283024 174.644299 178.920625 177.005575 182.2394 187.473276     3
# f2()   3.339221   3.374647   3.396874   3.410072   3.4257   3.441327     3
# f3() 170.981220 179.481681 182.520730 187.982142 188.2905 188.598829     3
# f4() 167.076081 173.427232 179.332300 179.778383 185.4604 191.142437     3

velox$extract 显然是性能赢家。其他三个版本看起来非常相似。

zonal 的处理时间主要用于对多边形图层进行栅格化，如果您用于多个指标，它会超过提取和聚合。从文件中读取栅格包含在 velox 表达式中，但不包含在其他版本中。

【讨论】：