标签: terra

我有一个脚本，我在其中读取存储为 .tif 的栅格：

f_treecover   <- raster('Landcover_data/treecover_res_100_q.tif')       

该脚本在几个月前运行良好，但现在我收到以下错误消息：

Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer",  : 
  Cannot create a RasterLayer object from this file.

我可以使用 terra 中的“rast”函数读取该文件：

f_treecover   <- rast('Landcover_data/treecover_res_100_q.tif')      
> f_treecover
class       : SpatRaster 
dimensions  : 1400, 3600, 1  (nrow, ncol, nlyr)
resolution  : 0.1, 0.1  (x, y)
extent      : -180, 180, -60, 80  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
source      : treecover_res_100_q.tif 
name        : treecover_res_100_q 

但是，脚本中的所有其他对象都是 RasterLayers ，因此我也需要该对象也是其中之一。我尝试将生成的 SpatRaster 强制转换为 Raster

raster(f_treecover)

，但这会导致与最初尝试读取 …

我想使用包aggregate中的函数terra R来聚合栅格，并使用分位数方法作为聚合函数。下面，我使用quantilefrom 函数R base使用本地包目录中的栅格来计算第 50 个百分位（即中位数）。我选择第 50 个百分位数与中位数进行比较，但我的目标确实是计算其他分位数......

library(terra)
# load elevation coming with the terra pakage
r <- rast( system.file("ex/elev.tif", package="terra") )

plot(r)

# number of iteration
n_it <- 20

# with a custom function
start_time <- Sys.time()
for (i in 1:n_it){
ra <- aggregate(r, 2 , fun = function(x) quantile(x, probs = .5, na.rm = T))
}
end_time <- Sys.time()

我的电脑大约花了 6秒做20次。

print(end_time-start_time)

时差 6.052727 秒

aggregate当我使用中值内置函数运行相同的运行时，大约花费了 执行同样的 20 …

我想使用 terra 将一些栅格堆叠在列表中。

这过去在光栅中非常容易工作，因为它可以列表stack()。然而，这在陆地上已经不可能了。请参阅下面的示例：

r1 <- raster(nrows = 1, ncols = 1, res = 0.5, xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5, vals = 0.1)
r2 <- raster(nrows = 1, ncols = 1, res = 0.5, xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5, vals = 0.2)
r_list <- list(r1, r2)
r_stack <- stack(r_list)

栅格堆栈中的结果：

class      : RasterStack 
dimensions : 6, 6, 36, 2  (nrow, ncol, ncell, nlayers)
resolution …

我正在使用 \xe2\x80\x9cterra\xe2\x80\x9d 包处理 SpatRaster 列表。我已经从脚本创建了栅格列表，并保存了 R 环境。但是，当我从另一个脚本加载 SpatRasters 列表时，遇到以下消息错误：

Error: external pointer is not valid\n

这是一个可重现的示例：

library(terra)\nx <- terra::rast(xmin=-110, xmax=-80, ymin=40, ymax=70, ncols=30, nrows=30)\nvalues(x) <- 1:ncell(x)\nr <- c(x, x, x, x)\nr <- list(r, r, r, r)\nsave(r, file = "test.Rdata")\n\nrm(list=ls(all=TRUE))\n\nload("test.Rdata")\nr\n#[[1]]\n#class       : SpatRaster \n#Error: external pointer is not valid\n

您能否提供解决此问题的指导？\xc2\xa0任何帮助将不胜感激。

我使用 terra 包从光栅文件中提取农药应用。在提取步骤中，terra 似乎删除了一些行。在此示例中，有 19 个警告，内容如下：

无法计算最小/最大，采样中未找到有效像素。（GDAL错误1）

在提取步骤之后，我剩下的行数减少了，并且似乎已经删除了一些行。

此步骤使用栅格效果很好（从这里使用栅格计算农药的平均施用量和总施用量，但数字没有相加），但不幸的是，考虑到我需要处理的文件量，速度很慢。

任何想法这个错误意味着什么以及如何解决它？

这是代码：

## Terra ----

data(wrld_simpl)

## Need to create a SparVector for terra
wrld_simpl = vect(wrld_simpl)

r <- terra::rast("https://raw.github.com/hansronald/Pesticide-data/master/APR_Soybean_Glyphosate_2015_L.tif")

## Remove 
r <- terra::clamp(r, lower=0, values=FALSE)

# area is in ha (values in raster are kg / ha per year)
a <- terra::area(r, sum=FALSE, mask=TRUE) * 0.0001

