Search code examples
rrasterspatial

projectRaster fails to change crs when applied to a list object in R


I want to stack 6 rasters in a list called allrasters but first must fix crs and extent inconsistencies. Here is my code attempt to set the second raster in list to the crs of the third raster in list:

projectRaster(allrasters[[2]], crs=crs(allrasters[[3]]))

However when I run this code and check, allrasters[[2]] is still proj.merc and nothing has changed...

Raster information:

crs(allrasters[[2]])
CRS arguments:
 +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0
+x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext
+no_defs

crs(allrasters[[3]])
CRS arguments:
 +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5
+x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 


Solution

  • I think you need a few steps:

    1. you need to get all of your rasters in the same projection
    2. you need to find the full extent of all rasters as if they were mosaicked together
    3. you need to resample your rasters so that they have the same extent, resolution, and projection
    4. you will stack your rasters.

    Here is an example I created with some fake data that should help you accomplish this:

    
    ##Loading necessary packages##
    library(raster)
    library(rgeos)
    library(tmaptools)
    
    #For reproducibility#
    set.seed(52)
    
    ##Creating fake rasters with different extents and projections##
    R1<-raster(nrow=100, ncol=100, xmn=44.52, xmx=45.1, ymn=-122.1, ymx=-121.2, crs=crs("+init=epsg:4267"))
    R2<-raster(nrow=100, ncol=100, xmn=44.49, xmx=45.8, ymn=-122.0, ymx=-121.3, crs=crs("+init=epsg:4326"))
    R3<-raster(nrow=100, ncol=100, xmn=44.48, xmx=45.1, ymn=-122.5, ymx=-121.5, crs=crs("+init=epsg:4979"))
    R4<-raster(nrow=100, ncol=100, xmn=44.55, xmx=45.6, ymn=-122.2, ymx=-121.0, crs=crs("+init=epsg:4269"))
    
    values(R1)<-rnorm(10000, 500, 10)
    values(R2)<-rnorm(10000, 1000, 60)
    values(R3)<-rnorm(10000, 300, 10)
    values(R4)<-rnorm(10000, 2500, 70)
    
    ##Creating a list of the rasters##
    tmp<-list(R1,R2,R3,R4)
    
    ##Looping to reproject the rasters all into the same projection##
    allras<-list()
    for (i in 1:length(tmp)){
      if(i==1){
        allras[[i]]<-tmp[[i]]
      }else{
        allras[[i]]<-projectRaster(tmp[[i]], crs=crs(tmp[[1]]))
      }
      
    }
    
    ##Creating a function to make a polygon of each raster's extent##
    fxn<-function(ras){
      bb<-bbox(ras)
      bbpoly<-bb_poly(bb)
      st_crs(bbpoly)<-crs(ras)
      return(as_Spatial(bbpoly))
    }
    ext<-lapply(allras, fxn)
    
    ##Aggregating and dissolving all extents to get the full extent of all rasters##
    full.ext<-aggregate(do.call(bind, ext), dissolve=TRUE)
    
    ##Creating a blank raster with the full extent, the desired final projection, and the desired resolution##
    blank<-raster(ext=extent(full.ext), nrow=allras[[1]]@nrows, ncol=allras[[1]]@ncols, crs=allras[[1]]@crs)
    
    ##Resampling all rasters in the list to the desired extent and resolution##
    rastostack<-lapply(allras, resample, y=blank)
    
    ##Stacking the rasters##
    Ras<-stack(rastostack)