Error: [[[<-,SpatVector] cannot add these values after erase() with specific geometries in a specific order

[erkent-carb]

I originally asked the question on stackoverflow here and in my testing since think it is possibly related to erase(). Thanks for your helpful advice on figuring out where the issue happens.

I have a SpatVect (this_weight) which I then erase parts from using two other SpatVects (I will call these masks 1 and 2). After the second erase() call, I am unable to add a new column (Error: [[[<-,SpatVector] cannot add these values), plot it, or use writeVector() to save it (Error: basic_string::_M_create). In the course of testing this I found that if I only used one mask, or if I reversed the order of the masking (i.e., erase mask 2 first and then mask 1) then I did not get these errors. It seems it is possibly an issue with my specific geometries, though they appear valid. The reprex below shows first the original ordering test (errors) and the reverse order test (no errors). The attached zip file contains the three necessary shapefiles used.

Any help figuring out if it is just these geometries or a bug are appreciated.

bug_testing_2025-01-16.zip

library(terra)
library(ggplot2)

io_dir <- "C:/Users/ekent/Desktop/bug_testing_2025-01-16/"

#################################################################################################################################################
# I get the error when running this code

# the SpatVect I am working with (a subset of another file obtained using query() with the filter argument)
this_weight <- vect(file.path(io_dir,"this_weight_original.shp"))
nrow(weight)
this_weight$test <- 1:nrow(this_weight)

# first mask
mask1 <- terra::vect(file.path(io_dir,"this_protected_lands_mask.shp"))
# now project the mask features to the crs of the weight file
mask1_projected <- terra::project(mask1, this_weight)
mask1_projected_valid <-  terra::makeValid(mask1_projected)

# erase the parts of the weight file features that fall inside the first mask
this_weight <- terra::erase(this_weight, mask1_projected_valid)
nrow(this_weight)
# adding new column here okay
this_weight$test2 <- 1:nrow(this_weight)

# geometry seems valid
any(!is.valid(this_weight))

# writing to file succeeds here
writeVector(this_weight, filename = file.path(io_dir,"this_weight_after_protected_lands_applied4.shp"))

# Plotting succeeds here
ggplot()+
  theme_bw() +
  geom_spatvector(data = this_weight)

# second mask
mask2 <- terra::vect(file.path(io_dir,"this_coast_mask2.shp"))
# now project the mask features to the crs of the weight file
mask2_projected <- terra::project(mask2, this_weight)
mask2_projected_valid <-  terra::makeValid(mask2_projected)

# erase the parts of the weight file features that fall inside the second mask
this_weight <- terra::erase(this_weight, mask2_projected_valid) # this line appears to be where things go wrong, though we don't get an error here
nrow(this_weight)

# geometry seems valid
any(!is.valid(this_weight))

# when trying to write at this point I get:
# Error: basic_string::_M_create
writeVector(this_weight, filename = file.path(io_dir,"this_weight_after_protected_lands_and_coast_applied4.shp"))

# I am also not able to plot after the last erase()
# Error in `$<-.data.frame`(`*tmp*`, geometry, value = c("010600000.........
ggplot()+
  theme_bw() +
  geom_spatvector(data = this_weight)

# this is the line when we have the original error: Error: [[[<-,SpatVector] cannot add these values
this_weight$test3 <- 1:nrow(this_weight)


#################################################################################################################################################
# if I reverse the order the masks are applied, I don't get the error anymore

# the SpatVect I am working with (a subset of another file obtained using query() with the filter argument)
this_weight <- vect(file.path(io_dir,"this_weight_original.shp"))
nrow(weight)
this_weight$test <- 1:nrow(this_weight)

# second mask (this time applied first)
mask2 <- terra::vect(file.path(io_dir,"this_coast_mask2.shp"))
# now project the mask features to the crs of the weight file
mask2_projected <- terra::project(mask2, this_weight)
mask2_projected_valid <-  terra::makeValid(mask2_projected)

# erase the parts of the weight file features that fall inside the second mask
this_weight <- terra::erase(this_weight, mask2_projected_valid) # this line appears to be where things go wrong, though we don't get an error here
nrow(this_weight)

# geometry seems valid
any(!is.valid(this_weight))

# writes fine
writeVector(this_weight, filename = file.path(io_dir,"this_weight_after_protected_lands_and_coast_applied5.shp"))

# plots fine
ggplot()+
  theme_bw() +
  geom_spatvector(data = this_weight)

# this is the line when we had the original error, but not in this order
this_weight$test3 <- 1:nrow(this_weight)

# first mask
mask1 <- terra::vect(file.path(io_dir,"this_protected_lands_mask.shp"))
# now project the mask features to the crs of the weight file
mask1_projected <- terra::project(mask1, this_weight)
mask1_projected_valid <-  terra::makeValid(mask1_projected)

# erase the parts of the weight file features that fall inside the first mask
this_weight <- terra::erase(this_weight, mask1_projected_valid)
nrow(this_weight)
# adding new column here okay
this_weight$test2 <- 1:nrow(this_weight)

# geometry seems valid
any(!is.valid(this_weight))

# writing to file succeeds here
writeVector(this_weight, filename = file.path(io_dir,"this_weight_after_protected_lands_applied6.shp"))

# Plotting succeeds here
ggplot()+
  theme_bw() +
  geom_spatvector(data = this_weight)

[rhijmans]

Thank you. I believe this has now been fixed. At least your example now works:

library(terra)
w1 <- vect("this_weight_original.shp")
m1 <- terra::vect("this_protected_lands_mask.shp")
m2 <- terra::vect("this_coast_mask2.shp")

m1v <- terra::project(m1, w1) |> terra::makeValid()
m2v <- terra::project(m2, w1) |> terra::makeValid() |> crop(ext(w1) + 10000)
values(w1) <- data.frame(w1=1:nrow(w1))

w2 <- terra::erase(w1, m1v)
w2$w2 <- 1:nrow(w2)
w3 <- terra::erase(w2, m2v) 
w3$w3 <- 1:nrow(w3)
w3
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 193, 3  (geometries, attributes)
# extent      : -222827.4, -156509.7, 86497.14, 161041  (xmin, xmax, ymin, ymax)
# coord. ref. : MM5_sphere 
# names       :    w1    w2    w3
# type        : <int> <int> <int>
# values      :     1     1     1
#                   2     2     2
#                   3     3     3

[erkent-carb]

Thanks very much. I confirmed that my example works now too.