我想在不手动定义范围的情况下裁剪一些栅格的无数据部分( 1中无数据为黑色的图像示例)。
任何的想法?
您可以使用trim
删除仅具有NA
值的外部行和列:
library(raster)
r <- raster(ncols=18,nrows=18)
r[39:49] <- 1
r[205] <- 6
s <- trim(r)
要更改其他值或更改其他值,NA
您可以使用reclassify
. 例如,要更改NA
为 0:
x <- reclassify(r, cbind(NA, 0))
[
为对象定义了子集和[<-
替换方法,raster
因此您可以简单地r[ r[] == 1 ] <- NA
摆脱1
nodata 值所在的值(NAvalue(r)
如果您不确定,请使用找出 R 认为您的 nodata 值应该是什么)。
请注意,您必须r[]
在[
subsetting 命令中使用才能访问这些值。这是一个工作示例...
# Make a raster from system file
logo1 <- raster(system.file("external/rlogo.grd", package="raster"))
# Copy to see difference
logo2 <- logo1
# Set all values in logo2 that are > 230 to be NA
logo2[ logo2[] > 230 ] <- NA
# Observe difference
par( mfrow = c( 1,2 ) )
plot(logo1)
plot(logo2)
我有 2 个略有不同的解决方案。第一个需要手动识别范围,但使用预定义的功能。第二个更自动,但更手工制作。
创建前 2 行为 NA 的可重现栅格
library(raster)
# Create a reproducible example
r1 <- raster(ncol=10, nrow=10)
# The first 2 rows are filled with NAs (no value)
r1[] <- c(rep(NA,20),21:100)
使用 drawExtent() 从绘制的图形中手动获取范围
plot(r1)
r1CropExtent <- drawExtent()
使用从图中选择的范围裁剪栅格
r2 <- crop(r1, r1CropExtent)
作比较
layout(matrix(1:2, nrow=1))
plot(r1)
plot(r2)
它标识栅格中仅具有 NA 值的行和列,并删除位于栅格边缘的行和列。然后它使用 计算范围extent()
。
将栅格转换为一个矩阵,用于标识值是否为 NA。
r1NaM <- is.na(as.matrix(r1))
查找未完全被 NA 填充的列和行
colNotNA <- which(colSums(r1NaM) != nrow(r1))
rowNotNA <- which(rowSums(r1NaM) != ncol(r1))
通过使用未完全由 NA 填充的第一个和最后一个列和行来查找新栅格的范围。用于crop()
裁剪新栅格。
r3Extent <- extent(r1, rowNotNA[1], rowNotNA[length(rowNotNA)],
colNotNA[1], colNotNA[length(colNotNA)])
r3 <- crop(r1, r3Extent)
绘制栅格以进行比较。
layout(matrix(1:2, nrow=1))
plot(r1)
plot(r3)
我已经根据 Marie 的回答编写了一个小函数来快速绘制裁剪的栅格。但是,如果栅格非常大,则可能存在内存问题,因为计算机可能没有足够的 RAM 来将栅格加载为矩阵。
因此,我编写了一个内存安全函数,如果计算机有足够的 RAM(因为它是最快的方法),它将使用 Marie 的方法,或者如果计算机没有足够的 RAM(它速度较慢但内存安全),则使用基于光栅函数的方法)。
这是功能:
plotCroppedRaster <- function(x, na.value = NA)
{
if(!is.na(na.value))
{
x[x == na.value] <- NA
}
if(canProcessInMemory(x, n = 2))
{
x.matrix <- is.na(as.matrix(x))
colNotNA <- which(colSums(x.matrix) != nrow(x))
rowNotNA <- which(rowSums(x.matrix) != ncol(x))
croppedExtent <- extent(x,
r1 = rowNotNA[1],
r2 = rowNotNA[length(rowNotNA)],
c1 = colNotNA[1],
c2 = colNotNA[length(colNotNA)])
plot(crop(x, croppedExtent))
} else
{
xNA <- is.na(x)
colNotNA <- which(colSums(xNA) != nrow(x))
rowNotNA <- which(rowSums(xNA) != ncol(x))
croppedExtent <- extent(x,
r1 = rowNotNA[1],
r2 = rowNotNA[length(rowNotNA)],
c1 = colNotNA[1],
c2 = colNotNA[length(colNotNA)])
plot(crop(x, croppedExtent))
}
}
例子 :
library(raster)
r1 <- raster(ncol=10, nrow=10)
r1[] <- c(rep(NA,20),21:100)
# Uncropped
plot(r1)
# Cropped
plotCroppedRaster(r1)
# If the no-data value is different, for example 0
r2 <- raster(ncol=10, nrow=10)
r2[] <- c(rep(0,20),21:100)
# Uncropped
plot(r2)
# Cropped
plotCroppedRaster(r2, na.value = 0)
如果您使用该rasterVis
软件包(2021 年 6 月 25 日之后的任何版本),它将自动裁剪NA
掉terra
's的值SpatRaster
rasterVis
从 GitHub安装开发版本
if (!require("librarian")) install.packages("librarian")
librarian::shelf(raster, terra, oscarperpinan/rastervis)
# Create a reproducible example
r1 <- raster(ncol = 10, nrow = 10)
# The first 2 rows are filled with NAs (no value)
r1[] <- c(rep(NA, 20), 21:100)
levelplot()
为了r1
rasterVis::levelplot(r1,
margin = list(axis = TRUE))
转换为terra's SpatRaster
然后再次使用levelplot()
r2 <- rast(r1)
rasterVis::levelplot(r2,
margin = list(axis = TRUE))
由reprex 包于 2021-06-26 创建 (v2.0.0 )