Tuesday, 7 April 2015

Combine multiple partially overlapping rasters into a single raster in R


I have about 1,000 rasters which are around 10km x 10km and partially overlap to cover a whole country.


I would like to convert these into a single national raster, and where raster overlap, add the cell values together.


Currently, I am using ArcGIS's workspace to new raster tool, but I would like to replicate the process in the R language.


I've tried a couple of things but usually end up with only the intersection of the rasters. Many of the solutions on StackExchange advocate cropping to a common area, which is the exact opposite of what I want.


IS there a general solution for aligning two rasters and then making a mosaic which has the combined extent of both rasters?


Can this method be scaled to work for thousands of rasters?


EDIT


> sessionInfo()

R version 3.3.2 (2016-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=English_United Kingdom.1252 LC_CTYPE=English_United Kingdom.1252 LC_MONETARY=English_United Kingdom.1252
[4] LC_NUMERIC=C LC_TIME=English_United Kingdom.1252

attached base packages:
[1] stats graphics grDevices utils datasets methods base


other attached packages:
[1] rgdal_1.2-5 raster_2.5-8 sp_1.2-4

loaded via a namespace (and not attached):
[1] tools_3.3.2 Rcpp_0.12.9 grid_3.3.2 lattice_0.20-34

> rasterOptions()
format : raster
datatype : FLT8S

overwrite : FALSE
progress : none
timer : FALSE
chunksize : 1e+07
maxmemory : 1e+08
tmpdir : D:/RTemp\RtmpgzISfL/raster/
tmptime : 168
setfileext : TRUE
tolerance : 0.5
standardnames : TRUE

warn depracat.: TRUE
header : none

Answer



The best solution to this is making a list of the rasters, then passing this to a function based on the apply family


The following code was pulled from a similar question wrapped into a function and should work for you


mosaicList <- function(rasList){

#Internal function to make a list of raster objects from list of files.
ListRasters <- function(list_names) {
raster_list <- list() # initialise the list of rasters

for (i in 1:(length(list_names))){
grd_name <- list_names[i] # list_names contains all the names of the images in .grd format
raster_file <- raster::raster(grd_name)
}
raster_list <- append(raster_list, raster_file) # update raster_list at each iteration
}

#convert every raster path to a raster object and create list of the results
raster.list <-sapply(rasList, FUN = ListRasters)


# edit settings of the raster list for use in do.call and mosaic
names(raster.list) <- NULL
#####This function deals with overlapping areas
raster.list$fun <- sum

#run do call to implement mosaic over the list of raster objects.
mos <- do.call(raster::mosaic, raster.list)

#set crs of output
crs(mos) <- crs(x = raster(rasList[1]))

return(mos)
}

To use this do something along the lines of:'


raster_files <- list.files(path ="folder_with_files",pattern = ".tif$",full.names = TRUE )

national_layer <- mosaicList(raster_files )

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...