segfault when masking raster with cropped SpatVector

[plantarum]

When cropping a raster with a SpatVector that has itself been cropped I'm getting a repeatable segfault.

MWE:

library(terra)
library(geodata)

## note the path needs to exist on your machine! ##
noam <- gadm(country = c("USA", "CAN"),
             resolution = 2, path = "~/data/gadm/")

mainland <- crop(noam, ext(-180, -30, 0, 70))

wc <- worldclim_global(path = "~/data", res = 10, var = "bio")

wcml <- crop(wc, mainland)
wcml <- mask(wc, mainland)

Crashes with the following output:

 *** caught segfault ***
address (nil), cause 'unknown'

Traceback:
 1: .External(list(name = "CppMethod__invoke_notvoid", address = <pointer: 0x556ffd9bb460>,     dll = list(name = "Rcpp", path = "/home/smithty/R/x86_64-pc-linux-gnu-library/4.4/Rcpp/libs/Rcpp.so",         dynamicLookup = TRUE, handle = <pointer: 0x556fff0828c0>,         info = <pointer: 0x556ffd783770>, forceSymbols = FALSE),     numParameters = -1L), <pointer: 0x557001c46030>, <pointer: 0x5570015a4d30>,     .pointer, ...)
 2: x@cpp$mask_vector(mask@cpp, inverse[1], updatevalue[1], touches[1],     opt)
 3: .local(x, mask, ...)
 4: mask(wc, mainland)
 5: mask(wc, mainland)

I have used is.valid to check mainland is a valid SpatVector.

I've tried different cropping. Setting y max to 84 or higher doesn't produce an error. The upper extent of noam is 83.1104, so maybe this has something to do with how the cropped polygons are handled?

System info:

Ubuntu 20.04.3 LTS Focal Fossa - Release amd64 (20210819)]
R version 4.4.1 (2024-06-14) -- "Race for Your Life"
terra: 1.7.71
terra::gdal(lib="all")

gdal	proj	geos
"3.0.4"	"6.3.1"	"3.8.0"

Thanks for your attention!

[namanpaul]

I have faced a similar error when working with SpatRaster re-projection.
Furthermore, when testing SpatExtent with is.valid(), I got unexpected results.

> rast_extent
SpatExtent : -204.084105293173, -13.5442262185562, 85.8658537594388, 85.9408192387974 (xmin, xmax, ymin, ymax)

> is.valid(rast_extent)
[1] TRUE

[Rapsodia86]

I do get: Error: basic_string::_M_create

wcml <- crop(wc, mainland)
wcml
class       : SpatRaster 
dimensions  : 307, 759, 19  (nrow, ncol, nlyr)
resolution  : 0.1666667, 0.1666667  (x, y)
extent      : -179.1667, -52.66667, 18.83333, 70  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source(s)   : memory
names       : wc2.1~bio_1, wc2.1~bio_2, wc2.1~bio_3, wc2.1~bio_4, wc2.1~bio_5, wc2.1~bio_6, ... 
min values  :   -17.89260,     1.00000,    9.131122,       0.000,      -1.722,    -39.3305, ... 
max values  :    27.14929,    21.14754,  100.000000,    1739.376,      44.687,     22.4000, ... 

wcml <- crop(wc, mainland,mask=TRUE)
Error: basic_string::_M_create

wcml <- mask(wc, mainland)
Error: basic_string::_M_create

sessionInfo()
R version 4.4.0 (2024-04-24 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)
Matrix products: default

locale:
[1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] geodata_0.5-9 terra_1.7-81 

loaded via a namespace (and not attached):
[1] compiler_4.4.0    parallelly_1.37.1 parallel_4.4.0    tools_4.4.0      
[5] Rcpp_1.0.12       codetools_0.2-20 

[plantarum]

@Rapsodia86 I just reran my example on my current version of R and got the same error as initially reported (see below).

Perhaps there is something different in how R on Windows detects and catches errors compared to Linux?

> wcml <- mask(wc, mainland)

 *** caught segfault ***
address (nil), cause 'unknown'

Traceback:
 1: .External(list(name = "CppMethod__invoke_notvoid", address = <pointer: 0x559f2c9dd130>,     dll = list(name = "Rcpp", path = "/home/smithty/R/x86_64-pc-linux-gnu-library/4.4/Rcpp/libs/Rcpp.so",         dynamicLookup = TRUE, handle = <pointer: 0x559f2cb7d140>,         info = <pointer: 0x559f2afa19b0>, forceSymbols = FALSE),     numParameters = -1L), <pointer: 0x559f2b43f650>, <pointer: 0x559f2ac1c140>,     .pointer, ...)
 2: x@ptr$mask_vector(mask@ptr, inverse[1], updatevalue[1], touches[1],     opt)
 3: .local(x, mask, ...)
 4: mask(wc, mainland)
 5: mask(wc, mainland)

My full sessioninfo():

R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 20.04.6 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_CA.UTF-8       LC_NUMERIC=C               LC_TIME=en_CA.UTF-8       
 [4] LC_COLLATE=en_CA.UTF-8     LC_MONETARY=en_CA.UTF-8    LC_MESSAGES=en_CA.UTF-8   
 [7] LC_PAPER=en_CA.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Toronto
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] geodata_0.6-2 terra_1.7-78 

loaded via a namespace (and not attached):
[1] compiler_4.4.1   Rcpp_1.0.13      codetools_0.2-20

[Rapsodia86]

