我有一个96个分类栅格的列表(每个栅格都有一个相关变量,即rcp、period和month),其值范围为1-16,我想计算每个栅格中每个类别所覆盖的面积,如果栅格中不存在该类别,则返回NA。
这是我现在创建的函数
a <- function(x){
rs <- x #categorical raster
b <- getValues(area(x, weights=FALSE))
b <- aggregate(b, by=list(getValues(x)), sum, na.rm=T)
b <- as.data.frame(b)
names(b) <- c("CLASS","AREA")
return(b)
}
问题是,它返回的数据帧仅包含现有光栅值,而不包含缺少的值。见下文:
CLASS AREA
1 1 145084.052
2 5 39425.336
3 6 37912.591
4 10 10089.541
5 11 3150.571
6 15 1451.912
7 16 4289.296
如何返回包含所有类别(1-16)的数据帧?并将所有输出合并为一个?列名应为rcp、period和month。
以下是我目前的代码:
mthLs <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
rcpLs <- c("rcp26", "rcp45", "rcp60", "rcp85")
periodLs <- c("2020_2049", "2040_2069")
for(rcp in rcpLs){
rcpDir <- paste0(iDir, "/", rcp)
for(period in periodLs){
for (mth in 1:12) {
rs <- raster(paste0(iDir, "/", rcp, "/", period, "/cng_", mth, ".tif", sep=""))
b <- a(rs) #using function above
}
im <- cbind(RCP=rep(rcp,times=nrow(b)), PERIOD=rep(period,times=nrow(b)), MONTH=rep(mth,times=nrow(b)), b)
}
write.csv(im, paste(oDir, "/output.csv", sep=""), quote=T, row.names=F)
}
示例数据
library(raster)
r <- raster(nc=5, nr=5)
set.seed(20171214)
r1 <- setValues(r, sample(16, ncell(r), replace=TRUE))
r2 <- setValues(r, sample(16, ncell(r), replace=TRUE))
r3 <- setValues(r, sample(16, ncell(r), replace=TRUE))
如果可以制作光栅堆栈(具有相同范围和分辨率的光栅层),则可以执行以下操作:
s <- stack(r1,r2,r3)
x <- freq(s)
m <- data.frame(value=1:16)
for (i in 1:length(x)) {
y <- x[[i]]
colnames(y)[2] <- paste0("v", i)
m <- merge(m, y, all.x=TRUE)
}
m
# value v1 v2 v3
#1 1 3 2 1
#2 2 NA 3 NA
#3 3 2 2 1
#4 4 1 3 2
#5 5 1 NA 2
#6 6 1 3 3
#7 7 NA 3 4
# ...
或者,在没有RasterStack的情况下,创建一个RasterLayer对象列表,例如
# z <- lapply(filenames, raster)
但是对于这个例子,我们可以更手动地执行:
z <- list(r1, r2, r3)
m <- data.frame(value=1:16)
for (i in 1:length(z)) {
y <- freq(z[[i]])
colnames(y)[2] <- paste0("v", i)
m <- merge(m, y, all.x=TRUE)
}
或者——也许是最简单的——使用文件名为ff的向量
m <- data.frame(value=1:16)
for (i in 1:length(ff)) {
y <- freq(raster(ff[i]))
colnames(y)[2] <- paste0("v", i)
m <- merge(m, y, all.x=TRUE)
}
这个可能会有帮助
txt <-
"CLASS,REA
1,145084.052
5,39425.336
6,37912.591
10,10089.541
11,3150.571
15,1451.912
16,4289.296"
b <- read.csv(textConnection(txt),sep = ",")
classes <- 1:16
merge(data.frame(CLASS = classes), b, by='CLASS', all.x=T)
结果:
CLASS REA
1 1 145084.052
2 2 NA
3 3 NA
4 4 NA
5 5 39425.336
6 6 37912.591
7 7 NA
8 8 NA
9 9 NA
10 10 10089.541
11 11 3150.571
12 12 NA
13 13 NA
14 14 NA
15 15 1451.912
16 16 4289.296
我有一个问题,我必须选择部分在一个圆内的所有正方形(想想像素)(即使圆只穿过正方形的一个小角,但如果它穿过角顶点中的一个就不是)。半径是像素大小的整数倍。 问题是圆心在像素之间,即在四个像素的角顶点上。 我只想访问每个像素一次。 例如,我想选择以下图像中的所有白色像素: 对于圆心位于像素中心的圆,这不是问题,我可以使用Bresenham算法的通常形式:
我如何获得土地覆盖类别1、2、3、4、5的总面积,给出以下示例: 一种方法可能是基于这些值对光栅进行子集划分,但是类似于 该网站提供了一个解决方案,但它需要为要分析的每个值创建新的光栅层。我宁愿不这样做,因为我的原始光栅包含许多不同的值,而且非常大。这里介绍了一种类似的方法,但它仅适用于小区域。
我有50多个需要裁剪的光栅文件(ASCII格式)。我已经以ASCII格式从ArcMap导出了遮罩,并将其加载到R中。如何使其适用于一行中的所有光栅,并以与之前相同的名称导出它们(当然是在不同的文件夹中,以避免覆盖)? 我知道光栅软件包中有裁剪功能,但到目前为止我从未使用过。我只是把它们堆放起来做进一步的栖息地分析。 到目前为止,我的代码:
下面的代码在我的图像上生成两个框。我正计划进一步分析这些框内的像素。 在下面的例子中,在红色方块的情况下,我不想继续下去,因为它的右上角有黑色像素。而我想继续在绿色方块的情况下,因为它没有一个黑色像素沿着它的边缘。
我已经为此挣扎了几个小时。我有一个包含177个多边形(即177个县)的shapefile(称为“shp”)。这个shapefile覆盖在光栅上。我的光栅(称为“ras”)由具有不同污染值的像素组成。 现在我想提取每个多边形的所有像素值及其出现次数。 这正是QGIS功能“分区直方图”所做的。但我想在R中做同样的事情。 我尝试了提取()函数,并设法获得了每个县的平均值,这已经是第一步,但我想制作像素分
在计算皮尔逊相关性时,下面的脚本对我来说也适用于相同的数据。我最近对其进行了调整,创建了一个协方差矩阵,以输入到pca中。我在论坛上读到,输入预先创建的协方差矩阵可能会避免记忆问题,但我的情况并非如此。运行协方差矩阵时,我会出现以下错误: 有人能提出一个更有效的方法来做到这一点,这样我就不会遇到内存问题了吗?如果我在计算协方差方面完全偏离了基础,那很好。PCA是我最终唯一需要的东西。我的数据是12