[GIS] Re-projecting a Raster in R to matching resolution, crs, and origin, but different extent

coordinate systemrrasterreprojection-mathematics

I have raster layers in a few different projections that I've loaded into R. I eventually plan on using the "merge" function of the raster package to bring these together as one surface, but to do that, all of the rasters need to have the same crs, resolution, and origin. Using projectRaster with "alignOnly=TRUE, says it will use "to or other parameters only to align the output (i.e. same origin and resolution), but use the projected extent from from," which sounds like what I need. However, this fails and gives me the warning "Warning in projectRaster(x, crs = projection(y)) : 'from' has no cell values." Does anyone have any suggestions on how to accomplish this? Is this a bug, or am I interpreting the documentation incorrectly?

Below is a small reproducible example of trying to reproject one raster with the GRS80 ellipsoid and 1.5m resolution to align with another raster of a different extent that uses the WGS84 ellipsoid and has 4m resolution.

#Create a raster with GRS80 crs and 1.5m resolution
r1<- raster(xmn=731458, xmx = 731458+100*1.5,  ymn=3116234, ymx= 3116234+100*1.5, 
nrow=100, ncol=100, vals=sample(1:50, size = 100*100, replace = TRUE), res=1.5, 
crs= "+proj=utm +zone=16 +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs ") 

#Create a raster with WGS84 crs and 4m resolution
r2<- raster(xmn=740000, xmx = 740000+100*4,  ymn=3000000, ymx= 3000000+100*4, 
nrow=100, ncol=100, vals=sample(1:50, size = 100*100, replace = TRUE), 
res=4, crs= "+proj=utm +zone=16 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs") 

#Reproject r1 to WGS with 4m resolution 
r1_res4_WGS84<- projectRaster(from = r1, crs = crs(r2), res= res(r2)) 

#Origins do not match after reprojecting    
origin(r1_res4_WGS84) #1.999999e+00 9.094551e-05
origin(r2) #0 0 

#Try reprojecting with alignOnly=TRUE        
r1_res4_WGS84_align<- projectRaster(from = r1, to = r2, alignOnly = TRUE) 
        # Warning in projectRaster(x, crs = projection(y)) :
        #   'from' has no cell values

Best Answer

After messing with it for a bit I was able to write up this function that seems to accomplish what I want.

#' 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
#' @import raster
#' @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)
}
Related Question