我在这篇文章中使用了解决方案中的代码:
How to pass multiple function in raster::extract()
functions_summ <- function(x, na.rm) c(mean = mean(x, na.rm=na.rm), min = min(x, na.rm=na.rm), max = max(x, na.rm=na.rm), sd = sd(x, na.rm = na.rm))
extraction <- extract( Stack_for_Classification,TrVal, fun=functions_summ, na.rm=TRUE, df = TRUE, sp = T) ## the sp = T. means that I am adding polygon values to the csv
nom <- sapply(TrVal@polygons, slot, "ID")
extraction <- data.frame(ID = nom, Value = extraction)
names(extraction) <- gsub(x = names(extraction), pattern = "\\Value.", replacement = "")
typeof(extraction)
View(extraction)Stack_for_Classification是堆叠栅格,其中包含我想要使用多边形提取的值。TrVal是包含多边形的shapefile。它有一个带有三个标识值的字段"Type“,我用它来根据提取的值来查看三种不同类型之间的差异。问题是,在前面的代码中,我得到了均值、最大值、最小值和后缀为"0,1,2“的sd。总共12。是否可以使用字段“类型”来命名行,例如"mean_Type1","max_Type2"..诸若此类?使用View(提取)后上传图片
发布于 2021-04-23 21:36:23
你可以这样做:
library(raster)
# simulate data -------------------------------------------------------------
r <- s <- t <- raster()
r[] <- 1:ncell(r)
s[] <- sample(1:10,ncell(s),replace = T)
t[] <- runif(ncell(t))
stacked <- stack(r,s,t)
#polys (from the SpatialPolygonsDataFrame example)
Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))
Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))
Sr3 = Polygon(cbind(c(4,4,5,10,4),c(5,3,2,5,5)))
Sr4 = Polygon(cbind(c(5,6,6,5,5),c(4,4,3,3,4)), hole = TRUE)
Srs1 = Polygons(list(Sr1), "s1")
Srs2 = Polygons(list(Sr2), "s2")
Srs3 = Polygons(list(Sr3), "s3")
Srs4 = Polygons(list(Sr4), "s4")
SpP = SpatialPolygons(list(Srs1,Srs2,Srs3,Srs4), 1:4)
plot(SpP, col = 1:4, pbg="white")
types <- letters[1:4]
polys <- SpatialPolygonsDataFrame(SpP,
data = data.frame(
types=types,
row.names = row.names(SpP)
))
#extract
plot(r)
plot(polys,add=T)
functions_summ <-
function(x, na.rm = T){
c(
mean = mean(x, na.rm = na.rm),
min = min(x, na.rm = na.rm),
max = max(x, na.rm = na.rm),
sd = sd(x, na.rm = na.rm)
)
}
ex <- extract(stacked,polys,fun=functions_summ,df=T)
rownames(ex) <-
paste0(gsub("\\..","",rownames(ex)),"_",rep(unique(polys@data$types),each=4))
#convert rownames to column
ex$var <- rownames(ex)但我认为不可能将此data.frame附加到多边形的@data槽中,因为您拥有的观察值比多边形特征更多。您必须将数据调整为宽格式才能合并这两个数据库
发布于 2021-04-24 08:46:01
如果"Type“是多边形的属性,则这些多边形表示为行。使用单层,您可以
library(terra)
v <- shapefile(system.file("external/lux.shp", package="raster"))
r <- raster(v, res=.1)
values(r) <- 1:ncell(r)
f <- function(x, na.rm = T) {
c(mean=mean(x, na.rm = na.rm),
min=min(x, na.rm = na.rm),
max=max(x, na.rm = na.rm),
sd=sd(x, na.rm = na.rm)
)
}
e <- extract(r, v, f)
x <- cbind(v, e)但这不适用于多个层(在RasterStack中);因此,您必须遍历这些层。
使用terra版本1.2-1 (当前为development version),您可以执行以下操作
library(terra)
#terra version 1.2.1
r <- rast(system.file("ex/elev.tif", package="terra"))
x <- c(r, r/2)
names(x) <- c("a", "b")
v <- vect(system.file("ex/lux.shp", package="terra"))
f <- function(x, na.rm = T) {
c(mean=mean(x, na.rm = na.rm),
min=min(x, na.rm = na.rm),
max=max(x, na.rm = na.rm),
sd=sd(x, na.rm = na.rm)
)
}
e <- extract(x, v, f)
head(e)
# ID a.mean a.min a.max a.sd b.mean b.min b.max b.sd
#1 1 467.1052 339 547 34.58480 233.5526 169.5 273.5 17.29240
#2 2 333.8629 195 514 68.03058 166.9315 97.5 257.0 34.01529
#3 3 377.3712 256 517 77.14169 188.6856 128.0 258.5 38.57085
#4 4 373.6000 213 520 82.79129 186.8000 106.5 260.0 41.39564
#5 5 418.6490 293 511 48.35488 209.3245 146.5 255.5 24.17744
#6 6 314.9969 164 403 49.02558 157.4985 82.0 201.5 24.51279https://stackoverflow.com/questions/66205192
复制相似问题