Troubleshooting work flow for 10m Sentinel data at large scales (contiguous USA) in [R] using [terra]

Question

I am working with ESRI's Sentinel landcover data at 10meter resoltuion: https://livingatlas.arcgis.com/landcoverexplorer/#mapCenter=-83.21%2C34.332%2C4&mode=step&timeExtent=2017%2C2021&year=2017&downloadMode=true

I'm working with these data in R with the package terra

My project extent is the entire continental USA. I have binned the USA into equal area hexagonal bins

(projection: proj<-"+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs")

of varying spatial scales and I need to intersect the landcover data with these bins (a polygon file) and extract landcover information (relative areas etc.) for each bin. Because the bins can overlap with two or more images depending on its position I am currently reprojecting all the sentinel images to the same equal area projection as my bins but this is taking a very long time given the resolution of the data. Because of my project's focus (insects) the higher the resolution of the images, the better. Additionally, because of the time-varying aspect of my data I need to do this for every year from 2017-2022.

This is just an example code with terra but this is where I am hitting a slowdown because one run of reprojecting 37 rasters takes anywhere from 12-15 hours on my computer:

library(terra)

#get directory
raster_dir<-"/Users/leoohyama/Downloads/USA_Sentinel_2022/"
file_Rast<-list.files(raster_dir)
setwd(raster_dir)

#set projection
proj<-"+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"

#lapply to reproject all raster images and write them to disk
lapply(file_Rast, function(x){
  rast1<-terra::rast(x)
  ras1p<-terra::project(rast1, proj)
  terra::writeRaster(ras1p, paste(x, "projected.tif", sep = "_"))
  rm(rast1)
  rm(ras1p)
})

Once I reproject the images per year I plan to group then in virtual rasters and mosaic them. Then I plan to run the intersections with my hex bins on the mosaic-ed raster. But I also know this will take forever (or I may run out of memory in terms of RAM).

However one other issue I am also realizing is that the reprojected rasters for just one year take up close to 40-50 Gigabytes of space. So the issue of memory usage is also apparent.

Is there a more efficient way to do this?

Answer 1

Don't reproject the rasters! Project the hexagons (which are a polygon vector data set) to the raster CRS. Then do the overlay. Then if you need, reproject the hexagons back to their original projection.

Generally, reprojecting vector data back-and-forth is lossless and you get back what you put in. Raster reprojection is lossy and for rasters of the size you have, very slow.

Raster reprojecting is usually only necessary when you have two rasters of different projections that you need to relate together, and even then sometimes you can avoid reprojections by considering one raster as points or polygons.

Answer 2

You do not provide any code, so it is hard to guess what the bottlenecks in your code might be. It might be helpful if you provide further details.

There is also an R package called mapme.biodiversity specifically designed for such use cases. Unfortunatley, the ESRI landcover product is currently not implemented, but there is the ESA Landcover product and also many other indicators. The code can easily be run in parallel achieving good performances for large spatial extents (there are people using it to e.g. extract forest cover for all of the Americas in reasonable time).

Disclaimer: I am the original author of the mapme.biodiversity package.