1

I am trying to do a point in polygon analysis This is my polygon

my_poly

class       : SpatVector 
geometry    : polygons 
dimensions  : 2791790, 35  (geometries, attributes)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
names       : basin_id lat_id lon_id l1     l2     l3    l4    l5     l6   l7 (and 25 more)
type        :    <num>  <num>  <num> <int> <int> <int> <int> <int> <int> <int>              
values      : 3.12e+09    654    407    40    54    43    55    53    50    53              
              8.12e+09    630    895    41    40    47    47    50    52    46              
              3.12e+09    598    267    37    37    47    39    42    50    45 

# this is my point 
grid_pts              
class       : SpatVector 
geometry    : points 
dimensions  : 9177, 1  (geometries, attributes)
extent      : -76.35417, -75.25417, 44.9625, 45.52917  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
names       :  cell
type        : <num>
values      :     1
                  2
                  3               
pts_basin <- terra::extract(my_poly, grid_pts) 

Error in x@ptr$relate_between(y@ptr, relation) : std::bad_alloc

The polygon covers the entire globe while my point shapefile covers a small area. I only need the basin_id for each point.

89_Simple
  • 3,393
  • 3
  • 39
  • 94
  • 1
    Sounds like you’re out of memory. – user438383 Aug 25 '22 at 12:45
  • 1
    Maybe buffer your grid_pts extent and use that to extract your basin roi from poly, then do the terra::extract. [and sometimes std::bad_alloc isn't out of memory](https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2019/p1404r1.html), but probably simpler to assume it does. – Chris Aug 25 '22 at 14:20
  • Could you comment when you've filed bug report at `terra`, interested in this issue and would like to follow. – Chris Aug 26 '22 at 13:11
  • Okay. How do I do it? I don't even know if this is a bug or my bad programming practice. If it is and if you would like me to file it, I can do it if you can share/show me how to do it – 89_Simple Aug 26 '22 at 13:45
  • If the pkg author says it is a bug, it is. And just like here, well, in addition to what's presented above, link to the data sets used, then code and result. `terra` is supposed to help prevent these issues by reporting meta-data, mapping larger than memory files on disk, and only using what's needed in an operation. It appears some more tuning would be required to prevent 'asking for too much', perhaps warn & suggest subsetting as below...[terra bug reports](https://%20https://github.com/rspatial/terra/issues), then new issue. Actually very useful and exciting to trip on a real bug. Chat? – Chris Aug 26 '22 at 14:41
  • Actually, my above Chat? is a dumb question (on my part) as I've never gotten there but for making way to many comments until the suggestion 'Continue in Chat` is triggered...But also, if you navigate over to the open issues, and read through the ones that note `comments` as a number on the right, you'll see a discussion typically of the kinds of info useful to skewer a bug. – Chris Aug 26 '22 at 14:50
  • did a git clone of [terra](https://github.com/rspatial/terra), then R CMD build terra, R CMD INSTALL terra, and have terra 1.6.11; suggest you do the same to continue with your interesting research. Were you basins from [HydroSHEDS](https://www.hydrosheds.org/)? – Chris Aug 27 '22 at 17:51
  • Yes my basin was from HydroSHEDS. https://data.hydrosheds.org/file/technical-documentation/HydroATLAS_TechDoc_v10_1.pdf – 89_Simple Aug 28 '22 at 01:24

1 Answers1

3

This clearly is an out-of-memory issue as there are 2791790 * 9177 ~ 25 billion combinations to be evaluated. I think you can avoid that by first subsetting the polygons. Like this:

Example data

library(terra)
pol <- vect(system.file("ex/lux.shp", package="terra"))
set.seed(1)
pts <- spatSample(as.polygons(pol, ext=T), 10)

Solution

i <- is.related(pol, pts, "covers")
extract(pol[i], pts)
#   id.y ID_1       NAME_1 ID_2       NAME_2 AREA   POP
#1     1    3   Luxembourg    8     Capellen  185 48187
#2     2    2 Grevenmacher   12 Grevenmacher  210 29828
#3     3    3   Luxembourg    8     Capellen  185 48187
#4     4    1     Diekirch    3      Redange  259 18664
#5     5    2 Grevenmacher   12 Grevenmacher  210 29828
#6     6   NA         <NA>   NA         <NA>   NA    NA
#7     7   NA         <NA>   NA         <NA>   NA    NA
#8     8    1     Diekirch    3      Redange  259 18664
#9     9    2 Grevenmacher    6   Echternach  188 18899
#10   10    1     Diekirch    3      Redange  259 18664

Same result as without subsetting

extract(pol, pts)
#   id.y ID_1       NAME_1 ID_2       NAME_2 AREA   POP
#1     1    3   Luxembourg    8     Capellen  185 48187
#2     2    2 Grevenmacher   12 Grevenmacher  210 29828
#3     3    3   Luxembourg    8     Capellen  185 48187
#4     4    1     Diekirch    3      Redange  259 18664
#5     5    2 Grevenmacher   12 Grevenmacher  210 29828
#6     6   NA         <NA>   NA         <NA>   NA    NA
#7     7   NA         <NA>   NA         <NA>   NA    NA
#8     8    1     Diekirch    3      Redange  259 18664
#9     9    2 Grevenmacher    6   Echternach  188 18899
#10   10    1     Diekirch    3      Redange  259 18664

This bug has now been fixed in terra version 1.6-11.

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63