3

I am trying to get crop area for admin units of whole African continent using the data available here at 30-m spatial resolution globally (raster value 1 = crops). I've tried both terra and raster packages in R and tried to parallelize the approach but couldn't succeed. A few posts stated that parallel computing is not supported for terra as of yet. However, I tried to use doParallel and foreach to parallelize the code but couldn't succeed. I have also made unsuccessful attempts to use raster with parallel package that I am not posting here for the sake of shortening my post.

So far, I have written the following code that works great for admin units of 1 country (166 polygons but overall area of the country is around 15 times smaller than that of actual shapefile) only but when I try to run it for my actual file containing 450+ polygons, my disk space (100 GB free) runs out quickly and I get Error: [classify] insufficient disk space (perhaps from temporary files?). Also, the process is too slow and it took around half an hour to do the calculations for 15 polygons for a single raster file only (while I have 450+ polygons in the shapefile that vary greatly in size). Can an appropriate method be used to do this through parallel computing or some other more efficient method that is memory and disk space friendly? As my AOI lies in different coverages, therefore, I am using vrt files (instead of actual mosaics) to combine all coverages into 1 vrt file for each year.

library(terra)
gfile <- ***<<list of vrt files path>>***
s <- ***<<SpatVector>>***
for(f1 in gfile){    
  r <- rast(f1)
  y <- substr(basename(f1), 5, 8) # Extract name of year
  for (i in 1:length(s)){        
    rm <- crop(r, s[i], mask = TRUE)
    if(sum((as.data.frame(rm)))){
      r1 <- ifel(rm == 1, 1, NA)
      if(sum((as.data.frame(r1)))){
        s[i,paste0('crop',y)] <- sum(as.data.frame(cellSize(r1)))/10000
      }else{
        s[i,paste0('crop',y)] <- 0
      }
    }else{
      s[i,paste0('crop',y)] <- 0
  }
  }
}
writeVector(s, outshp, overwrite=TRUE)
1
  • This isn't going to run in parallel, but you might want to see how exactextractr::exact_extract(crops, admins, 'frac', coverage_area = TRUE) performs. It would require the GitHub version of exactextractr.
    – dbaston
    May 30, 2022 at 1:05

1 Answer 1

5

You could write a function like this

fun <- function(rastfile, vectfile) {
    library(terra)
    r <- rast(rastfile)
    s <- vect(vectfile)
    rs <- crop(r, s, mask = TRUE)
    if (minmax(rs)[2] > 0) {
        rs <- cellSize(rs, unit="ha") * rs
        global(rs, "sum", na.rm=TRUE) |> unlist()
    } else {
        0
    }
}

And run that in parallel. The function can be sent to the different nodes, and as it uses filenames there are no serialization issues. You should use a filename (of the vrt) for the raster data, but you could alternatively send a wraped SpatVector for each polygon.

The bottleneck is the masking (the rasterization of the polygons), but the approach I show should be faster and safer than what you had by using minmax instead of as.data.frame and because * is faster than ifel.

If it is still too slow, then aggregating the input (with "sum") would help a lot, and would not change the results much if there are many cells in a polygon.

Here is more discussion on parallel computing with terra

3
  • I really appreciate all this help. I'll apply the method you shared later today and see how it goes before I mark this answer as accepted. I did try to convert my code into function and parallelize it after going through the same link you shared but It all makes more sense now and I hope I'll now be able to get it done. May 30, 2022 at 0:15
  • BTW you meant if (minmax(rs)[2] > 0) instead of x, I guess. Is that right? May 30, 2022 at 0:15
  • that's correct, and fixed now May 30, 2022 at 0:31

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.