## Get the total area that pesticide has been applied to
tot_area <- terra::extract(a, wrld_simpl, fun = sum, na.rm = TRUE) …

我正在尝试在 terra 包中按行号和列号对栅格进行子集化。显然，这在栅格中很容易，至少没有地理范围和 crs： 使用行/列索引对栅格进行子集化。但我无法让它在陆地上工作。必须有一个简单的方法。terra::subset 仅选择栅格的图层。

期待有人问为什么：在对高程栅格进行采样并计算坡度和坡向之前，我用行和列填充了栅格，这依赖于相邻像元。现在我需要去掉那些填充的行和列。

library(terra)
EXT <- c( -108, -105, 39, 42 )
R <- rast( extent=EXT, ncol=14, nrow=14, crs="epsg:4326" )
R[] <- 1:ncell(R)

# Now try to strip off the outer 2 rows and columns
crop( x=R, y=ext( 3, 12, 3, 12 ) )
# Error: [crop] extents do not overlap

# Normal R-style subsetting also does not work,
# just gives values of that subset
R[ 3:12, 3:12 ]

st_buffer我对from packagesf和bufferfrom terrapackage on R的用法感到困惑。

我正在围绕点创建缓冲区，并且在使用时sf::st_buffer我知道我需要在参数中使用缓冲区的半径dist，但在使用时terra::buffer我需要在参数中指定缓冲区的直径width。

它是否正确？

我\xc2\xb4ve在这里检查了sf::st_buffer，参数dist： https: //rdrr.io/r/stats/dist.html

在这里terra::buffer，参数width：https ://rdrr.io/cran/terra/man/width.html

我想使用 terra 包进行多边形分析中的一个点。我有一组点，我想提取它们落在哪些多边形中。下面的示例数据显示了单个多边形

library(terra)
crdref <- "+proj=longlat +datum=WGS84"
  
longitude <- c(-116.7, -120.4, -116.7, -113.5, -115.5, -120.8, -119.5, -113.7, -113.7, -110.7)
latitude <- c(45.3, 42.6, 38.9, 42.1, 35.7, 38.9, 36.2, 39, 41.6, 36.9)
lonlat <- cbind(longitude, latitude)
pts <- vect(lonlat, crs=crdref)
  
lon <- c(-116.8, -114.2, -112.9, -111.9, -114.2, -115.4, -117.7)
lat <- c(41.3, 42.9, 42.4, 39.8, 37.6, 38.3, 37.6)
lonlat <- cbind(id=1, part=1, lon, lat)
pols <- vect(lonlat, type="polygons", crs=crdref)
  
plot(pols, border='blue', col='yellow', lwd=3)
points(pts, col='red', pch=20, cex=3)

terra::extract(pts, pols)
       id.y …

我想将SpatVector对象转换为数据框以在 ggplot2 中使用。

pkgs <- c("geodata", "raster", "ggplot2", "tidy")
lapply(pkgs, require, character.only = TRUE)

boundary_GB <- geodata::gadm(country = "GB", path = tempdir(), resolution = 2, level = 1)

我目前的方法需要很长时间：

boundary_GB_df <- broom::tidy(methods::as(boundary_GB, "Spatial"))

剧情：

ggplot(data = boundary_GB_df, mapping = aes(x = long, y = lat, group = group)) +
  geom_polygon(fill = NA, colour = "black")

我对 SpatVector 对象没有经验，有更快的方法吗？

我知道 tidyterra 包（即 tidyterra::geom_spatvector()）。

谢谢

问题

我有空间线（作为一个terra::vect()对象），我想与栅格（terra::rast()）相交，以便通过每条线相交的 id 来识别单元格。

作为输出，我想要一个data.frame带有列line_id和的长表单cell_id（参见下面的示例）。

可能的解决方案

我更喜欢一个terra解决方案，但我也喜欢使用sf、stars、raster等包的解决方案。

对栅格进行矢量化是一种选择，但我的现实生活数据非常大 - 我可能很容易遇到性能问题。所以我想摆脱栅格的矢量化。

我想到了 terra::intersect() , but it does not take SpatRasters as an argument. And terra::rasterizeGeom()` 不返回 ids，但计算长度等......

例子

library("terra")

# Create example lines
v <- vect(system.file("ex/lux.shp", package="terra"))
lns <- as.lines(v)[c(1,7), "ID_1"] # filter for two lines only & only keep ID
names(lns) <- "line_id"
lns
# class       : SpatVector 
# geometry    : lines …