0

Simply trying to use a vector (.shp) to mask a SpatRaster using terra::mask; get the following error

>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)

Error: [mask] cannot create dataset

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"
user217532
  • 1
  • 1
  • Welcome to SO. Can you provide these two datasets through dropbox or Google Drive? Please visit [How to make a great R reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) – UseR10085 Jan 05 '23 at 04:45
  • You bet Shapefile (https://drive.google.com/drive/folders/1pHsohVKwl5K0d5wk8U_U0_O0NRAnEpa8?usp=share_link) LCC84 is being uploaded to google drive, but estimated 25hrs (on WIFI) available directly here... (https://opendata.nfis.org/downloads/forest_change/CA_forest_VLCE2_1984.zip) – user217532 Jan 05 '23 at 22:15

2 Answers2

0

You can use the following code

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)
UseR10085
  • 7,120
  • 3
  • 24
  • 54
0

It works for me with the data you provided

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)

enter image description here

First cropping seems logical here because the entire dataset has ~25 billion cells. I also tried

mask(r, pv)

It took a while to run, but it worked. If it does not for you, my guess would be that you may not have sufficient disk-space in your temp folder. See terraOptions() for the location of the temp folder.

Also, you do the equivalent of

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

But that makes no sense, as your raster data do not have the EPSG:4269 coordinate reference system (lon/lat NAD83).

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • I do get the same error when I try rasterize step alone. > show(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 – user217532 Jan 05 '23 at 22:09
  • I should also mention, terra::crop is working just fine; only mask throws the error. – user217532 Jan 05 '23 at 22:13
  • The comments are not a good place for this. The idea is that you update your question,. Can you also add `crs(LCC84, proj=TRUE)` ? – Robert Hijmans Jan 05 '23 at 23:30
  • Thanks for letting me know that; I appreciate it. I updated/cleaned up the question. I could get it to work with the crop/Mask=T code you had above - so that's great! Not sure if I need to keep fiddling to get the mask function to work for me... – user217532 Jan 06 '23 at 02:44
  • `mask` worked for me, see my expanded answer – Robert Hijmans Jan 06 '23 at 04:51