时间戳

首页 > 解决方案 > 如何更好地编写用于在 R 中提取栅格数据的代码。也许是一个循环？

r - 如何更好地编写用于在 R 中提取栅格数据的代码。也许是一个循环？

问题描述

我仍然对 R 感到困惑。我有一个想法，可以使用以下脚本将一些 ArcView 函数迁移到 R 中。我希望能够让我的 df 帮助填写。我是否可以通过为底部关注的每个指标（范围、avg1、avg2）编写长管道然后通过循环运行它来解决这个问题？制作我的哑数据没有意义，因为我没有可用的哑栅格数据。

library(rgdal)
library(sp)
library(raster)

df <- data.frame(location = paste('location', 1:4), 
             latitude = c(35.123+(1:4)), longitude = c(-120.123+(1:4)))

#steps 创建 shapefile

#加载栅格数据

a <- raster(x= 'data/raster/a.tif')#datum=NAD83
b <- raster(x= 'data/raster/b.tif')#datum=NAD83
c <- raster(x= 'data/raster/c.tif')#datum=NAD83

#提取自

location1_ws_a <- raster::extract(a,location, df= TRUE)
location1_ws_b <- raster::extract(b,location, df= TRUE)
location1_ws_c <- raster::extract(c,location, df= TRUE)

range <- max(location1_ws_b)-min(location1_ws_b)
avg1 <- mean(location1$a)
avg2 <- mean(location1$c)

output <- data.frame(location,
                 range,
                 avg1,
                 avg2)

最后，我想要一个包含列名称 c(location, range, avg1, avg2) 的表，其中包含行 location1、location2、location3 及其结果非常感谢！！！

标签： rr-rasterrgdal

解决方案

以下是如何在示例中包含栅格数据

library(raster)
df <- data.frame(location = paste('location', 1:4), 
             latitude = c(35.123+(1:4)), longitude = c(-120.123+(1:4)))
             
r <- raster(ymn=36, ymx=40, xmn=-120, xmx=-116, nrow=100, ncol=100)
set.seed(0)
a <- setValues(r, sample(ncell(r)))
b <- setValues(r, sample(ncell(r)))
c <- setValues(r, sample(ncell(r)))
s <- stack(a, b, d)
names(s) <- c("a", "b", "c")

使用提取物。请注意，顺序是经度/纬度，并且在您的示例数据中这些是相反的。因此使用df[,3:2]来纠正这个问题。

e <- extract(s, df[,3:2])

并绑定

out <- cbind(df, e)
out

#    location latitude longitude    a    b    c
#1 location 1   36.123  -119.123 3276 3575 5866
#2 location 2   37.123  -118.123 6785 4712 2944
#3 location 3   38.123  -117.123 9586 7439 7887
#4 location 4   39.123  -116.123 5874 9826 3570

从这里开始，我不太确定你想要什么。也许下面？

out$avg1 <- mean(out$a)
out$avg2 <- mean(out$c)
out$range <- max(out$b) - min(out$b)

推荐阅读