How to create a raster brick with rasters of different extents?
I am new in R so this question is very basic but I have been struggling with it and could not find a solution that worked. I want to create a raster brick from some landsat images of the same area.
They were downloaded in HDF-EOS format, and I used the Modis Reprojection Tool to convert them to .tif.
The resulting rasters have the same projection, but differ in their extent, resolution and origin.
I tried several approaches, summarized here below:
defining a subset extent manually and subsetting all the rasters. Then trying to make a brick with the subset rasters
Resampling the rasters, to give them the same number of columns and rows. Ideally that would ensure the raster cells are aligned and can be put into a raster brick. This option created a brick where rasters had no values, they were empty.
I am wondering what is the concept I should follow to correct the extent. Would it be correct (and efficient) to create an empty raster that I would fill in later with the values of the imported landsat image? Can you see where I am making a mistake? In case it is relevant, I am working on a Mac OSX Version 10.9.1, and using rgdal version 0.8-14
Any help will be very appreciated!
Thankyou
I add here the code I have been using:
# .tif files have been creating using the Modis Reprojection Tool. Input
# files used for this Tool was LANDSAT HDF-EOS imagery.
library(raster)
library(rgdal)
setwd()=getwd()
# Download the files from dropbox:
dl_from_dropbox <- function(x, key) {
require(RCurl)
bin <- getBinaryURL(paste0("https://dl.dropboxusercontent.com/s/", key, "/", x),
ssl.verifypeer = FALSE)
con <- file(x, open = "wb")
writeBin(bin, con)
close(con)
message(noquote(paste(x, "read into", getwd())))
}
dl_from_dropbox("lndsr.LT52210611985245CUB00-vi.NDVI.tif", "qb1bap9rghwivwy")
dl_from_dropbox("lndsr.LT52210611985309CUB00-vi.NDVI.tif", "sbhcffotirwnnc6")
dl_from_dropbox("lndsr.LT52210611987283CUB00-vi.NDVI.tif", "2zrkoo00ngigfzm")
# Create three rasters
tif1 <- "lndsr.LT52210611985245CUB00-vi.NDVI.tif"
tif2 <- "lndsr.LT52210611985309CUB00-vi.NDVI.tif"
tif3 <- "lndsr.LT52210611987283CUB00-vi.NDVI.tif"
r1 <- raster(tif1, values=TRUE)
r2 <- raster(tif2, band=1, values=TRUE)
r3 <- raster(tif3, band=1, values=TRUE)
### Display their properties
# projection is identical for the three rasters
projection(r1)
# "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
projection(r2)
projection(r3)
# Extents are different
extent(r1)
# class : Extent
# xmin : -45.85728
# xmax : -43.76855
# ymin : -2.388705
# ymax : -0.5181549
extent(r2)
# class : Extent
# xmin : -45.87077
# xmax : -43.78204
# ymin : -2.388727
# ymax : -0.5208711
extent(r3)
# class : Extent
# xmin : -45.81952
# xmax : -43.7173
# ymin : -2.405129
# ymax : -0.5154312
# origin differs for all
origin(r1)
# 5.644590e-05 -8.588605e-05
origin(r2)
# 0.0001122091 -0.0001045107
origin(r3)
# 6.949976e-05 -5.895945e-05
# resolution differs for r2
res(r1)
# 0.0002696872 0.0002696872
res(r2)
# 0.0002696875 0.0002696875
res(r3)
# 0.0002696872 0.0002696872
## Try different approaches to create a raster brick
# a- define a subset extent, and subset all the rasters
plot(r1, main="layer1 NDVI")
de <- drawExtent(show=TRUE, col="red")
de
# class : Extent
# xmin : -45.36159
# xmax : -45.30108
# ymin : -2.002435
# ymax : -1.949501
e <- extent(-45.36159,-45.30108,-2.002435,-1.949501)
# Crop each raster with this extent
r1c <- crop(r1,e)
r2c <- crop(r2,e)
r3c <- crop(r3,e)
# Make raster brick
rb_a <- brick(r1c,r2c,r3c)
# Error in compareRaster(x) : different extent
# b- Resample each raster
s <- raster(nrow=6926, ncol=7735) # smallest nrow and ncol among r1,r2 and r3
r1_res <- resample(r1,s, method="ngb")
r2_res <- resample(r2,s, method="ngb")
r3_res <- resample(r3,s, method="ngb")
# Resampling gives for the three rasters the following message:
# Warning message:
# In .local(x, y, ...) :
# you are resampling y a raster with a much larger cell size,
# perhaps you should use "aggregate" first
# Make raster brick
rb_c <- brick(r1, r2, r3)
# Error in compareRaster(x) : different extent
Answer
here are some things to help you out. Since I don't have your .tif files just some hints. Have you checked the extent on your raster s? It's the size of the world, with just those columns its cells are extremely large. So you have to add an extent to your raster before resampling it. From your info I did something like this:
# create an extent that includes all your data
e<-extent(-46, -43, -2, -0.6)
# create a raster with that extent, and the number of rows and colums to achive a
# similar resolution as you had before, you might have to do some math here....
# as crs, use the same crs as in your rasters before, from the crs slot
s<-raster(e, nrows=7000, ncols=7800, crs=r1@crs)
# use this raster to reproject your original raster (since your using the same crs,
# resample should work fine
r1<-resample(r1, s, method="ngb")
Happy Holidays, Ben
PS a better way to find your desired extent & resolution:
# dummy extent from your rasters, instead use lapply(raster list, extent)
a<-extent(-45.85728, -43.76855, -2.388705, -0.5181549)
b<-extent(-45.87077, -43.78204, -2.388727, -0.5208711)
c<-extent(-45.81952 ,-43.7173 , -2.405129 ,-0.5154312)
extent_list<-list(a, b, c)
# make a matrix out of it, each column represents a raster, rows the values
extent_list<-lapply(extent_list, as.matrix)
matrix_extent<-matrix(unlist(extent_list), ncol=length(extent_list)
rownames(matrix_extent)<-c("xmin", "ymin", "xmax", "ymax")
# create an extent with the extrem values of your extent
best_extent<-extent(min(matrix_extent[1,]), max(matrix_extent[3,]),
min(matrix_extent[2,]), max(matrix_extent[4,]))
# the range of your extent in degrees
ranges<-apply(as.matrix(best_extent), 1, diff)
# the resolution of your raster (pick one) or add a desired resolution
reso<-res(r1)
# deviding the range by your desired resolution gives you the number of rows and columns
nrow_ncol<-ranges/reso
# create your raster with the following
s<-raster(best_extent, nrows=nrow_ncol[2], ncols=nrow_ncol[1], crs=r1@crs)