4

I am trying to project a large raster to UTM (meters), but I keep crashing my laptop (16GB) due to a memory shortage. It seems like using projectRaster tends to deplete the machine's memory, and some solutions have been suggested in questions like this one. However, even using the suggested solution my laptop keeps crashing. The raster in question is a mosaic that I create from 4 DEM tiles downloaded from USGS with a resolution of 1/3 arc second (~10m) that correspond to the city of Milwaukee, WI. This mosaic is 396.8 Mb. The process goes as follows and can be run in any machine to download the rasters:

library(raster)

rasterOptions(memfrac = 0.3) # the solution suggested in several questions to limit access to memory by the raster package.

import the four rasters via their individual link.

r1 <- raster("https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n43w089/USGS_1_n43w089.tif")
r2 <- raster("https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n44w088/USGS_1_n44w088.tif")
r3 <- raster("https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n43w088/USGS_1_n43w088.tif")
r4 <- raster("https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n44w089/USGS_1_n44w089.tif")

# create a list with the four rasters and merge them into one single raster layer
rl <- list(r1, r2, r3 , r4)
mosaic <- do.call(merge, rl)

rm(r1,r2,r3,r4, rl) # remove unnecessary files in memory to save space

# project the mosaic layer to UTM zone 16N (m) -- THIS IS WHERE THE COMPUTER CRASHES
mosaic <- projectRaster(mosaic, crs = "+proj=utm +zone=16 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")

In this other question, rasterOptions is used to set chunksize = 1e+04 and maxmemory to 1e+0.6. I tried doing this and my machine didnt complete the line mosaic <- do.call(merge, rl) in the entire night.

I have considered the option of projecting each raster (r1, r2, r3, r4) separately before creaing the mosaic, but I am weary of the issues that this may cause at the time of merging them together, and I would very much rather find a memory efficient way to handle the reprojection than just avoiding it with a by-pass.

4
  • 2
    It is probably better to use GDAL for this task. If I have time later I can write out a complete answer but have a look at gdalUtils and the function gdalwarp, that should handle the the task without any issues
    – JonasV
    Jan 7 '21 at 8:51
  • @JonasV That seems to work like a charm in less than 4 minutes! but it created a new issue. I used gdalwarp("my_raster_file", dstfile = "path_to_save_new_raster", r = "bilinear", t_srs = "+proj=utm +zone=16 +ellps=GRS80 +datum=NAD83 +units=m +no_defs", overwrite = TRUE). The issue is that the resolution of the raster changes from ~10m to 27m. I guess I will have to define the new resolution in the line? If you see any other red flag in this command or something I am missing, please let me know. Thank you for your help. Jan 7 '21 at 9:51
  • I added tr = c(10,10) to the function and it returned the resolution that I needed. Cheers! Jan 7 '21 at 9:54
  • 2
    Nice to see that it works! I have only one note: It is preferable to use ESPG Codes wherever possible since PROJ strings are on their way out. So instead you could use t-srs = "EPSG:26916". Right now it would probably not make a difference but it's a good habit to pick up for the future.
    – JonasV
    Jan 7 '21 at 13:35
4

I am posting an answer based on the comment by JonasV, which worked perfectly, and may be of help for other users that encounter the (common) problem of projecting large raster datasets in R. Basically, the solucion has been to use GDAL instead of the raster package:

library(gdalUtils)
library (raster)

r1 <- raster("https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n43w089/USGS_1_n43w089.tif")
r2 <- raster("https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n44w088/USGS_1_n44w088.tif")
r3 <- raster("https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n43w088/USGS_1_n43w088.tif")
r4 <- raster("https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n44w089/USGS_1_n44w089.tif")


rl <- list(r1, r2, r3 , r4)
mosaic <- do.call(merge, rl)

writeRaster(mosaic, "path_to_mosaic_raster")

gdalwarp("path_to_mosaic_raster", 
         dstfile = "path_to_save_projected_raster", 
         r = "bilinear", 
         t_srs = "EPSG:26916", 
         overwrite = TRUE,
         tr = c(10,10)) ## to ensure that the projected raster has the desired resolution
1
  • 1
    Man this package is fast, projectRaster took me a whole day to reproject a big file and it came to a point where my C drive started running short of disk space. But, gdalwarp did the whole thing in a few minutes with no memory and storage problems. Thank you for posting this question. May 25 '21 at 7:09

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.