6

I've noticed that running extract takes very different amounts of time depending on whether or not the raster is in memory, on disk in the Native format (.grd), or on disk in another format (e.g., GeoTiff). The following MWE demonstrates:

library(rgdal)
library(raster)
# Make a big raster
rst <- raster(ncol=10000, nrow=10000, xmn=-1000, xmx=1000, ymn=-100, ymx=900)
values(rst) <- runif(ncell(rst))

# Make a big set of points
ext <- extent(rst)
points <- data.frame(lng=runif(100000, ext@xmin, ext@xmax),
                     lat=runif(100000, ext@ymin, ext@ymax))

# Save raster to disk as grd (Native) and tif (GeoTIff)
writeRaster(rst, "tmp_rast.grd")
writeRaster(rst, "tmp_rast.tif")

# Load rasters from disc
rst.grd <- raster("tmp_rast.grd")
rst.tif <- raster("tmp_rast.tif")

# Test extract time for rasters in memory, from native format, and from tif
system.time(test <- extract(rst, points)) # <1s on my machine
system.time(test <- extract(rst.grd, points)) # 6s on my machine
system.time(test <- extract(rst.tif, points)) # 88s on my machine

This isn't too much of a surprise, given that the raster docs note that the Native format is binary. But I'm wondering if there are any useful workarounds here if the input file is in GeoTIFF format (outputted from gdalwarp), assuming the raster is too big to load in memory. Is there a fast way to convert it to binary on the fly?

2
  • 3
    You can writeRaster(rst.tif, "rst_native.grd") to convert the file format. Extract is slow for tiled files, it's a known issue with a pending fix. You can use gdalinfo to get information about whether your file is tiled.
    – mdsumner
    Nov 23, 2016 at 1:30
  • @mdsumner Thanks, that does save time. For the example above, resaving takes another 12s, making that method about 70s faster overall. I don't believe my file is tiled, though.
    – pbaylis
    Nov 26, 2016 at 19:36

2 Answers 2

2

Even though it is in beta (on CRAN) I would recommend migrating to the terra package. The terra package is from the same developer as raster and is considered an eventual replacement for raster. The functions/code have mostly been moved over to C++ and exhibit massive performance and speed gains, including in extract. One thing to note is that vector (points/polygons/lines) objects need to be coerced into the vect class. The raster/stack/brick classes have been replaced by a single rast object class (regardless the number of bands). I have been having good luck with nested "on the fly" coersion eg., extract(x, vect(y)).

2
  • Thanks, I didn't know about terra. I have been use exactextractr for polygonal extraction recently, but (weirdly?) it doesn't allow point extraction right now. I'll start playing around with terra. I wonder if the author has any plans to allow it interface with sf, which has become my (and many others') standard shapefile manipulation tool.
    – pbaylis
    Oct 3, 2020 at 21:46
  • 2
    @pbaylis not sure what you mean by "interface with sf". Currently, coercion (eg., terra::vect(x)) as well as some generalizable functions (eg., terra::ext(x)) is supported for sf class objects. So, for raster extract using an sf class object one would just use terra::extract(x, terra::vect(sf_object)) Oct 3, 2020 at 23:37
0

If you're looking for some speed increases, you could try the velox package. It requires your raster data to be converted into a velox format, and your points to be in a SPATIALPOINTS or sf format. Here's an example building on the code you have provided:

sf_points <- sf::st_as_sf(points, coords = c("lng", "lat"))
vx <- velox::velox(rst)
system.time(test <- vx$extract_points(sp = sf_points)) # 0.06s on my machine
1
  • 1
    At the time of this post, velox was a great alternative. However, it is now a dead package and no longer on CRAN. It only compiles from GitHub on certain platforms. For polygon extraction, I would reccomend the exactextractr package. Also see my post on the terra package. Oct 2, 2020 at 17:52

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.