1

I'm new on R.

I was listing .tif files from a folder to R. My final product is going to be a value mean of all the files in only one image. So I start stacking all of them, then brick the stack, and finally make a mean from the brick.

My problem starts when I stack the files, because the stack changes the original dimensions, extent and crs.

I need to maintain the values from the original files, so what can I do?

Example:

list160927<-list.files("D:/ImageR/160927", pattern=".tif", full.names=TRUE)
stack160927<-stack(list160927)
brick160927<-brick(stack160927)
mean160927<-mean(brick160927)

the original extent and csr from each file are:

dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 

when I stack all the files the

class       : RasterStack 
dimensions  : 1006, 804, 808824, 5  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 0, 804, 0, 1006  (xmin, xmax, ymin, ymax)
coord. ref. : NA 

Also, I tried to set the extent (the same that are in the original files) to the result brick but that didn't respond, and the brick maintain the change extent and csr.

  • try changing stack to raster::stack? I'm not sure what other packages use that command besides utils, but raster::stack shouldn't be capable of doing that. If one of those listed files isn't the same extent, it should fail to stack in the first place. Also: change your list.files pattern to "\\.tif$' so you don't accidentally list sidecar files like file.tif.aux.xml. – obrl_soil Apr 2 '17 at 0:17
  • I don't think your report is correct, can you include the code to get the extent and crs? You only need run stack(list160927); raster(list160927[1]) and share the output here. – mdsumner Apr 2 '17 at 1:50
1

Check the files you are using! Check that all your raster files have the same CRS and projection. I think the error should be there. I give you one very simple, reproducible and commented example below. Check this code:

# Load libraries

library("raster")
library("mapview") # (not required, but for nice interactive plots)

# Make some example data
# Same extent and projection but different cells values

r1 <- raster(matrix(data = rep(1, 100), nrow = 10, ncol = 10)) # raster1
extent(r1) <- c(-54, -52, -37, -35)
projection(r1) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

r2 <- raster(matrix(data = rep(2, 100), nrow = 10, ncol = 10)) # raster 2
extent(r2) <- c(-54, -52, -37, -35)
projection(r2) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

r3 <- raster(matrix(data = rep(3, 100), nrow = 10, ncol = 10)) # raster 3
extent(r3) <- c(-54, -52, -37, -35)
projection(r3) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

# Plot example data

mapview(r1, col.regions = "#EF2929") + mapview(r2, col.regions = "#AD7FA8") + mapview(r3, col.regions = "#729FCF") # stack

mapviewRasters

# Write rasters to working directory

writeRaster(r1, "r1.tif", overwrite = TRUE)
writeRaster(r2, "r2.tif", overwrite = TRUE)
writeRaster(r3, "r3.tif", overwrite = TRUE)

# List of files

listOfRasterFiles <- list.files(getwd(), pattern = "[r][0-9].tif", full.names = TRUE) # (look for a pattern that have "r", then one number between 0 and 9 and finally ".tif" string)

# One file properties

r1 <- raster("r1.tif")

# class       : RasterLayer 
# dimensions  : 10, 10, 100  (nrow, ncol, ncell)
# resolution  : 0.2, 0.2  (x, y)
# extent      : -54, -52, -37, -35  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
# data source : /home/glopez/r1.tif 
# names       : r1 
# values      : 1, 1  (min, max)

# Stack
rasterStack <- stack(listOfRasterFiles)

mapview(stack) # plot stack

# Stack properties

# class       : RasterStack 
# dimensions  : 10, 10, 100, 3  (nrow, ncol, ncell, nlayers)
# resolution  : 0.2, 0.2  (x, y)
# extent      : -54, -52, -37, -35  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
# names       : r1, r2, r3 
# min values  :  1,  2,  3 
# max values  :  1,  2,  3 

# Brick
rasterBrick <- brick(rasterStack)

# class       : RasterBrick 
# dimensions  : 10, 10, 100, 3  (nrow, ncol, ncell, nlayers)
# resolution  : 0.2, 0.2  (x, y)
# extent      : -54, -52, -37, -35  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
# data source : in memory
# names       : r1, r2, r3 
# min values  :  1,  2,  3 
# max values  :  1,  2,  3 

mapview(rasterBrick)

# Mean over layers
mean <- mean(rasterStack)

# class       : RasterLayer 
# dimensions  : 10, 10, 100  (nrow, ncol, ncell)
# resolution  : 0.2, 0.2  (x, y)
# extent      : -54, -52, -37, -35  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
# data source : in memory
# names       : layer 
# values      : 2, 2  (min, max)

Note: The final value of the mean function is 2 because is the mean of (1 + 2 + 3) / 3. That's because all the rasters have the same repeated value in all the cells: r1 have 100 cells with value 1, r2 have 100 cells with value 2, and so on.

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.