Changing the raster spatial extent in R

user2807119 picture user2807119 · Jul 31, 2017 · Viewed 9k times · Source

I have two rasters and I want to make the spatial extent of one to another. Then save it as a new raster. I used following code. However, I cannot save the 2013 images with new spatial extent as a new raster. Any guidance is greatly appreciated.

raster_2013 <- raster("avgt2013.tif")
extent(raster_2013)
class       : Extent 
xmin        : 112.91 
xmax        : 153.64 
ymin        : -43.75 
ymax        : -9 
> res(raster_2013)
[1] 0.01 0.01
> 
> raster_2015 <- raster("avgt2015.tif")
> extent(raster_2015)
class       : Extent 
xmin        : 112 
xmax        : 154 
ymin        : -44 
ymax        : -9 
> res(raster_2015)
[1] 0.01 0.01
> 
> e <- extent(112, 154, -44, -9)
> 
> ex = extent(raster_2015)
> r2 = crop(raster_2013, ex)
> 
> 
> new_2013 <- alignExtent(e, raster_2013, snap='near')
> str(new_2013)
Formal class 'Extent' [package "raster"] with 4 slots
  ..@ xmin: num 112
  ..@ xmax: num 154
  ..@ ymin: num -44
  ..@ ymax: num -9
> 
> rc <- crop(raster_2013, e, snap='near')
> extent(rc)
class       : Extent 
xmin        : 112.91 
xmax        : 153.64 
ymin        : -43.75 
ymax        : -9 

Answer

ztl picture ztl · Aug 1, 2017

First, please make a simple reproducible example to ask a question.

library(raster)

set.seed(11)
raster_2013 = raster(ext=extent(112.91, 153.64, -43.75, -9), res=c(0.01, 0.01))
raster_2013[] = rnorm(ncell(raster_2013))
raster_2015 = raster(ext=extent(112, 154, -44, -9), res=c(0.01, 0.01))
raster_2015[] = rnorm(ncell(raster_2015))

Then, there are several issues with your code. In your case, alignExtent is useless since the two rasters have the same resolution and their extents correspond with regards to this resolution.

If your goal is to give the extent of raster_2015 to raster_2013, you need to realize that extent(raster_2015) is shorter (smaller) with respect to xmin, but larger or equal elsewhere. So cropping alone will just affect xmin of raster_2013. You first need to extend and second to crop in order to have the exact same extent:

new_2013 <- crop(extend(raster_2013, raster_2015), raster_2015)
all.equal(extent(raster_2015), extent(new_2013))
#[1] TRUE

As @Geo-sp mentions, you can also resample raster_2013, but you would typically use this if the rwo rasters are not aligned (and be aware that it would, in such case, result in modified data due to the interpolation). Here, since they are, it would give the same result as crop(extend()), but it would be much slower and more resource-consuming:

system.time(new_2013 <- crop(extend(raster_2013, raster_2015), raster_2015))
#  user  system elapsed 
# 0.676   0.036   0.712 
system.time(new_2013_res <- resample(raster_2013, raster_2015))
#   user  system elapsed 
# 10.324   0.536  10.869
all.equal(new_2013, new_2013_res)
# [1] TRUE

Finally, in order to save it, well... you can use writeRaster, as reading the documentation would have lead you to ;-)

writeRaster(new_2013, "raster_2013_extent2015.grd")