I think the problem is with the shapefile of Canada, which after cropping ends up with 8812 geometries, while originally has 13. Using bbox should not increase the number of geometries. But it does, maybe due to islands? See further down.

The original vector "noam" has 64 geometries (polygons: 51 USA + 13 CAN), and using it, mask() works:
wcml <- mask(wc, noam) , while the cropped vector "mainland" has 8863 (8812 CAN + 51 USA) geometries.

noam_us <- subset(noam, noam$GID_0 == "USA") 
noam_can <- subset(noam, noam$GID_0 == "CAN") 
> crop(noam_us, ext(-180, -30, 0, 70))
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 51, 11  (geometries, attributes) ### looks ok
 extent      : -179.1506, -66.9489, 18.9099, 70  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 
 names       :   GID_1 GID_0       COUNTRY  NAME_1 VARNAME_1 NL_NAME_1 TYPE_1 ENGTYPE_1
 type        :   <chr> <chr>         <chr>   <chr>     <chr>     <chr>  <chr>     <chr>
 values      : USA.1_1   USA United States Alabama   AL|Ala.        NA  State     State
               USA.2_1   USA United States  Alaska AK|Alaska        NA  State     State
               USA.3_1   USA United States Arizona  AZ|Ariz.        NA  State     State
  CC_1 HASC_1 ISO_1
 <chr>  <chr> <chr>
    NA  US.AL US-AL
    NA  US.AK US-AK
    NA  US.AZ US-AZ

> crop(noam_can, ext(-180, -30, 0, 70))
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 8812, 11  (geometries, attributes) ## !!!
 extent      : -141.0069, -52.6189, 41.6769, 70  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 
 names       :   GID_1 GID_0 COUNTRY           NAME_1       VARNAME_1 NL_NAME_1   TYPE_1
 type        :   <chr> <chr>   <chr>            <chr>           <chr>     <chr>    <chr>
 values      : CAN.1_1   CAN  Canada          Alberta              NA        NA Province
               CAN.2_1   CAN  Canada British Columbia Colombie brita~        NA Province
               CAN.3_1   CAN  Canada         Manitoba              NA        NA Province
 ENGTYPE_1  CC_1 HASC_1 ISO_1
     <chr> <chr>  <chr> <chr>
  Province    48  CA.AB CA-AB
  Province    59  CA.BC CA-BC
  Province    46  CA.MB CA-MB

##BUT

> length(crop(noam_can, ext(-180, -30, 0, 70))$GID_0)
[1] 13

And you can not select for the "noam_can" beyond 13, although it shows 8812 geometries!

> crop(noam_can, ext(-180, -30, 0, 70))[13]
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 1, 11  (geometries, attributes)
 extent      : -123.8646, -123.8458, 69.3583, 69.3688  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 
 names       :    GID_1 GID_0 COUNTRY NAME_1       VARNAME_1 NL_NAME_1     TYPE_1 ENGTYPE_1
 type        :    <chr> <chr>   <chr>  <chr>           <chr>     <chr>      <chr>     <chr>
 values      : CAN.13_1   CAN  Canada  Yukon Yukon Territor~        NA Territoire Territory
  CC_1 HASC_1 ISO_1
 <chr>  <chr> <chr>
    60  CA.YT CA-YT
> crop(noam_can, ext(-180, -30, 0, 70))[14]
Error: basic_string::_M_create
> crop(noam_can, ext(-180, -30, 0, 70))[15] ## R session aborted

And look here into a number of geometries depending on the ext() of cropping:

> crop(noam_can, ext(-85, -75, 50, 60))
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 3, 11  (geometries, attributes)
 extent      : -85, -75, 50, 60  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 
 names       :    GID_1 GID_0 COUNTRY  NAME_1    VARNAME_1 NL_NAME_1     TYPE_1 ENGTYPE_1
 type        :    <chr> <chr>   <chr>   <chr>        <chr>     <chr>      <chr>     <chr>
 values      :  CAN.8_1   CAN  Canada Nunavut           NA        NA Territoire Territory
                CAN.9_1   CAN  Canada Ontario Upper Canada        NA   Province  Province
               CAN.11_1   CAN  Canada  Québec Lower Canada        NA   Province  Province
  CC_1 HASC_1 ISO_1
 <chr>  <chr> <chr>
    62  CA.NU CA-NU
    35  CA.ON CA-ON
    24  CA.QC    NA
> plot(crop(noam_can, ext(-85, -75, 50, 60)))

> crop(noam_can, ext(-80, -70, 50, 60))
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 2001, 11  (geometries, attributes)
 extent      : -80, -70, 50, 60  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 
 names       :    GID_1 GID_0 COUNTRY  NAME_1    VARNAME_1 NL_NAME_1     TYPE_1 ENGTYPE_1
 type        :    <chr> <chr>   <chr>   <chr>        <chr>     <chr>      <chr>     <chr>
 values      :  CAN.8_1   CAN  Canada Nunavut           NA        NA Territoire Territory
                CAN.9_1   CAN  Canada Ontario Upper Canada        NA   Province  Province
               CAN.11_1   CAN  Canada  Québec Lower Canada        NA   Province  Province
  CC_1 HASC_1 ISO_1
 <chr>  <chr> <chr>
    62  CA.NU CA-NU
    35  CA.ON CA-ON
    24  CA.QC    NA
> plot(crop(noam_can, ext(-80, -70, 50, 60)))

Hi @rhijmans, do you know what is going on with the cropping of Canada gadm shapefile? Thanks!