重采样栅格
Resample raster
我正在尝试将具有高分辨率(25 米)和分类数据(1 到 13)的森林覆盖栅格重新采样为具有较低分辨率(~ 1 公里)的新 RasterLayer
。我的想法是将森林覆盖数据与其他低分辨率栅格数据结合起来:
我尝试了 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
如您所见,新栅格具有所需的分辨率,但也有很多零。我想这是正常的,因为我在 resample
.
中使用了“ngb”选项
第二个策略是使用 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 的系数:"The number of cells is the number of cells of x divided by fact*fact
(when fact is a single number)." 尽管如此,生成的栅格数据似乎没有与我的相同的行数和列数其他低分辨率栅格。
我从来没有使用过这种高分辨率的数据,而且我的计算能力也有限(其中一些命令可以使用 clusterR
并行化,但有时它们与非并行化命令花费的时间相同,特别是因为它们不适用于最近邻计算)。
我缺乏想法;也许我可以尝试 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。这样做可以让您降低分辨率,同时保留尽可能多的信息。 N.B 如果您需要与 %
不同的汇总统计数据,通过更改 aggregate
中的 fun =
函数,应该很容易使以下方法适应您想要的任何统计数据。 =16=]
以下方法非常快(在我的笔记本电脑上只需几秒钟即可处理具有 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
最近邻重采样也不是一个好主意,因为结果是任意的。
这是您建议的逐层方法 。
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)
现在你有比例了。如果你想要一个 class 你可以做
y <- which.max(x)
在这种情况下,另一种方法是聚合 classes
ag <- aggregate(r_hr, fact, modal)
agx <- resample(ag, r_lr, method='ngb')
注意agx
和y
是一样的。但它们都可能有问题,因为您可能有 5 个 classes,每个大约 20%,这使得选择一个获胜者相当不合理。
我正在尝试将具有高分辨率(25 米)和分类数据(1 到 13)的森林覆盖栅格重新采样为具有较低分辨率(~ 1 公里)的新 RasterLayer
。我的想法是将森林覆盖数据与其他低分辨率栅格数据结合起来:
我尝试了
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
如您所见,新栅格具有所需的分辨率,但也有很多零。我想这是正常的,因为我在
resample
. 中使用了“ngb”选项
第二个策略是使用
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 的系数:"The number of cells is the number of cells of x divided byfact*fact
(when fact is a single number)." 尽管如此,生成的栅格数据似乎没有与我的相同的行数和列数其他低分辨率栅格。
我从来没有使用过这种高分辨率的数据,而且我的计算能力也有限(其中一些命令可以使用 clusterR
并行化,但有时它们与非并行化命令花费的时间相同,特别是因为它们不适用于最近邻计算)。
我缺乏想法;也许我可以尝试 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。这样做可以让您降低分辨率,同时保留尽可能多的信息。 N.B 如果您需要与 %
不同的汇总统计数据,通过更改 aggregate
中的 fun =
函数,应该很容易使以下方法适应您想要的任何统计数据。 =16=]
以下方法非常快(在我的笔记本电脑上只需几秒钟即可处理具有 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
最近邻重采样也不是一个好主意,因为结果是任意的。
这是您建议的逐层方法
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)
现在你有比例了。如果你想要一个 class 你可以做
y <- which.max(x)
在这种情况下,另一种方法是聚合 classes
ag <- aggregate(r_hr, fact, modal)
agx <- resample(ag, r_lr, method='ngb')
注意agx
和y
是一样的。但它们都可能有问题,因为您可能有 5 个 classes,每个大约 20%,这使得选择一个获胜者相当不合理。