[R-sig-Geo] problems to stack() rasters due to different pixel resolution and raster extent

[Lionel Hertzog]

Dear Gabriele,

You should not use the 'res' function to set the raster resolution once 
it has values in it, to achieve this you need to use whether the 
'aggregate' or 'disaggregate' function to change the resolution since 
the cell values need to be re-computed.

In your case it appears that the difference in extent and resolution 
comes from one row (9813 vs 9812). I would recreate a new raster 
removing this extra row, below are some sample code:


#create two raster with one row difference leading to small variation in 
resolution and extent
e1<-extent(c(0,10,0,10))
r1<-raster(nrows=10,ncols=10,ext=e1)
values(r1)<-rnorm(100,0,1)
e2<-extent(c(0,10,0,11))
r2<-raster(nrows=11,ncols=10,ext=e2)
values(r2)<-rnorm(110,0,1)

#create a third raster having the dimension of r1 and taking the data 
from the first 10 rows of r2
r3<-r1
values(r3)<-r2[1:100]
stack(r1,r3)

There might be other more elegant ways to achieve this but this is how I 
would do it.

Yours,
Lionel


On 12/09/2013 12:59, Gabriele Cozzi wrote:
> Dear list,
>
> I was given three raster files with 'slope', 'aspect' and 'altitude' of my
> study site. I realised that one of the three rasters (namely 'altitude')
> has a slightly different raster cell resolution, which results in a
> slightly different raster extent and does not allow me to raster::stack()
> the three layers to create a single stack object.
> I tried to change the raster resolution using raster::res() for the raster
> 'altitude' but, when I subsequently cropped the three rasters with an
> object of class extent, I still got different extents between "slope and
> aspect" vs. "altitude".
>
> Any advice on how to solve this cell resolution incongruence is very
> welcome.
> Here below I report the output to make things more understandable.
>
> Thanks in advance,
> Gabriele
>
>
>
> #----------------
> #----------------
>
> aspect   <- raster("~/Envir/aspect/prj.adf")
> slope     <- raster("~/Envir/slope/prj.adf")
> altitude  <- raster("~/Envir/20130508094222_992463356.tif")
>
> aspect ( N.B. 'slope' returns exactly the same attributes!)
>     #class       : RasterLayer
>     #dimensions  : 9813, 8428, 82703964  (nrow, ncol, ncell)
>     #resolution  : 27.28532, 27.28532  (x, y)
>     #extent      : 135446.9, 365407.6, 4345919, 4613669  (xmin, xmax, ymin,
> ymax)
>     #coord. ref. : +proj=utm +zone=38 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0
> +units=m    +no_defs
>
> altitude
>     #class       : RasterLayer
>     #dimensions  : 9812, 8428, 82695536  (nrow, ncol, ncell)
>     #resolution  : 27.28532, 27.28677  (x, y)
>     #extent      : 135446.9, 365407.6, 4345919, 4613656  (xmin, xmax, ymin,
> ymax)
>     #coord. ref. : +proj=utm +zone=38 +datum=WGS84 +units=m +no_defs
> +ellps=WGS84 +towgs84=0,0,0
>
> ex <- extent(animal_locations)
> ex
>     #class       : Extent
>     #xmin        : 257473.1
>     #xmax        : 312499.1
>     #ymin        : 4452821
>     #ymax        : 4570717
>
> Before using the function  raster.crop <- crop('raster_file', ex) to crop
> all rasters to the same extent, I tried to change the resolution of
> 'altitude' using
>
> res(altitude) <- 27.28532
>
> altitude
>     #class       : RasterLayer
>     #dimensions  : 9813, 8428, 82703964  (nrow, ncol, ncell)
>     #resolution  : 27.28532, 27.28532  (x, y)
>     #extent      : 135446.9, 365407.6, 4345905, 4613656  (xmin, xmax, ymin,
> ymax)
>     #coord. ref. : +proj=utm +zone=38 +datum=WGS84 +units=m +no_defs
> +ellps=WGS84 +towgs84=0,0,0
>
> When I then run crop() for aspect and the new defined altitude I however
> become different extents (i.e. ymin and ymax) despite now the dimensions
> and resolutions are the same
>
> crop(aspect, ex)
>     #class       : RasterLayer
>     #dimensions  : 5054, 2749, 13893446  (nrow, ncol, ncell)
>     #resolution  : 27.28532, 27.28532  (x, y)
>     #extent      : 247480.5, 322487.8, 4442809, 4580709  (xmin, xmax, ymin,
> ymax)
>     #coord. ref. : +proj=utm +zone=38 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0
> +units=m +no_defs
>     #names       : prj
>     #values      : -1, 359.8158  (min, max)
>
> crop(altitude, ex)
>     #class       : RasterLayer
>     #dimensions  : 5054, 2749, 13893446  (nrow, ncol, ncell)
>     #resolution  : 27.28532, 27.28532  (x, y)
>     #extent      : 247480.5, 322487.8, 4442823, 4580723  (xmin, xmax, ymin,
> ymax)
>     #coord. ref. : +proj=utm +zone=38 +datum=WGS84 +units=m +no_defs
> +ellps=WGS84 +towgs84=0,0,0
>
> It also appears to me that by redefining res(altitude) I loose the actual
> values assigned to each pixel. Is there an argument to assign to the new
> pixels the value of the centre of the pixel?
>
>
>