rpoint-in-polygonterra

point in polygon std::bad_alloc error in terra

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.

Solution

  • 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.