Cyclic bug for calculating accumulated cost rasters in gdistance R

Dear HPC R-users,

I’m attempting to run a loop process on shared-bigmem. In the code below you can find my R script and my sbatch script. I ran this code on a smaller scale locally on my macOS Big Sur version 11.4 and it succeeded.

The goal of the script is to calculate the accumulated cost of traversing a landscape. This means I overlay coordinates/points on top of a so called transition layer. The expected output are approximately 4000 raster files, because that’s the number of points/coordinates I have and I want to calculate a separate cost layer for each.

Unfortunately on Baobab something goes wrong. Once in every two coordinates I get two empty rasters (values are inf until -inf). Since the process seems cyclic, I’m thinking it has something to do with the parallelization. I can’t figure out what it is. I’ve tried solving the issue for the last two weeks. The exact same R code produces the right output on my laptop.

In the code box below, you can see that once in every two files, the rasters have exactly the same size, namely 11344427. This file size corresponds to an empty raster.

I tried running the code on 2 cpus with each 220gb and on 1cpu with 480gb. Both produced the same bug.

Any help or guidance would be much appreciated!

SBATCH script

#SBATCH -J traveltime
#SBATCH -e traveltime.e%j
#SBATCH -o traveltime.o%j
#SBATCH -n 1
#SBATCH --cpus-per-task 1
#SBATCH --mem=480000
#SBATCH -p shared-bigmem
#SBATCH -t 12:00:00

module load GCC/10.2.0  OpenMPI/4.0.5
module load GDAL/3.2.1
module load PROJ/7.2.1
module load R/4.0.4

srun Rscript  01_accumulated_cost.R

R script

# libraries

# data to load in
facilities <- st_read("/home/users/h/hierink/philippines/data/raw/vHealth_facilities/healthfacilities_final.shp")
friction_layer <- raster("/home/users/h/hierink/philippines/data/raw/rFriction/landcover_friction.tif")
#transition_gc <- readRDS("/home/users/h/hierink/philippines/data/processed/transition_philippines.RDS")
print("data loaded")

# reproject raster
projected_friction <- projectRaster(friction_layer, crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

# transition layer divided by resolution to bring it to cost per m
friction_layer <- friction_layer/res(friction_layer)

# create transition layer for each country
transition <- transition(friction_layer, function(x) 1/mean(x), 8)
transition_gc <- geoCorrection(transition)

# facility transformation
facilities <- st_transform(facilities, crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
facilities <- as(facilities, "Spatial")
points <- coordinates(facilities)

# to try with the facilities
facilities_try <- points[1:4,]

# set the number of cores in the sbatch script
print("step registerparallel")

# print number of workers

# write for loop to add all values in raster with each other
system.time(foreach(i = 1:nrow(facilities_try), .packages= c("raster", "sf", "gdistance")) %dopar% {
  transition <- transition_gc
  facility <- facilities_try[i,]
  traveltime <- accCost(transition, facility)
  writeRaster(traveltime, paste0("/home/users/h/hierink/philippines/data/processed/rTravel_time/travel_time_", i, ".tif"), overwrite = T)})


-rw-r--r-- 1 hierink hpc_users 11344427 Jun 17 15:16 travel_time_10.tif
-rw-r--r-- 1 hierink hpc_users 11344427 Jun 17 15:19 travel_time_11.tif
-rw-r--r-- 1 hierink hpc_users 11344427 Jun 17 15:22 travel_time_12.tif
-rw-r--r-- 1 hierink hpc_users 11344427 Jun 17 15:25 travel_time_13.tif
-rw-r--r-- 1 hierink hpc_users 11344427 Jun 17 15:28 travel_time_14.tif
-rw-r--r-- 1 hierink hpc_users 11344427 Jun 17 15:32 travel_time_15.tif
-rw-r--r-- 1 hierink hpc_users 11344427 Jun 17 15:35 travel_time_16.tif
-rw-r--r-- 1 hierink hpc_users 11344427 Jun 17 15:38 travel_time_17.tif
-rw-r--r-- 1 hierink hpc_users 11851011 Jun 17 14:47 travel_time_1.tif
-rw-r--r-- 1 hierink hpc_users 11344427 Jun 17 14:50 travel_time_2.tif
-rw-r--r-- 1 hierink hpc_users 11344427 Jun 17 14:53 travel_time_3.tif
-rw-r--r-- 1 hierink hpc_users 11849362 Jun 17 14:56 travel_time_4.tif
-rw-r--r-- 1 hierink hpc_users 11849695 Jun 17 15:00 travel_time_5.tif
-rw-r--r-- 1 hierink hpc_users 11344427 Jun 17 15:03 travel_time_6.tif
-rw-r--r-- 1 hierink hpc_users 11344427 Jun 17 15:06 travel_time_7.tif
-rw-r--r-- 1 hierink hpc_users 41927453 Jun 17 15:09 travel_time_8.tif
-rw-r--r-- 1 hierink hpc_users 11344427 Jun 17 15:12 travel_time_9.tif


do you have the same issue with an older R version? Are you using the same R version on your laptop?


Hi Yann,

Locally I run R version 4.0.3 and on the cluster I run 4.0.4. I will try it today with R 3.6 on the cluster.


We do have 4.0.0 on the cluster, try with this one maybe?