标签: terra

当我加载 rgdal 包时，我一直注意到这条消息：

“请注意，rgdal 将于 2023 年底退役，计划尽早使用 GDAL 和 PROJ 过渡到 sf/stars/terra 功能。”

我不知道这意味着什么。这是否意味着有一些 sf/stars/terra 函数使用 rgdal 而其他函数则不使用？或者这是否意味着 sf/stars/terra 根本不使用 rgdal，并且我应该使用这些软件包而不是其他依赖 rgdal 的软件包？

哪些依赖 rgdal 的常见软件包会受到影响？

我在完全理解 terra:extract 时遇到问题。我希望提取管理 GADM 多边形的平均栅格值。我的栅格每个国家/地区只有一个值。我希望特定国家/地区内的每个行政多边形都具有相同的值，并且某些包含某些国家/地区边界的多边形会分配面积加权平均值。不幸的是，我当前的脚本并非如此。raster::extract 似乎给出了合理的结果，但不是 terra:extract （请参阅下面的示例代码 - 提供具有不同值的输出）。有人可以根据下面的代码解释一下为什么吗？非常感谢。

## libraries
library(terra)
library(raster)

#===============================================    
## sample example - provides results as expected (1.333, that is (2*0.5+1*1)/1.5)

# sample raster and SpatialPolygons
r <- raster(ncol=2, nrow=3, xmn= 0, ymn= 0, xmx = 30,ymx = 30)
r[] <- c(2, 2, 2, 1, NA, NA)
cds <- rbind(c(7.5,0), c(7.5,20), c(30, 20),c(30,10))
library(sp)
p = Polygon(cds)
ps = Polygons(list(p),1)
sps = SpatialPolygons(list(ps))
plot(r)
plot(sps, add=T)

