I have two rather big spatial polygon datasets for 250m x 250m gridded data. Although they use an identical grid, they do use different grid identifiers. I'd like to join them to have a lookup between the different identifiers.

I can't share the data and I don't think I could create fake data that results in the same issue. So here are some examples:

demographics <- st_read("../data/RTTK/rttk_250m_2019.shp")
## Reading layer `rttk_250m_2019' from data source `C:\Users\***\rttk_250m_2019.shp' using driver `ESRI Shapefile'
## Simple feature collection with 632961 features and 108 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 83750 ymin: 6629000 xmax: 728000 ymax: 7776250
## Projected CRS: ETRS89_TM35FIN_E_N

size in R workspace = 1GB

  landcover <- st_read("../data/Corine/T11_grd_e_250m.shp")
  ## Reading layer `T11_grd_e_250m' from data source `C:\Users\***\T11_grd_e_250m.shp' using driver `ESRI Shapefile'
  ## Simple feature collection with 6262152 features and 7 fields
  ## Geometry type: POLYGON
  ## Dimension:     XY
  ## Bounding box:  xmin: 61500 ymin: 6605750 xmax: 733000 ymax: 7776500
  ## Projected CRS: ETRS89 / TM35FIN(E,N)

size in R workspace = 5.1GB

  ## 5.15 GB

*Using centroids of the demographics polygons makes the demographics object a bit smaller and makes it easier to join (always lies within the landcover polygon).

d_cent <- st_centroid(demographics)

... and predictably I end up with this error message:

test <- landcover %>%
  join = st_intersects)
## Error: cannot allocate vector of size ***

Here's a very small example of the data grid. The blue numbers are the grid ID's from the demographics file (smaller object, as polygons without buildings are excluded), the black numbers are grid ID's from the landcover file.

enter image description here

There's no variable that could be used for subsetting the landcover and as you can see from the image above, there's no immediate apparent logic in the grid ID assignment:

## [1] "ID"         "GRIDCODE"   "x"          "y"          "xyind"      "Shape_Leng"
## [7] "Shape_Area" "geometry"

I've already switched to Microsoft R Open, which does help a bit, but it's not enough.
(There's no research money to buy a better laptop. QGIS also regularly crashes.)

System Info:

Microsoft R Open 4.0.2 R version 4.0.2 (2020-06-22) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19042)

Processor Intel(R) Core(TM) i7-6700 CPU @ 3.40GHz 3.41 GHz Installed RAM 16.0 GB (15.9 GB usable) System type 64-bit operating system, x64-based processor

  • 1
    I used to work with large data (up to 5GB) and the old increasing memory limit trick worked just fine for me though it took some time to process( Windows 64 bit 8G RAM) ###checking memory size/limit gc() memory.size() memory.limit() #increase the limit memory.limit(size=56000)
    – Es_a
    Jun 17, 2021 at 7:05
  • @Esmaeel I've been using gc(), but memory.limit() doesn't seem to do much (RStudio) Jun 17, 2021 at 7:31
  • 1
    I'm not sure if it'll work in your case, It suppose to trick R about the size of your memory by allocating something like 50000, it did not help make things faster for sure, but it solved the same error in my case
    – Es_a
    Jun 17, 2021 at 7:44
  • 2
    You mentioned that the grids are identical in terms of their geometry? If so, you can calculate X and Y columns with the centroid coordinates for each of the grids, then do an ordinary (non-spatial) join with by=c("X", "Y"). (Possibly rounding the values to make sure they match.) This should be much more efficient than any kind of spatial join. Jun 17, 2021 at 7:44
  • 1
    @MichaelDorman that's a really interesting idea - thanks for that Jun 17, 2021 at 8:08

1 Answer 1


I don't have a solution for the whole dataset, but it is possible to run a subset without too much problems as follows:

  1. set memory.limit (it worked after rebooting my machine)
  1. Subset the data to the absolute minimum you want to use
# Helsinki municipalities
demo.HR <- demographics %>%
  filter(kunta %in% c('049','091','092','235'))
# or use st_crop
  1. Calculate centroids
d_cent <- st_centroid(demo.HR)
  1. Join grids
test <- landcover %>%
          st_join(d_cent, join = st_intersects)

# write to geopackage if you need to use the data in the future

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.