r光栅包zonal函数中的值重复
嗨，光栅战士！
经过数月的数据处理，我对结果感到头痛。我对R和空间分析很陌生，但在我的学习过程中很开心。
这里是我的问题：一旦我将zonal函数应用于一组5个光栅对象（.tif），不知何故，我会在我假设为零的区域中得到重复的值。此外，某些值会显示在其他单元格中。
这里我的代码：
folder <- "my_directory"
#for one country
ET<-raster("my_directory/file_ET.tif")
ED<-raster("my_directory/file_ED.tif")
SH<-raster("my_directory/file_SH.tif")
SC<-raster("my_directory/file_SC.tif")
WH<-raster("my_directory/file_WH.tif")
#zones
lowbound<- seq(0,cellStats(ED,"max", na.rm=TRUE),1000)
upbound<-lowbound+1000
band<-upbound/1000
rclmat <- as.matrix(cbind(lowbound,upbound,band)) #to be used on reclassify
rclmat<-rclmat[1:length(rclmat[,1])-1,]
edb<- reclassify(ED, rclmat)
# zonal statistics # I already tried this: na.rm=FALSE, !anyDuplicated(sum_ET)) to avoid duplication – did not work
sum_ET<-zonal(ET,edb,"sum") #, na.rm=FALSE, !anyDuplicated(sum_ET))
stats<-cbind(band,upbound,sum_ET[,2], sum_SH[,2], sum_SC[,2], sum_WH[,2],perc_SH,perc_SC,perc_WH,perc_Tot)
以下是我的结果：
upbound ET
1 1000 10523272.53
2 2000 5156046.27
3 3000 5053895.54
4 4000 4796505.3
5 5000 4392162.97
6 6000 4156065.87
31 31000 10523272.53
32 32000 5156046.27
33 33000 5053895.54
34 34000 4796505.3
35 35000 4392162.97
36 36000 4156065.87
从31到36与从1到6的值相同。
这里是我的一位同事的结果，我正在与他进行比较
upbound ET
1 1000 10523272.53
2 2000 5156046.27
3 3000 5053895.54
4 4000 4796505.3
5 5000 4392162.97
6 6000 4156065.87
31 31000 0
32 32000 21247.26
33 33000 0
34 34000 0
35 35000 0
36 36000 23877.74
正如您所看到的，我得到了重复的值。
输入可以在这里找到ET、ED文件
任何帮助都非常感谢非常感谢
迭戈
共有1个答案
首先，让我用一个简单的例子来解释代码中发生了什么：
# Create "zone" variable
zone=1:10
# Create ET dataframe, that does not have some zones:
dET <- data.frame(zone=c(1:5,7,8), ET=runif(7, min=1, max=10))
dET
# zone ET
# 1 1 5.776128
# 2 2 9.067579
# 3 3 9.874737
# 4 4 7.846662
# 5 5 3.588964
# 6 7 3.843509
# 7 8 8.916714
#merge zone vector and the second column from dET dataframe:
# As you can see some values in the second columnd of the result are repeated to match
# the length of the zone vector and the value that originally belonged to 7th zone was moved into 6th zone
# since there was no 6th zone in the dataframe
cbind(zone, dET[,2])
# zone
# [1,] 1 5.776128
# [2,] 2 9.067579
# [3,] 3 9.874737
# [4,] 4 7.846662
# [5,] 5 3.588964
# [6,] 6 3.843509
# [7,] 7 8.916714
# [8,] 8 5.776128
# [9,] 9 9.067579
# [10,] 10 9.874737
#Instead you should use merge and then fill the missing values with zeros:
result <- merge(data.frame(zone=zone), dET, by="zone", all=TRUE)
result
# zone ET
# 1 1 5.776128
# 2 2 9.067579
# 3 3 9.874737
# 4 4 7.846662
# 5 5 3.588964
# 6 6 NA
# 7 7 3.843509
# 8 8 8.916714
# 9 9 NA
# 10 10 NA
result$ET[is.na(result$ET)] <- 0
# zone ET
# 1 1 5.776128
# 2 2 9.067579
# 3 3 9.874737
# 4 4 7.846662
# 5 5 3.588964
# 6 6 0.000000
# 7 7 3.843509
# 8 8 8.916714
# 9 9 0.000000
# 10 10 0.000000
下面是如何更正代码。在这里，我只使用了您的两个文件ED和ET。所有其他文件都应以相同的方式进行修改：
library(raster)
library(rgdal)
ED<-raster("file_ED.tif")
ET<-raster("file_ET.tif")
lowbound<- seq(0,cellStats(ED,"max", na.rm=TRUE),1000)
upbound<-lowbound+1000
band<-upbound/1000
rclmat <- as.matrix(cbind(lowbound,upbound,band)) #to be used on reclassify
# Commented the following line since cellStats(ED,"max", na.rm=TRUE)=35125.57
# so you should not remove the values in the range 35K-36K
#rclmat<-rclmat[1:length(rclmat[,1])-1,]
ener_dens_band<- reclassify(ED, rclmat)
sum_ET<-zonal(ET,ener_dens_band,"sum")
# Let's look at the tail of the result:
tail(sum_ET)
# zone sum
# [25,] 26 27011.60
# [26,] 27 53905.17
# [27,] 28 18490.15
# [28,] 29 19322.07
# [29,] 32 21247.26
# [30,] 36 23877.74
# As you can see some zones are missing: 30, 31, 33-35
# So they need to be added
sum_ET <- merge(data.frame(zone=1:36), sum_ET, by="zone", all=TRUE)
sum_ET$sum[is.na(sum_ET$sum) ] <- 0
tail(sum_ET, n=10)
# zone sum
# 27 27 53905.17
# 28 28 18490.15
# 29 29 19322.07
# 30 30 0.00
# 31 31 0.00
# 32 32 21247.26
# 33 33 0.00
# 34 34 0.00
# 35 35 0.00
# 36 36 23877.74
您需要对其他列重复相同的操作，然后您可以使用cbind将它们全部放入一个数据帧中。
-
下面的代码在我的图像上生成两个框。我正计划进一步分析这些框内的像素。 在下面的例子中，在红色方块的情况下，我不想继续下去，因为它的右上角有黑色像素。而我想继续在绿色方块的情况下，因为它没有一个黑色像素沿着它的边缘。
-
我正在尝试在R中设置一个randomForest，以便根据其他光栅图像对光栅图像进行分类。我的训练数据是一个完全填充的光栅图像，我想训练许多其他光栅，以尝试基于初始光栅创建光栅输出。代码示例如下： <代码>rf1 ...其中，是我的光栅格式的实际已知值，而到是我想用来预测trainingRaster1是什么的其他光栅图像。我知道您将使用向量或点的训练类来训练一系列光栅，但在我的情况下，我希望使用光
-
我有一个带有空间坐标和一个变量的矩阵数据。空间分辨率为1000米。 我想将其转换为光栅格式。 我使用下面的代码来完成它。但我得到的决心与我得到的不一样。有没有更好的方法可以用我的真实数据获得相同的分辨率？
-
我想知道是否有办法为光栅图层对象创建分区统计数据，特别是R中给定单元值（例如土地使用类别）的计数，而无需重新分类整个光栅。解决方案应具有内存效率，以便处理大型光栅文件，即不需要将值提取到R中的矩阵中。 下面是一个到目前为止我是如何处理它的示例。在这种情况下，我将原始光栅重新分类为仅保留1的利息值和所有其他值的缺失值。 我提出的解决方案创建了冗余数据和额外的处理步骤，以实现我的初始目标。我认为类似于
-
我已经为此挣扎了几个小时。我有一个包含177个多边形（即177个县）的shapefile（称为“shp”）。这个shapefile覆盖在光栅上。我的光栅（称为“ras”）由具有不同污染值的像素组成。 现在我想提取每个多边形的所有像素值及其出现次数。 这正是QGIS功能“分区直方图”所做的。但我想在R中做同样的事情。 我尝试了提取（）函数，并设法获得了每个县的平均值，这已经是第一步，但我想制作像素分
-
我正在尝试以0.0417度的分辨率将一个人工光发射的全球光栅从经度/纬度重新投影到贝尔曼等面积（EPSG:6933）。由于在重投影期间对像素进行插值时，城市地区周围的数据出现峰值，因此整个图层的数据丢失率约为15%。 我尝试将光栅转换为空间点数据框，重新投影空间点数据框，然后使用使用“projectraster”功能创建的光栅作为模板光栅进行光栅化（我认为模板光栅的尺寸、范围和分辨率可能是问题所在