# test raster package
test1 <- raster::extract(r , …

我正在与 rTerra 合作，并且对来自 LANDFIRE 的 CONUS 历史扰动数据集有疑问： https: //landfire.gov/version_download.php（HDist 是名称）。总结一下我想要做的事情，我想获取这个数据集，裁剪并投影到我的范围内，然后获取单元格的值并将它们分离为图层。因此，我想要一层用于表示严重性，一个用于表示干扰类型，等等。历史干扰数据将这些内容全部包含在一个属性表中。在terra中，这个属性表是在类别下设置的，这带来了很多问题。我没有遇到裁剪或重新投影的问题，它正在进入值并将类别分成图层。我有以下代码

library(terra)


setwd("your pathway to historical disturbance tif here")

h1 <- terra::rast("LC16_HDst_200.tif")  #read in the Hdist tif
h2 <- terra::project(h1, "EPSG:5070", method = "near") #project it using nearest neighbor
h3 <- crop(h2, ext([xmin,xmax,ymin,ymax]) #crop to the extent
h3

然后给出我想要的范围和投影的输出，但主要焦点是类别

categories  : Count, HDIST_ID, DISTCODE_V, DIST_TYPE, TYPE_CONFI, SEVERITY, SEV_CONFID, HDIST_CAT, FDIST, R, G, B

所以我了解到，对于这些类型的数据集，值存储在这些类别下。如果我用 绘图，plot(h3) 我只能得到计数类别的第一行。为了切换该类别我可以使用

activeCat(h3) <- 4
h3

我会得到

name       :   DIST_TYPE
min value  :   Clearcut
max …

基本上就是标题。我知道您可以使用 rast() 读取栅格文件夹，但我只想堆叠两个单独读取的栅格。谢谢

我正在使用rasterpackage 并尝试切换到，terra但由于某些我不明白的原因，无法重现与和等包并行工作时terra的相同操作。这是一个可重现的示例。rastersnowfallfuture.apply

library(terra)
r <- rast()
r[] <- 1:ncell(r)
m <- rast()
m[] <- c(rep(1,ncell(m)/5),rep(2,ncell(m)/5),rep(3,ncell(m)/5),rep(4,ncell(m)/5),rep(5,ncell(m)/5))
ms <- separate(m,other=NA)
plot(ms)
mymask <- function(ind){
  tipo <- tipo_tav[ind]
  mask <- ms[[ind]]
  
  masked <-
    terra::mask(
      r,
      mask
    )
  
  richard <- function(x){
    k <-0.2
    v <-0.3
    a <-200
    y0 <-2
    y <- k/v*x*(1-((x/a)^v))+y0
    return(y)
  }
  pred <- richard(masked)
  pred <- clamp(pred,lower=0)
  return(pred)
}
#the sequential usage works fine, faster than the `raster` counterpart
system.time(x <- mymask(1))#0.03 …

我刚刚注意到terra::cellSize()生成的单元格面积估计值与 生成的值不匹配raster::area()。

首先，为什么这两种方法不能提供相同的答案？第二，哪个估计最准确？请参阅下面的示例。

library(raster)
#> Loading required package: sp
library(terra)
#> terra version 1.3.4

# make test raster with raster::raster()
a <- raster::raster(ncols = 100, nrows = 100,
            xmn = -84, xmx = -83, 
            ymn = 42, ymx = 43)

# make test raster with terra::rast()
b <- terra::rast(ncols = 100, nrows = 100, 
          xmin = -84, xmax = -83, 
          ymin = 42, ymax = 43)

# calculate cell areas (km2)
a_area <- raster::area(a) # km …

我尝试使用“xyz”阅读风格构建一个 2 层的 SpatRast。它使用 3 列作为 rast 函数的输入，但我收到一条包含 4 列的警告消息：

> rast(as.matrix(data.frame(x=c(1,1,2,2),y=c(1,2,1,2),z1=1:4)),type="xyz")
class       : SpatRaster 
dimensions  : 2, 2, 1  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 0.5, 2.5, 0.5, 2.5  (xmin, xmax, ymin, ymax)
coord. ref. :  
data source : memory 
names       : z1 
min values  :  1 
max values  :  4

> r=rast(as.matrix(data.frame(x=c(1,1,2,2),y=c(1,2,1,2),z1=1:4,z2=5:8)),type="xyz")
Warning message:
In v[cells] <- xyz[, 3:d[2]] :
 number of items to replace is not a multiple of replacement length

知道为什么吗？

我有一个关于通过最大重叠进行多边形栅格化的问题，即分配与栅格单元具有最高面积重叠的多边形的值。

现实世界的练习是在 R 中对土壤 ID 的多边形进行栅格化，以便生成相对较低分辨率的土壤属性图作为模型输入。

问题在于rasterize()terra 包（以及类似的 star' st_rasterize()）的函数从包含像元中点的多边形分配像元值。如果栅格单元包含多个多边形，我宁愿选择栅格单元中面积覆盖率最高的多边形（土壤 ID）的值。

这是一个独立的小示例，使用 terra 将我的问题可视化。

library(terra)

f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
r <- rast(v, ncols = 3, nrow = 3)
rcc <- vect(xyFromCell(r, cell = 1:ncell(r)))

x <- rasterize(v, r, field = "NAME_2")
plot(x)
lines(r, col = "light gray")
lines(v)
points(rcc)

大多数情况下，包含单元中心的多边形似乎也具有最高的面积份额。然而，在某些情况下（顶行，第三个单元格），情况并非如此。与多边形相比，单元越大，问题似乎变得越严重。因此，我可以从高分辨率栅格开始，然后使用聚合函数（例如模式）重新采样到所需的（较低）分辨率。但是，也许有人有更有效的想法？

感谢您的帮助！

global()我正在尝试使用terra 包中的函数来计算 spatRaster 中的非 NA 值。"isNA"除了和之外，所有函数（mean、max、sd 等）似乎都可以工作"notNA"。对于这两个函数，它返回此错误：Error in fun(values(x[[i]]), ...) : could not find function "fun"，这与它为拼写错误/不存在的函数返回的错误相同。

r <- rast(ncols=10, nrows=10)
values(r) <- c(1:(ncell(r)-1),NA) # Add one NA value
global(r, fun="mean", na.rm=TRUE) # works
global(r, fun="notNA") # error
global(r, fun="notAfunction") # error

有趣的是，在查看文档 ( ?global) 时，函数描述中对 NA 函数进行了命名，但并未专门列为 的参数选项fun。

那么可以global()计算 NA/非 NA 吗？NA 函数名称是否正确？

编辑：大地版本：1.4.22

我尝试使用 getData （来自栅格包）读取 DEM 栅格，然后将 RasterLayer 转换为 SpatRaster （terra 包）。第一步成功了，但第二步失败了。

library(raster)
library(terra)
 
(alt <- getData('alt', country='DEU', mask=T))
#class      : RasterLayer 
#dimensions : 960, 1116, 1071360  (nrow, ncol, ncell)
#resolution : 0.008333333, 0.008333333  (x, y)
#extent     : 5.8, 15.1, 47.2, 55.2  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : D:/dummy/DEU_msk_alt.grd 
#names      : DEU_msk_alt 
#values     : -179, 2701  (min, max)

altT <- rast(alt)            
# rast is supposed to be able to read RasterLayer, but it fails. Why?
# Error …