Wednesday 2 April 2014

Merge .ASC grids with R

A couple of years ago I found online a script to merge several .asc grids into a single file in R.
I do not remember where I found it but if you have the same problem, the script is the following:

 setwd("c:/temp")  
 library(rgdal)  
 library(raster) 

 
 # make a list of file names, perhaps like this:  
 f <-list.files(pattern = ".asc")  


 # turn these into a list of RasterLayer objects  
 r <- lapply(f, raster) 

 
 # as you have the arguments as a list call 'merge' with 'do.call'  
 x <- do.call("merge",r) 

 
 #Write Ascii Grid  
 writeRaster(x,"DTM_10K_combine.asc")  

It is a simple and yet very affective script.
To use this script you need to put all the .asc grids into the working directory, the script will take all the file with extension .asc in the folder, turn them into raster layers and then merge them together and save the combined file.

NOTE:
if there are other file with ".asc" in their name, the function
list.files(pattern = ".asc") 

will consider them and it may create errors later on. For example, if you are using ArcGIS to visualize your file, it will create pyramids files that will have the same exact name of the ASCII grid and another extension.
I need to delete these file and keep only the original .asc for this script and the following to work properly.



A problem I found with this script is that if the raster grids are not properly aligned it will not work.
The function merge from the raster package has a work around for this eventuality; using the option tolerance, it is possible to merge two grids that are not aligned.
This morning, for example, I had to merge several ASCII grids in order to form the DTM shown below:


The standard script will not work in this case, so I created a loop to use the tolerance option.
This is the whole script to use in with non-aligned grids:

 setwd("c:/temp")  
 library(rgdal)  
 library(raster)  
   
 # make a list of file names, perhaps like this:  
 f <-list.files(pattern = ".asc")  
   
 # turn these into a list of RasterLayer objects  
 r <- lapply(f, raster)  
   
   
 ##Approach to follow when the asc files are not aligned  
 for(i in 2:length(r)){  
 x<-merge(x=r[[1]],y=r[[i]],tolerance=5000,overlap=T)  
 r[[1]]<-x  
 }  
   
 #Write Ascii Grid  
 writeRaster(r[[1]],"DTM_10K_combine.asc")  

The loop merge the first ASCII grid with all the other iteratively, re-saving the first grid with the newly created one. This way I was able to create the DTM in the image above.

1 comment:

  1. Hi, nice blog and useful post. However, I don´t understand what the tolerance parameter means, and how to choose the best number (besides trial and error). Any hint would be awesome. R is not very helpful :/
    Thanks

    ReplyDelete

Note: only a member of this blog may post a comment.