4

I have sf data with single point geometries (points) and a 23,000 x 18,000 elevation RasterLayer with a single band, loaded from a geotiff (elev). They have the same CRS. I want to extract the elevation values to the points.

When I try running extract(elev, points), since my raster is very large I get:

Error: cannot allocate vector of size 3.1 Gb

I followed advice here R - multicore approach to extract raster values using spatial points and tried to parallelize, but I only have a single raster, so chances are the work won't get divided. Still I try the following:

library(snowfall)
sfInit(parallel=TRUE, cpus=parallel:::detectCores()-1)
sfLibrary(raster)
sfLibrary(sp)

# Run parallelized 'extract' function and stop cluster
e.df <- sfSapply(elev, raster::extract, y=points)
sfStop()

I get the same result:

Error: cannot allocate vector of size 3.1 Gb

All the examples I've found so far, such as the above, store multiple rasters in a list. How do I parallelize extraction from a single raster like this? Or is there a way to make use of sparse matrices and memory management as noted here?

points covering extent of elev raster

I'm on Windows 10 running RStudio 1.1.419 (64 bit). 8GB RAM. Intel Core i5-6200U CPU @ 2.3GHz 2.4GHz

  • 1
    There is something not tracking here. I have extract ~1M points from a 2GB single-band raster and did not even touch the top of memory on my laptop. Most raster functions operate out of memory so do not read in a raster that large. The memory consumption should be around the raster block read in and the resulting vector associated with extract. For polygons I could easily see memory issues but just not with points. There is seemingly something else going on here. I would recommend trying the velox package to see if you get a memory error. – Jeffrey Evans Feb 26 at 15:44
1

You can split your image to tiles like 100*100 tiles, then loop tiles to extract points. Or you can try database solution with building spatial index. When you load your image to spatial database, it also can tile your image. PostGIS is open source, in PostGIS just a SQL query can do this job, the function is ST_Value.

[EDIT]You can loop records(here are tiles) using PL/pgSQL like below:

CREATE OR REPLACE FUNCTION work_on_tiles() AS
$BODY$
BEGIN
    FOR tile IN
        SELECT rast FROM public.elev WHERE rid > 0
    LOOP
        -- can do some processing here
        var = ST_Value(rast)
    END LOOP;
    RETURN;
END
$BODY$
LANGUAGE plpgsql;

SELECT * FROM work_on_tiles();
  • 1
    Thanks! I tried using a version of the raster I loaded in postgis, but I'm not sure how it works. raster2pgsql -I -C -s 26985 -t 182x226 "~/elev.tif" public.elev | psql -d postgres. After reading the resulting elev table from my psql connection, I don't know how to loop over tiles, so I'm probably only getting one – Wassadamo Feb 26 at 8:09
  • @Wassadamo you have almost worked it out. You can use PLpgsql language to loop table records(here is tiles) to apply ST_Value. – dorbodwolf Feb 27 at 8:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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