尝试使用带有Terra的SpatVector对SpatRaster进行掩膜时,出现“错误:[蒙版]无法创建数据集”

vpfxa7rd  于 2023-01-06  发布在  其他
关注(0)|答案(2)|浏览(161)

简单地尝试使用矢量(.shp)来屏蔽使用terra::mask的SpatRaster;得到以下错误

>Error: \[mask\] cannot create dataset

LCC84 <- rast("C:/Users_forest_VLCE2_1984.tif") 
vec <- vect("C:/Users/Land_Management_Units.shp")
vec_proj <- project(vec, LCC84)
LCC84_masked <- terra::mask(LCC84, vec_proj)

错误:[mask]无法创建数据集

vec
#class       : SpatVector 
#geometry    : polygons 
#dimensions  : 1, 8  (geometries, attributes)
#extent      : -117.3165, -115.1691, 50.70613, 52.27127  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat NAD83 (EPSG:4269) 

LCC84 
#class      : SpatRaster 
#dimensions : 128340, 193936, 1 (nrow, ncol, nlyr) 
#resolution : 30, 30 (x, y) 
#extent     : -2660911, 3157169, -851351.9, 2998848 (xmin, xmax, ymin, ymax) 
#coord. ref.: Lambert_Conformal_Conic_2SP 
#source     : CA_forest_VLCE2_1984.tif 
#name       : CA_forest_VLCE2_1984

crs(LCC84, proj=TRUE)
[1] "+proj=lcc +lat_0=49 +lon_0=-95 +lat_1=49 +lat_2=77 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"
czq61nw1

czq61nw11#

您可以使用以下代码

library(terra)
library(sf)

#Read the data
LCC84 <- rast("C:/Users_forest_VLCE2_1984.tif")
vec <- st_read("C:/Users/Land_Management_Units.shp")

#Convert the crs of shapefile
vec_proj <- sf::st_transform(vec, crs(LCC84))

#Masking the raster using the shapefile
LCC84_masked <- terra::mask(LCC84, vec_proj)
xxls0lw8

xxls0lw82#

你提供的数据对我有用

library(terra)
#terra 1.6.51

v <- vect("Extent_BNP_Extact.shp")
r <- rast("CA_forest_VLCE2_1984.tif")
pv <- project(v, r)
z <- crop(r, pv, mask=T)

r
#class       : SpatRaster 
#dimensions  : 128340, 193936, 1  (nrow, ncol, nlyr)
#resolution  : 30, 30  (x, y)
#extent      : -2660911, 3157169, -851351.9, 2998848  (xmin, xmax, ymin, ymax)
#coord. ref. : Lambert_Conformal_Conic_2SP 
#source      : CA_forest_VLCE2_1984.tif 
#name        : CA_forest_VLCE2_1984 

v
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 1, 2  (geometries, attributes)
# extent      : 474065.5, 635666.6, 5613645, 5798288  (xmin, xmax, ymin, ymax)
# source      : Extent_BNP_Extact.shp
# coord. ref. : NAD83 / UTM zone 11N (EPSG:26911) 
# names       : Shape_Leng Shape_Area
# type        :      <num>      <num>
# values      :  6.744e+05  2.828e+10

plot(z)

第一次裁剪看起来很合理,因为整个数据集有250亿个细胞。

mask(r, pv)

运行它花了一段时间,但它工作了。如果它不适合你,我猜你可能没有足够的磁盘空间在你的临时文件夹。查看terraOptions()的临时文件夹的位置。
同时,您还可以执行以下操作

pv <- project(v, "EPSG:4269")

但这毫无意义,因为您的光栅数据没有EPSG:4269坐标参考系(lon/lat NAD 83)。

相关问题