Geographic Information Systems Asked by Quechua on December 23, 2020
I have a list of rasters with the same projection and resolution but different origins and I need to merge them in a single raster file.
Here’s my code
l <- list.files('path', full.names=TRUE)
lst <- lapply(l, raster)
r_merged <- do.call(raster::merge, lst)
Error in compareRaster(x, extent = FALSE, rowcol = FALSE, orig = TRUE,
: different origin
Here’s what I’d like to get. Blue squares are original rasters, I want to get a single raster as the red square.
I would suggest that you use the stack
function from the raster
package. For instance, as the example in the documentation shows:
# file with one layer
fn <- system.file("external/test.grd", package="raster")
s <- stack(fn, fn)
r <- raster(fn)
s <- stack(r, fn)
nlayers(s)
Answered by dmci on December 23, 2020
I ran into a similar issue and developed this code to accomplish this task. The first function "reproject_align_raster" is called within "combine_rasters." combine_rasters takes a list of rasters and will align them all to a common grid and then merge them.
#' Reprojects/resamples and aligns a raster
#'
#' Reprojects/resamples and aligns a raster by matching a raster a raster to a specified origin, resolution, and coordinate reference system, or that of a reference raster. Useful for preparing adjacent areas before using raster::merge.
#' @param rast raster to be reprojected or resampled
#' @param ref_rast reference raster with desired properties (Alternatively can supply desired_origin, desired_res, and desired_crs)
#' @param desired_origin desired origin of output raster as a vector with length 2 (x,y)
#' @param desired_res desired resolution of output raster. Either an integer or a vector of length 2 (x,y)
#' @param desired_crs desired coordinate reference system of output raster (CRS class)
#' @param method resampling method. Either "bilinear" for bilinear interpolation (the default), or "ngb" for using the nearest neighbor
#' @importFrom raster crs
#' @importFrom raster extent
#' @importFrom raster origin
#' @importFrom raster projectExtent
#' @importFrom raster raster
#' @importFrom raster resample
#' @importFrom raster projectRaster
#' @export
reproject_align_raster<- function(rast, ref_rast=NULL, desired_origin, desired_res, desired_crs, method= "bilinear"){
if (!is.null(ref_rast)) {
desired_origin<- origin(ref_rast) #Desired origin
desired_res<- res(ref_rast) #Desired resolution
desired_crs<- crs(ref_rast) #Desired crs
} #Set parameters based on ref rast if it was supplied
if(length(desired_res)==1){
desired_res<- rep(desired_res,2)}
if(identical(crs(rast), desired_crs) & identical(origin(rast), desired_origin) & identical(desired_res, res(rast))){
message("raster was already aligned")
return(rast)} #Raster already aligned
if(identical(crs(rast), desired_crs)){
rast_orig_extent<- extent(rast)} else{
rast_orig_extent<- extent(projectExtent(object = rast, crs = desired_crs))} #reproject extent if crs is not the same
var1<- floor((rast_orig_extent@xmin - desired_origin[1])/desired_res[1])
new_xmin<-desired_origin[1]+ desired_res[1]*var1 #Calculate new minimum x value for extent
var2<- floor((rast_orig_extent@ymin - desired_origin[2])/desired_res[2])
new_ymin<-desired_origin[2]+ desired_res[2]*var2 #Calculate new minimum y value for extent
n_cols<- ceiling((rast_orig_extent@xmax-new_xmin)/desired_res[1]) #number of cols to be in output raster
n_rows<- ceiling((rast_orig_extent@ymax-new_ymin)/desired_res[2]) #number of rows to be in output raster
new_xmax<- new_xmin+(n_cols*desired_res[1]) #Calculate new max x value for extent
new_ymax<- new_ymin+(n_rows*desired_res[2]) #Calculate new max y value for extent
rast_new_template<- raster(xmn=new_xmin, xmx =new_xmax, ymn=new_ymin, ymx= new_ymax, res=desired_res, crs= desired_crs) #Create a blank template raster to fill with desired properties
if(!identical(desired_origin,origin(rast_new_template))){
message("desired origin does not match output origin")
stop()} #Throw error if origin doesn't match
if(identical(crs(rast),desired_crs)){
rast_new<- raster::resample(x=rast, y=rast_new_template, method = method)} else{
rast_new<- projectRaster(from=rast, to=rast_new_template, method = method)} #Use projectRaster if crs doesn't match and resample if they do
if(!identical(desired_origin,origin(rast_new))){
message("desired origin does not match output origin")
stop()} #Throw error if origin doesn't match
return(rast_new)
}
#' Combines multiple rasters into one
#'
#' Combines multiple rasters into one by using merge after first reprojecting or resampling and aligning rasters by matching them up with a a specified origin, resolution, and coordinate reference system, or that of a reference raster.
#' @param raster_list a list of rasters
#' @param rast raster to be reprojected or resampled
#' @param ref_rast reference raster with desired properties (Alternatively can supply desired_origin, desired_res, and desired_crs)
#' @param desired_origin desired origin of output raster as a vector with length 2 (x,y)
#' @param desired_res desired resolution of output raster. Either an integer or a vector of length 2 (x,y)
#' @param desired_crs desired coordinate reference system of output raster (CRS class)
#' @param method resampling method. Either "bilinear" for bilinear interpolation (the default), or "ngb" for using the nearest neighbor
#' @param display_progress Logical specifying whether or not to indicate progress
#' @importFrom raster merge
#' @export
combine_rasters<- function(raster_list, ref_rast=NULL, desired_origin, desired_res, desired_crs, method= "bilinear", display_progress=TRUE){
raster_list2<- vector("list", length = length(raster_list))
for (i in 1:length(raster_list)) {
if(display_progress){
print(paste("Reprojecting", as.character(i), "of", as.character(length(raster_list))))}
raster_list2[[i]]<- reproject_align_raster(raster_list[[i]], ref_rast=ref_rast, desired_origin=desired_origin, desired_res=desired_res, desired_crs=desired_crs, method= method)}
for (j in 1:length(raster_list2)) {
if(display_progress){
print(paste("Combining raster", as.character(j), "of", as.character(length(raster_list2))))}
if(j==1){
combined_raster<- raster_list2[[j]]} else{
combined_raster<- raster::merge(combined_raster, raster_list2[[j]])}}
return(combined_raster)
}
[1]: https://gis.stackexchange.com/questions/333464/re-projecting-a-raster-in-r-to-matching-resolution-crs-and-origin-but-differe/335605
Answered by ailich on December 23, 2020
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP