r语言 - 在光栅砖/堆栈的图层上绘制的更好方法



我正在尝试绘制特定点的光栅砖中的所有值。这是为了为特定像素的遥感数据创建光谱图。

我可以通过多种方式做到这一点,但它们非常笨重且缓慢(见下面的示例(。这很慢,主要是因为将大型栅格文件转换为矩阵会占用大量内存。

有没有更好的方法来使用baseR或整洁的诗句来做到这一点,或者在其中一个光栅/遥感包中内置方法来做到这一点?

下面是一个可重现的示例:

library (raster)
library (rgdal)
library (sp)
library (tidyr)
library (ggplot2)
library (dplyr)

##############################
### dummy data
##############################

coord1 <- c(50, 47, 45)
coord2 <- c(50, 51, 49)
frame <- data.frame(coord1, coord2)
coords <- coordinates(frame)
x = list(1, 2, 3, 4, 5)
y = list(1, 2, 3, 4, 5)
for (i in 1:length(x)) { # create clusters based on coords
set.seed(i+100)
x[[i]] <- rnorm(5000, mean = coords[, 1], sd = 1)
y[[i]] <- rnorm(5000, mean = coords[, 2], sd = 1)
}
obs <- data.frame(x, y)
names(obs) <- c('x1', 'x2', 'x3', 'x4', 'x5', 'y1', 'y2', 'y3', 'y4', 'y5')
coordinates(obs) <- c('x1', 'y1') # this creates a spatial points data frame
# create blank raster of study area
ex <- extent(c(45, 50, 48, 52))
target <- raster(ex, resolution = 0.5)
# create raster brick
r_list <- list()
for (i in 1:ncol(obs)) {
r_list[[i]] <- rasterize(obs, target, fun = 'count')
}
obs_frequency <- brick(r_list)

这是一个可能但缓慢的解决方案

############################
### Example Solution
############################
vals_all <- obs_frequency[, ] ### this is very slow ###
vals_all <- data.frame(vals_all)
### sample values
points <- sample_n(vals_all, size = 5)
points <- t(points)
points <- data.frame(points)
points_tidy <- gather(points)
points_tidy$xval <- rep(1:8, times = 5)

### plot
ggplot(points_tidy, aes(x = xval, y = value)) + geom_line(aes(colour = key)) + theme_classic()

我使用raster::extract函数找到了更好的解决方案。这直接对值进行采样,避免将整个光栅砖变成破坏内存的矩阵。

值得注意的是,在这里,使用砖块比使用光栅堆栈快得多。

############################
### Extract values and plot 
############################
### extract values
points <- c(49, 50, 51) #arbitrary points
pointvals <- raster::extract(obs_frequency, points) ##### USE THE RASTER::EXTRACT FUNCTION
### manipulate data structure
pointvals <- data.frame(t(pointvals))
points_tidy <- gather(pointvals)
points_tidy$xval <- rep(1:8, times = 3)
### plot
ggplot(points_tidy, aes(x = xval, y = value)) + geom_line(aes(colour = key)) + theme_classic()

最新更新