1

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")

Example of my hexbin polygon

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?

2 Answers 2

2

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.

2
  • Got it! But what about polygons that overlap with mulitple rasters? How would I solve that issue? Just identify partial overlapping ones and combine them downstream?
    – Leo Ohyama
    Oct 19, 2023 at 16:43
  • 2
    Yes, maybe compute bounding rectangles for each raster, work out intersections, then do the overlay op once for hexes in one rectangle and N times for any in more. Then add areas (or compute percentages). Can also run parallel over rasters easily too...
    – Spacedman
    Oct 19, 2023 at 18:17
0

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.

2
  • Hi. I added some code (I have yet to get to the mosaic stage just yet). But I hope this helps better understand the problem. Your package looks useful but unfortunately I am still pretty tied to the Sentinel 10m rasters
    – Leo Ohyama
    Oct 19, 2023 at 13:34
  • It is usually more performant to project the vector geometries to the CRS of the rasters before extracting. On another note, mapme.biodiversity is extensible so you could add the ESRI Landcover resource and/or open an issue requesting the addition. Oct 19, 2023 at 13:42

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.