segfault when masking raster with cropped SpatVector
plantarum opened this issue · 6 comments
plantarum commented
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 commented
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 commented
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 commented
@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 commented
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!
rhijmans commented
I believe this now works as expected.
A work around would have been to do
wcml <- mask(wc, noam)