4

It might look like this question is duplicated, but I'm asking about data extraction by points, rather than by polygons, and I couldn't find hints about point extraction. Therefore, please bare with me.

I am working with a massive amount of climate files whose values need to be extracted and converted to data frames in order to serve as a database for a Shiny application.

Specifically, I need to extract daily temperature values from a raster stack for 655 locations (points) and store those values in a data frame. Given the amount of data, I estimate the final data frame to be around 30 GB in size.

Since it is such a big amount of data, I am looking to use a multicore (parallel) approach that benefits from a quad core (8 threads) system. I also consider to process these data on an Amazon EC2 instance, so I need to speed up the code as much as I can in order to prevent wasting time (i.e. money).

I know that the raster package provides the multicore function beginCluster that might help me, but according to the documentation it only works when extracting data using polygons, not points.

Please take a look at the code that reproduces my procedure using a simplified version of my raster stack:

library(raster)

# Create date sequence
idx <- seq(as.Date("2010/1/1"), as.Date("2099/12/31"), by = "day")

# Create raster stack and assign dates
# WARNING: raster stack will be ~ 400 MB in size
r <- raster(ncol=5, nrow=5)
s <- stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)))))
s <- setZ(s, idx)

# Create random spatial points
pts <- SpatialPoints(cbind(x=runif(655, -180, 180),
                           y=runif(655, -90, 90)),
                     proj4string=CRS(projection(s)))

# Extract values to a data frame
e.df <- extract(s, pts, df=TRUE)

# Fix resulting data frame
DF <- data.frame(t(e.df[-1])); row.names(DF)=NULL
DF <- cbind(getZ(s), DF)
names(DF)[1] <- "Date";
names(DF)[2:length(names(DF))] <- paste0('Point ', 1:length(pts))

In the code above, clearly the slowest part of the process is the data extraction itself.

Is there any way to improve the speed of this data extraction procedure?

Any way to implement it under a multicore function?

  • 1
    I don't see where your 30Gb comes from. 89 years x 365 days x 655 stations is 21 million. You need to be handling about 1000 bytes for each of those units to get 30Gb, where one temperature measurement is 8 bytes. 160-250Mb maybe, which is easily stored in RAM. – Spacedman Aug 29 '17 at 10:10
  • 1
    This is "trivially parallel" over layers in your raster stack and over points in your spatial points. Use functions from the parallel package. – Spacedman Aug 29 '17 at 10:15
  • @Spacedman my actual raster stack is 100x95 at 0.5 deg resolution with 54750 "layers" (i.e. time slices). Also, I have 32 files like that that I need to extract data from (different climate scenarios and models). I just thought it was too big of an object for users to reproduce here. – thiagoveloso Aug 29 '17 at 10:46
  • @Spacedman I would also appreciate any more detailed suggestion regarding the parallel package. Not very familiar with multicore processing here... – thiagoveloso Aug 29 '17 at 10:48
  • 1
    Suggest you find some R parallel tutorials, things like this cran.r-project.org/web/packages/doParallel/vignettes/… might get you started. – Spacedman Aug 29 '17 at 10:55
5

I ended up using an approach based on the snowfall package. It is quite simple, works really good and the point extraction function is as fast as the number of cores that you can use. The approach I used was inspired by this post, and here is my reproducible example:

library(raster)
library(snowfall)

# Create date sequence
idx <- seq(as.Date("2010/1/1"), as.Date("2099/12/31"), by = "day")

# Create raster stack and assign dates
# WARNING: raster stack will be ~ 400 MB in size
r <- raster(ncol=5, nrow=5)
s <- stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)))))
s <- setZ(s, idx)

# Create random spatial points
pts <- SpatialPoints(cbind(x=runif(655, -180, 180),
                           y=runif(655, -90, 90)),
                     proj4string=CRS(projection(s)))

# Extract values to a data frame - multicore approach
# First, convert raster stack to list of single raster layers
s.list <- unstack(s)
names(s.list) <- names(s)

# Now, create a R cluster using all the machine cores minus one
sfInit(parallel=TRUE, cpus=parallel:::detectCores()-1)

# Load the required packages inside the cluster
sfLibrary(raster)
sfLibrary(sp)

# Run parallelized 'extract' function and stop cluster
e.df <- sfSapply(s.list, extract, y=pts)
sfStop()

# Fix resulting data frame
DF <- data.frame(t(e.df)); row.names(DF)=NULL
DF <- cbind(getZ(s), DF)
names(DF)[1] <- "Date";
names(DF)[2:length(names(DF))] <- paste0('Point ', 1:length(pts))

# Check resulting data frame
head(DF)

I hope it can serve as a future reference for people interested in doing the same.

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.