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

Robert J. Hijmans r.hijmans at gmail.com
Mon Sep 16 08:07:56 CEST 2013


Gabriele,

That is not the best solution in this case; resampling should be a
function of last resort. As Lionel pointed out, aspect has one more
column than altitude, but has otherwise the same parameters.
Therefore, the approach should be using the crop function (which you
used, but after first making a mistake with res).

asp = crop(aspect, altitude)
slp = crop(slope, altitude)

OR the opposite route; add the row to the altitude raster:

alt = extend(altitude, aspect)

Robert

On Thu, Sep 12, 2013 at 8:45 AM, Gabriele Cozzi <gab.cozzi at gmail.com> wrote:
> Lionel, Paco,
>
> thanks for the quick answers.
>
> Going through your suggestions I found that the function *resample {raster}
> *that transfers values between non matching Raster* objects (in terms of
> origin and resolution) does exactly what I need.
>
> *resample* is actually embedded in the function *spatial_sync_raster* that
> I had problems to run. When I tried to run it I got the following error
> message:
>
> Error in spatial_sync_raster(alt, slo, method = "ngb", verbose = FALSE) :
>   could not find function "expand"
>
> *expand* is not part of the package climstats and I am not quite sure where
> to find it...
>
> But as said, *resample()* solved my problem.
>
> Thanks!
>
>
> Gabriele
>
>
>
> On Thu, Sep 12, 2013 at 1:33 PM, Francisco Rodriguez Sanchez <
> f.rodriguez.sanc at gmail.com> wrote:
>
>> Hi Gabriele,
>>
>> Try the spatial_sync_raster function in package climstats:
>> https://r-forge.r-project.org/**projects/climstats/<https://r-forge.r-project.org/projects/climstats/>.
>> It will align the projection, extent and resolution of your rasters all at
>> once
>>
>> Hope it helps,
>>
>> Paco
>>
>>
>> El 12/09/2013 11:59, Gabriele Cozzi escribió:
>>
>>  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?
>>>
>>>
>>>
>>>
>> --
>> Dr Francisco Rodriguez-Sanchez
>> Forest Ecology and Conservation Group
>> Department of Plant Sciences
>> University of Cambridge
>> Downing Street
>> Cambridge CB2 3EA
>> United Kingdom
>> http://sites.google.com/site/**rodriguezsanchezf<http://sites.google.com/site/rodriguezsanchezf>
>>
>>
>> ______________________________**_________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/**listinfo/r-sig-geo<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>>
>
>
>
> --
> Gabriele Cozzi
> Postdoctoral Research Associate
> Population Ecology Research Group
> http://www.popecol.org
>
> Zurich University
> Institute of Evolutionary Biology and Environmental Studies
> Winterthurerstr. 190
> 8057 Zurich - Switzerland
> E-mail: gabriele.cozzi at uzh.ch
> Phone: +41(0)44 635 49 12
> Fax: +41(0)16355711
> http://www.ieu.uzh.ch
>
>         [[alternative HTML version deleted]]
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>



More information about the R-sig-Geo mailing list