rspatial / terra

R package for spatial data handling https://rspatial.github.io/terra/reference/terra-package.html
GNU General Public License v3.0
543 stars 90 forks source link

R crash when transforming to polygons #1625

Closed latot closed 4 weeks ago

latot commented 4 weeks ago

Hi, here something weird:

Using this package for raster source: http://github.com/hypertidy/sds

terra::as.polygons(terra::vect(sds::CGAZ(), query = sds::CGAZ_sql("Chile")))
 class       : SpatVector 
 geometry    : polygons 
 dimensions  : 135, 1  (geometries, attributes)
 extent      : -109.4555, -66.42422, -55.91544, -17.49834  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 
/usr/lib/gcc/x86_64-pc-linux-gnu/13/include/g++-v13/bits/stl_vector.h:1128: std::vector<_Tp, _Alloc>::reference std::vector<_Tp, _Alloc>::operator[](size_type) [with _Tp = std::__cxx11::basic_string<char>; _Alloc = std::allocator<std::__cxx11::basic_string<char> >; reference = std::__cxx11::basic_string<char>&; size_type = long unsigned int]: Assertion '__n < this->size()' failed.
Aborted (dumped core)

There is more weird things, I checked this on stackoverflow:

a <- terra::as.polygons(terra::vect(sds::CGAZ(), query = sds::CGAZ_sql("Chile")))
b <- sf::st_as_sf(a)
Error in `$<-.data.frame`(`*tmp*`, geometry, value = c(".........wkb.........)", : 
  replacement has 135 rows, data has 1
sessionInfo()
R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Gentoo Linux

Matrix products: default
BLAS:   /usr/lib64/libblas.so.3.12.0 
LAPACK: /usr/lib64/liblapack.so.3.12.0

locale:
 [1] LC_CTYPE=es_CL.utf8       LC_NUMERIC=C             
 [3] LC_TIME=es_CL.utf8        LC_COLLATE=es_CL.utf8    
 [5] LC_MONETARY=es_CL.utf8    LC_MESSAGES=es_CL.utf8   
 [7] LC_PAPER=es_CL.utf8       LC_NAME=C                
 [9] LC_ADDRESS=C              LC_TELEPHONE=C           
[11] LC_MEASUREMENT=es_CL.utf8 LC_IDENTIFICATION=C      

time zone: Chile/Continental
tzcode source: system (glibc)

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

other attached packages:
[1] sds_0.0.1.9006 terra_1.8-0   

loaded via a namespace (and not attached):
[1] compiler_4.4.1    Rcpp_1.0.13       codetools_0.2-20  countrycode_1.6.0
kadyb commented 4 weeks ago

It seems to "work" on Windows (no crash, but there is a message with error):

terra::as.polygons(terra::vect(sds::CGAZ(), query = sds::CGAZ_sql("Chile")))
class       : SpatVector 
geometry    : polygons 
dimensions  : 135, 1  (geometries, attributes)
extent      : -109.4555, -66.42422, -55.91544, -17.49834  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
Error in (function (cond)  : 
  error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': basic_string::_M_create
Session info ``` R version 4.4.1 (2024-06-14 ucrt) Platform: x86_64-w64-mingw32/x64 Running under: Windows 10 x64 (build 19045) Matrix products: default attached base packages: [1] stats graphics grDevices utils datasets methods base loaded via a namespace (and not attached): [1] compiler_4.4.1 tools_4.4.1 sds_0.0.1.9006 Rcpp_1.0.13 codetools_0.2-20 [6] countrycode_1.6.0 terra_1.8-0 gdal proj geos "3.8.4" "9.3.1" "3.12.1" ```
mdsumner commented 4 weeks ago

fwiw the output of CGAZ() and CGAZ_sql("Chile" ) is just

dsn <- "/vsizip//vsicurl/https://github.com/wmgeolab/geoBoundaries/raw/main/releaseData/CGAZ/geoBoundariesCGAZ_ADM0.zip"
qu <- "SELECT shapeGroup FROM geoBoundariesCGAZ_ADM0 WHERE shapeGroup IN ('CHL')"
vect(dsn, query = qu)

just as a workaround for @latot you don't need as.polygons, this is a shapefile in a zip read with a query

kadyb commented 4 weeks ago

I think @latot used as.polygons() because he was trying to convert from MULTIPOLYGON to POLYGON.

latot commented 4 weeks ago

Hi, I was trying originally get the polygons as sf from the raster, and in the middle I found this issue.

https://gis.stackexchange.com/questions/452160/fast-way-to-convert-raster-to-polygon-shapefile-in-r

At least to me, without as.polygons() it does not crash.

@mdsumner I'm looking to work with each raster cell as a sf polygon.

rhijmans commented 4 weeks ago

Fixed with https://github.com/rspatial/terra/commit/a0c09d9d35e120d5bffeb9fbf4cb2ce16eb3c65c

Your example is not with raster data. Treating a raster as polygons is almost never necessary, and almost always very inefficient

rhijmans commented 4 weeks ago

to convert from MULTIPOLYGON to POLYGON you can use disagg