top*_*cat 5 gis r resolution raster
我正在尝试将具有高分辨率(25米)和分类数据(1到13)的森林覆盖栅格重新采样为RasterLayer具有较低分辨率(~1 km)的新光栅.我的想法是将森林覆盖数据与其他低分辨率栅格数据结合起来:
我试过raster::resample(),但由于数据是绝对的,我丢失了很多信息:
summary(as.factor(df$loss_year_mosaic_30m))
0 1 2 3 4 5 6 7 8 9 10 11 12 13
3777691 65 101 50 151 145 159 295 291 134 102 126 104 91
如您所见,新栅格具有所需的分辨率,但也有很多零.我认为这是正常的,因为我使用了'ngb'选项resample.
第二种策略是使用,raster::aggregate()但我发现难以定义因子整数,因为分辨率的变化不是直截了当的(如分辨率的两倍或相似).
我的高分辨率栅格具有以下分辨率,我希望它将其聚合0.008333333, 0.008333333 (x, y)到相同程度的分辨率.
loss_year
class : RasterLayer
dimensions : 70503, 59566, 4199581698 (nrow, ncol, ncell)
resolution : 0.00025, 0.00025 (x, y)
extent : -81.73875, -66.84725, -4.2285, 13.39725 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : /Volumes/LaCie/Deforestacion/Hansen/loss_year_mosaic_30m.tif
names : loss_year_mosaic_30m
values : 0, 13 (min, max)
在aggregate帮助描述之后,我尝试了一个~33.33的因子:"单元格的数量是x除以的单元格数fact*fact(当事实是单个数字时)." 尽管如此,生成的栅格数据似乎与我的其他低分辨率栅格具有相同的行数和列数.
我从来没有使用过这种高分辨率的数据,而且我的计算也是有限的(这些命令中的一些可以使用并行化clusterR,但有时它们比非并行化命令花费的时间相同,特别是因为它们不能用于最近的neighboor计算).
我缺乏想法; 也许我可以尝试layerize获得一个计数栅格,但我必须'聚集'并且factor问题出现了.由于这个过程花了我几天的时间来处理,我确实想知道创建较低分辨率栅格的最有效方法,而不会丢失太多信息
可重现的例子如下:
r_hr <- raster(nrow=70, ncol=70) #High resolution raster with categorical data
set.seed(0)
r_hr[] <- round(runif(1:ncell(r_hr), 1, 5))
r_lr <- raster(nrow=6, ncol=6) #Low resolution raster
第一个策略:信息丢失
r <- resample(r_hr, r_lr, method = "ngb") #The raster data is categorical
第二种策略:难以定义总因子
r <- aggregate(r_hr, factor) #How to define a factor to get exactly the same number of cells of h_lr?
另外一个选项: layerize
r_brick <- layerize(r_hr)
aggregate(r_brick, factor) #How to define factor to coincide with the r_lr dimensions?
谢谢你的帮助!
将土地覆盖图聚合为 %cover 层是非常标准的做法。也就是说，你的目标应该是生成 13 层，每一层都类似于该网格单元中的 %cover。这样做可以降低分辨率，同时保留尽可能多的信息。NB如果你需要比不同的汇总统计%，应该很容易适应以下方法给您想要的任何统计数据，通过改变fun =函数aggregate。
以下方法非常快（在我的笔记本电脑上处理具有 1 亿个像元的栅格只需几秒钟）：
首先，让我们创建一些虚拟栅格来使用
Nhr <- 1e4 # resolution of high-res raster
Nlr <- 333 # resolution of low-res raster
r.hr <- raster(ncols=Nhr, nrows=Nhr)
r.lr <- raster(ncols=Nlr, nrows=Nlr)
r.hr[] <- sample(1:13, Nhr^2, replace=T)
现在，我们首先将高分辨率栅格聚合到与低分辨率几乎相同的分辨率（到最接近的整数单元格）。每个结果图层都包含该像元内原始栅格值为 N 的区域部分。
Nratio <- as.integer(Nhr/Nlr) # ratio of high to low resolutions, to nearest integer value for aggregation
layer1 <- aggregate(r.hr, Nratio, fun=function(x, na.rm=T) {mean(x==1, na.rm=na.rm)})
layer2 <- aggregate(r.hr, Nratio, fun=function(x, na.rm=T) {mean(x==2, na.rm=na.rm)})
最后，将低分辨率栅格重新采样到所需的分辨率
layer1 <- resample(layer1, r.lr, method = "ngb")
layer2 <- resample(layer2, r.lr, method = "ngb")
对每一层重复，并将您的层构建到堆栈或多波段栅格中
r_hr <- raster(nrow=70, ncol=70) #High resolution raster with categorical data
set.seed(0)
r_hr[] <- round(runif(1:ncell(r_hr), 1, 5))
r_lr <- raster(nrow=6, ncol=6)
r_hr
#class : RasterLayer
#dimensions : 70, 70, 4900 (nrow, ncol, ncell)
#resolution : 5.142857, 2.571429 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
#data source : in memory
#names : layer
#values : 1, 5 (min, max)
r_lr
#class : RasterLayer
#dimensions : 6, 6, 36 (nrow, ncol, ncell)
#resolution : 60, 30 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
直接聚合是不可能的,因为70/6不是整数.
dim(r_hr)[1:2] / dim(r_lr)[1:2]
#[1] 11.66667 11.66667
最近邻重新采样也不是一个好主意,因为结果是任意的.
这是你建议的逐层方法,dww也已经显示了.
b <- layerize(r_hr)
fact <- round(dim(r_hr)[1:2] / dim(r_lr)[1:2])
a <- aggregate(b, fact)
x <- resample(a, r_lr)
现在你有比例.如果你想要一门课,你可以做
y <- which.max(x)
在这种情况下,另一种方法是聚合类
ag <- aggregate(r_hr, fact, modal)
agx <- resample(ag, r_lr, method='ngb')
请注意agx并且y是相同的.但是它们都可能有问题,因为你可能有5个类的单元格,每个类别大约有20%,因此选择一个胜利者是不合理的.