2

I have a RasterStack that I want to add two rasters to. I am using this code

library(rgdal)
library(raster)

#original rasterstack
rasterfiles   <- list.files("C:\\path_to_rasters", "*.tif$", full.names = TRUE)

#stack the bands
stacked <- raster::stack(rasterfiles)

#read in new rasters
aspect <- raster("C:\\Aspect.tif")
slope <- raster("C:\\Slope.tif")


#create a shapefile extent of the stacked raster
e <- extent(stacked)

# coerce to a SpatialPolygons object
p <- as(e, 'SpatialPolygons') 

#assign crs
projection(p) <- crs(stacked)

#function to mask the raster
mask_raster <- function(raster, shape){
      outr <- crop(raster, extent(shape))
      outr <- mask(outr, p)
      return(outr)
}

#apply function
slope <- mask_raster(raster= slope, shape = p)

aspect <- mask_raster(raster= aspect, shape = p)



#add new layers to stack
stacked <- addLayer(stacked, aspect)
stacked <- addLayer(stacked, slope)

This works completely fine for aspect but for slope when I try to add the layer I get the following error:

Error in compareRaster(c(x, rasters)) : different extent

The crs of all rasters are identical.

Before trying to mask this is what is in each object:

print(stacked) returns:

class       : RasterStack 
dimensions  : 5629, 253, 1424137, 428  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 480524, 480777, 5220075, 5225704  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=14 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names       :   Band_1,  Band_10, Band_100, Band_101, Band_102, Band_103, Band_104, Band_105, Band_106, Band_107, Band_108, Band_109,  Band_11, Band_110, Band_111, ... 
min values  :     0.00,    91.00,     0.00,     0.00,     0.00,     0.00,     0.00,     0.00,     0.00,     0.00,     0.00,     0.00,    80.00,     0.00,     0.00, ... 
max values  :  1914.00,  6061.00,  8805.00,  8830.00,  8745.00,  8455.00,  8230.00,  8317.00,  8092.00,  7945.00,  8124.00,  8122.00,  6160.00,  7642.00,  7222.00, ... 

print(aspect) returns:

class       : RasterLayer 
dimensions  : 21152, 22837, 483048224  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 475140, 497977, 5212857, 5234009  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=14 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : C:\NEON_2\lidar\Aspect.tif 
names       : Aspect 
values      : 0, 360  (min, max)

print(slope) returns:

class       : RasterLayer 
dimensions  : 23239, 29315, 681251285  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 471887, 501202, 5211828, 5235067  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=14 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : C:\NEON_2\lidar\Slope.tif 
names       : Slope 
values      : 0, 65.33745  (min, max)

After I mask each object contains this:

print(aspect) returns:

class       : RasterLayer 
dimensions  : 5629, 253, 1424137  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 480524, 480777, 5220075, 5225704  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=14 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : Aspect 
values      : 0, 360  (min, max)

print(slope) returns:

class       : RasterLayer 
dimensions  : 5629, 253, 1424137  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 480524, 480777, 5220075, 5225704  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=14 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : Slope 
values      : 0.0285432, 39.09481  (min, max)
  • Could you please edit your post to include the error? – Aaron Jun 5 '17 at 0:57
  • If i'm reading the code above correctly, your masking polygon will be the same size as the largest raster you have, so your masking function...extends the size of the raster? is that intentional? Might need to post the object summaries of both your inputs and outputs. – obrl_soil Jun 5 '17 at 1:52
  • the extent in stacked is smaller than slope or aspect I want to make slope and aspect fit to stacked so that I can add them as a band – Stefano Potter Jun 5 '17 at 2:08
  • 1
    Ok, per your GIS recomendation, the cells in slope are offset from those in aspect. aspect perfectly overlays stacked. This must be why, I will have to figure out a way to snap slope to the other rasters – Stefano Potter Jun 5 '17 at 4:48
  • 1
    resample is fine, but ngb is only appropriate for categorical data; use bilinear for these datasets. For aspect, be aware that either method can produce some wierdness since the data is circular. The better approach is to resample the source DEM and then recalculate aspect, or use trigonometry to convert aspect to a linear scale (northness or eastness) before resampling. – obrl_soil Jun 5 '17 at 10:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.