[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 17:38:17 CEST 2013
Gabrielle you are right, if the resolution is in fact different -- I
did not look carefully enough at the decimals -- then crop/extend
won't work and you will need resample.
Another approach would be to use the "terrain" function in raster and
compute slope and aspect from the elevation data.
Robert
On Sun, Sep 15, 2013 at 11:50 PM, Gabriele Cozzi <gab.cozzi at gmail.com> wrote:
> Dear Robert,
>
> thank you for going through my post again. Really appreciated
>
> I just tried what you suggested but neither crop() nor extend() do fix my
> problem; the extent and the resolution of forest remains different from
> slope and aspect
>
> Indeed aspect and slope have one more column than altitude, but I believe
> this is the result of the resolution being different in the original rasters
>
> res(aspect)
> res(slope)
> # 27.28532 27.28532
>
> res(altitude)
> # 27.28532 27.28677
>
> I believe this is the reason why cropping (or extending) does not work
> because a raster cannot be cropped using an extent that does not allow for a
> finite number of rows?
> I think that cropping aspect with altitude would result in a decimal number
> nrow(aspect) = xxx.y and a raster row can not take decimals?
>
> Does this make sense?
>
> Any further suggestion is appreciated
>
> I think I will also try to find out why, to begin with, the 3 raster layers
> that I was given have different resolution. They all come from the same
> satellite imagery so they should have the same resolution.
>
> Thanks in advance,
>
> Gabriele
>
>
>
> On Mon, Sep 16, 2013 at 8:07 AM, Robert J. Hijmans <r.hijmans at gmail.com>
> wrote:
>>
>> 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
>> >
>
>
>
>
> --
> 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
>
More information about the R-sig-Geo
mailing list