Grokbase Groups R r-sig-geo July 2013
FAQ

[R-sig-Geo] Quantiles (10th and 90th) of a Raster Stack?

Zia Uddin Ahmed
Jul 3, 2013 at 6:19 pm
I am trying to calculate 10- and 90-th quantiles of a raster stack . When I use following I always get quartile values from 0 to 100% with 25% interval. I am interested only 10th and 90-th quantiles of a set of raster. Any suggestion?
Thanks
Zia


library(raster)
fn <- system.file("external/test.grd", package="raster")
r <- raster(fn)
r2 <- raster(fn)+runif(ncell(r))
r3 <- raster(fn)+runif(ncell(r))
r4 <- raster(fn)+runif(ncell(r))
s <- stack(r, r2,r3, r4)


Q10 <- calc(s, fun = stats::quantile,probs = 0.10,na.rm=TRUE)
summary(Q10)
summary(Q10)
                0% 25% 50% 75% 100%
Min. 128.4340 128.4760 128.5845 128.7812 129.0879
1st Qu. 293.2325 293.4515 293.6626 293.9346 294.1769
Median 371.4120 371.6039 371.6643 371.6944 371.7912
3rd Qu. 499.8195 500.0833 500.3306 500.4950 500.5551
Max. 1805.7800 1805.8013 1805.8857 1806.0099 1806.1505
NA's 6097.0000 6097.0000 6097.0000 6097.0000 6097.0000
plot(Q10)


Q90 <- calc(s, fun = stats::quantile,probs = 0.90,na.rm=TRUE)
summary (Q90)
summary(Q90)
                0% 25% 50% 75% 100%
Min. 128.4340 128.4760 128.5845 128.7812 129.0879
1st Qu. 293.2325 293.4515 293.6626 293.9346 294.1769
Median 371.4120 371.6039 371.6643 371.6944 371.7912
3rd Qu. 499.8195 500.0833 500.3306 500.4950 500.5551
Max. 1805.7800 1805.8013 1805.8857 1806.0099 1806.1505
NA's 6097.0000 6097.0000 6097.0000 6097.0000 6097.0000


-----Original Message-----
From: r-sig-geo-bounces at r-project.org [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of Zia Uddin Ahmed
Sent: Wednesday, July 03, 2013 1:18 PM
To: Oscar Perpi??n Lamigueiro; Andrew Vitale
Cc: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] Help: Individual plots all raster in a raster stack in loop


Thank you so much.
  It works with following code:


saveGIF(
print(spplot(s, layout=c(1, 1))),
height = 500, width = 350, interval = .3, outdir = getwd())


-----Original Message-----
From: Oscar Perpi??n Lamigueiro [mailto:oscar.perpinan at gmail.com]
Sent: Wednesday, July 03, 2013 2:56 AM
To: Andrew Vitale; Zia Uddin Ahmed
Cc: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] Help: Individual plots all raster in a raster stack in loop


Hello,


Instead of using lapply, you can define layout=c(1, 1):


library(raster)


fn <- system.file("external/test.grd", package="raster") r <- raster(fn)
r2 <- raster(fn)+runif(ncell(r))
r3 <- raster(fn)+runif(ncell(r))
r4 <- raster(fn)+runif(ncell(r))
s <- stack(r, r2,r3, r4)


spplot(s, layout=c(1, 1))




On the other hand, you may be interested in this stackoverflow Q&A about animated plots:
http://stackoverflow.com/questions/1298100/creating-a-movie-from-a-series-of-plots-in-r


Last, there is a FAQ about the use of lattice/ggplot2 inside loops:
http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-do-lattice_002ftrellis-graphics-not-work_003f


Best,


Oscar.


Andrew Vitale <vit...@...com> writes:

I was about to give up when I looked at your original email again.
Looks like you just need to wrap it with print():

lapply(1:length(s), function(x){x11(); print(spplot(s[[x]], main > names(s[[x]])))})

Does anyone on the list know why?

On Tue, Jul 2, 2013 at 1:38 PM, Zia Uddin Ahmed wrote:

Hi Andrew,****

It works fine with plot function. But like to use ?spplot? function.
Because like to customize plot following ways to create animated gif plot.
****

However, when I use spplot function in your code, I have got plots
with NULL values. Thanks Zia****

** **
lapply(1:nlayers(s), function(x){spplot(s[[x]], main >> names(s[[x]]));x11()})****
[[1]]****

NULL****

** **

[[2]]****

NULL****

** **

[[3]]****

NULL****

** **

[[4]]****

NULL****

** **

# ****

bd<-readShapePoly("BD_boundary.shp")****

thana <- list("sp.lines", as(bd, "SpatialLines"), col="grey",
lwd=.7,lty=3)
****

at.NDVI <- seq(-3, +3, by = 0.5)****

rgb.palette <-
colorRampPalette(c("firebrick4","khaki1","green4"))****

** **

# Plot****

** **

spplot(s$May_2012, main = "",****

sp.layout=list(thana),at=at.NDVI,****

par.settings=list(axis.line=list(col="grey60",lwd=0.5)),****


colorkey=list(space="right",width=1.8,at=1:13,labels=list(cex=1.2,at>> 1:13,
****

labels=c("-3.0","","","-1.5","","", "0.0", "", "", "1.5", "", "",
"3.0"))),****

col.regions=rgb.palette(60))****

** **

# Animated GIF (This function did not work!)****

saveGIF(****

for (i in 1:ns){****

print(spplot(s[], main = list(label=paste("May_",i),cex=1.5),****

sp.layout=list(thana), at=at.NDVI,****

par.settings=list(axis.line=list(col="grey25",lwd=0.5)),****


colorkey=list(space="right",width=1.8,at=1:21,labels=list(cex=1.2,at>> 1:21,
****

labels=c("-3.0","","","-1.5","","", "0.0", "", "", "1.5", "", "",
"3.0"))),****

col.regions=rgb.palette(60)))}, ****

height = 500, width = 350, interval = .3, outdir = getwd())****

}****

** **

*From:* Andrew Vitale [mailto:vitale232 at gmail.com]
*Sent:* Tuesday, July 02, 2013 4:24 PM

*To:* Zia Uddin Ahmed
*Cc:* r-sig-geo at r-project.org
*Subject:* Re: [R-sig-Geo] Help: Individual plots all raster in a
raster stack in loop****

** **

Zia,****

** **

I like to use the plot function from the raster package.****

** **

?raster:::plot****

** **

To continue your example, I'm able to plot the individual layers of a
stack using the lapply() function. I'm not sure if there is a
simpler way to do it, but this way works for me:****

** **

library(raster)****

library(sp)****

library(gstat)****

** **

fn <- system.file("external/test.grd", package="raster")****

r <- raster(fn)****

r2 <- raster(fn)+runif(ncell(r))****

r3 <- raster(fn)+runif(ncell(r))****

r4 <- raster(fn)+runif(ncell(r))****

s <- stack(r, r2,r3, r4)****

** **

# Plot all rasters in a raster stack****

** **

lapply(1:nlayers(s), function(x){plot(s[[x]], main >> names(s[[x]]));x11()})
****

** **

** **

####And if you want the layers all in one plot, just use****

** **

plot(s)****

** **

** **

** **

** **

** **

** **

On Tue, Jul 2, 2013 at 10:47 AM, Zia Uddin Ahmed wrote:
****

I like to plot (individual plot all raster objects in raster stack in
a loop. Help will be appreciated.
Thanks
Zia

library(raster)
library(sp)
library(gstat)
fn <- system.file("external/test.grd", package="raster") r <-
raster(fn)
r2 <- raster(fn)+runif(ncell(r))
r3 <- raster(fn)+runif(ncell(r))
r4 <- raster(fn)+runif(ncell(r))
s <- stack(r, r2,r3, r4)
names(s)
# [1] "test.1" "test.2" "test.3" "test.4"

# Plot all rasters in a raster stack
print(spplot(s$test.1))
print(spplot(s$test.2))
print(spplot(s$test.3))
print(spplot(s$test.4))

From: Andrew Vitale [mailto:vitale232 at gmail.com]
Sent: Monday, July 01, 2013 3:13 PM
To: Zia Uddin Ahmed
Cc: Roman Lu?trik; r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] Error creating a standard deviation map from
raster stack

Did the sample code I sent not run on your system? That's strange.


On my machine:


library(raster)
fn <- system.file("external/test.grd", package="raster") r <-
raster(fn)
r2 <- raster(fn)+runif(ncell(r))
s <- stack(r, r2)
nlayers(s)
m<-mean(s)

std<-calc(s, fun = sd)

fun = function(x)
{ sqrt(var(x))
}
sd <- calc(s, fun)

par(mfrow=c(1,2));plot(sd,main = 'sd'); plot(std, main = 'std')

identical(sd, std)
#[1] TRUE
sessionInfo()
R version 3.0.1 (2013-05-16)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United
States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C [5]
LC_TIME=English_United States.1252

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] raster_2.1-37 sp_1.0-9

loaded via a namespace (and not attached):
[1] grid_3.0.1 lattice_0.20-15

On Mon, Jul 1, 2013 at 10:14 AM, Zia Uddin Ahmed <zua3 at cornell.edu<mailto:
zua3 at cornell.edu>> wrote:
Thanks you so much it works with following code. Zia

fun = function(x)
{ sqrt(var(x))
}
sd <- calc(s, fun)
plot(sd)

From: Andrew Vitale
[mailto:vitale232 at gmail.com<mailto:vitale232 at gmail.com
]
Sent: Monday, July 01, 2013 12:45 PM
To: Zia Uddin Ahmed
Cc: Roman Lu?trik;
r-sig-geo at r-project.org<mailto:r-sig-geo@r-project.org>

Subject: Re: [R-sig-Geo] Error creating a standard deviation map from
raster stack

Check out ?calc

library(raster)
fn <- system.file("external/test.grd", package="raster") r <-
raster(fn)
r2 <- raster(fn)+runif(ncell(r))
s <- stack(r, r2)
nlayers(s)
m<-mean(s)
std<-calc(s, fun = sd)


On Mon, Jul 1, 2013 at 9:30 AM, Zia Uddin Ahmed <zua3 at cornell.edu<mailto:
zua3 at cornell.edu>> wrote:
In my case, object ?s? is a set of 22 geo-tif files.
Thanks
zia

From: romunov at gmail.com[mailto:
romunov at gmail.com<mai...@...com>] On Behalf Of Roman
Lu?trik
Sent: Monday, July 01, 2013 12:09 PM
To: Zia Uddin Ahmed
Cc: r-sig-geo at r-project.org<mailto:r-sig-geo@r-project.org>
Subject: Re: [R-sig-Geo] Error creating a standard deviation map from
raster stack

What's the class/structure of object `s`?

Cheers,
Roman
On Mon, Jul 1, 2013 at 4:38 PM, Zia Uddin Ahmed <zua3 at cornell.edu<mailto:
zua3 at cornell.edu><mailto:zua3 at cornell.edu<mai...@...edu>>>
wrote:
I am trying to calculate standard deviation map from a raster stack.
I have got following error . Any suggestion?
std<-sd(s)
Error in as.double(x) :
cannot coerce type 'S4' to vector of type 'double'

Thanks
Zia

fn <- system.file("external/test.grd", package="raster") s <-
stack(fn, fn) r <- raster(fn) s <- stack(r, fn)
nlayers(s)
m<-mean(s)
std<-sd(s)
sessionInfo()
R version 3.0.1 (2013-05-16)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=English_United States.1252 [2] LC_CTYPE=English_United
States.1252 [3] LC_MONETARY=English_United States.1252 [4]
LC_NUMERIC=C [5] LC_TIME=English_United States.1252

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] gstat_1.0-16 rgdal_0.8-10 raster_2.1-37 sp_1.0-9

loaded via a namespace (and not attached):
[1] grid_3.0.1 intervals_0.14.0 lattice_0.20-15 spacetime_1.0-4
[5] xts_0.9-3 zoo_1.7-9
-----Original Message-----
From: r-sig-geo-bounces at r-project.org<mailto:
r-sig-geo-bounces at r-project.org><mailto:r-sig-geo-bounces at r-project.o
rg [mailto:
r-sig-geo-bounces at r-project.org<mailto:r-sig-geo-bounces at r-project.or
g
<mailto:r-sig-geo-bounces at r-project.org<mailto:
r-sig-geo-bounces at r-project.org>>] On Behalf Of ASANTOS
Sent: Sunday, June 23, 2013 4:08 PM
To: Sarah Goslee;
r-sig-geo at r-project.org<mailto:r-sig-geo at r-project.org
<mailto:r-sig-geo at r-project.org<mailto:r-sig-geo@r-project.org>>
Subject: Re: [R-sig-Geo] Problem with NDVI difference with subset
image

Hi Dra. Goslee,

I work in OS linux and I like very much your lssub
function, despites your notice. I think that the problems are in
linux system,
because:
landc<- stack(c("stackIm1.sample.tif","tackIm2.sample.tif"))
Error in .local(.Object, ...) :
`/home/asantos/Documentos/Sensoriamento remoto percevejo
bronzeador/Contornos/tackIm2.sample.tif' does not exist in the file
system, and is not recognised as a supported dataset name.

Error in .rasterObjectFromFile(x, band = band, objecttype >> "RasterLayer",
:
Cannot create a RasterLayer object from this file. (file does not
exist)


When I try to look the geotiff created
(stackIm1.sample.tif) there are something wrong, because the pictures
doesn't open (complete example below).

Thanks for your attention,

Alexandre

require(raster)
require(sp)
require(rgdal)
require(landsat)
#
#Create raster
r <- raster(nc00, nr00)
set.seed(20130622)
stackIm1 <- stack(lapply(1, function(x) setValues(r,
round(runif(ncell(r))* 255))))## Simulation red band
stackIm2 <- stack(lapply(1, function(x) setValues(r,
round(runif(ncell(r))* 255))))## Simulation nir

# define projection system
r.geo <- CRS("+proj=utm +zone# +south +datum=WGS84 +units=m
+no_defs") # geographical datum WGS84
proj4string(stackIm1) <- r.geo
proj4string(stackIm2) <- r.geo
#

#Create geotiff
writeRaster(stackIm1, filename="stackIm1.tif",
format="GTiff",overwrite=TRUE)
writeRaster(stackIm2, filename="stackIm2.tif",
format="GTiff",overwrite=TRUE)
#

#Subset a geotiff image 50 x 50 pixels
stackIm1.sample<-lssub("stackIm1.tif", "stackIm1.sample.tif", centerx
= 0, centery = 0, widthx = 50, widthy = 50)
stackIm2.sample<-lssub("stackIm1.tif", "sstackIm2.sample.tif",
centerx = 0, centery = 0, widthx = 50, widthy = 50) #
landc<- stack(c("stackIm1.sample.tif","tackIm2.sample.tif"))




Em 23/06/2013 12:59, Sarah Goslee escreveu:
Hi,

What happens?

Do you get an error? Where?

What is your sessionInfo()?

As it says in the lssub() help, this function was written for a
particular purpose, is only known to work on linux, and may not be
widely applicable.

It's a "use at your own risk" kind of function, and you may well be
better off using one of the many other methods available for
subsetting raster data, which would save you the whole "export as
GeoTIFF, subset, reimport" sequence.

Sarah

On Sun, Jun 23, 2013 at 11:36 AM, Alexandre Santos
<alexandresantosbr at yahoo.com.br<mailto:alexandresantosbr at yahoo.com.b
r
<mailto:alexandresantosbr at yahoo.com.br<mailto:
alexandresantosbr at yahoo.com.br>>> wrote:
Dear Members,


I'm having trouble calculating NDVI difference, first the images
Geotiff can not view using Ubuntu 4.12 64-bit and therefore can not
get to the NDVI, follow an example:
require(raster)
require(sp)
require(rgdal)
require(landsat)
#
#Create raster
r <- raster(nc00, nr00)
set.seed(20130622)
stackIm1 <- stack(lapply(1, function(x) setValues(r,
round(runif(ncell(r))* 255))))## Simulation red band
stackIm2 <- stack(lapply(1, function(x) setValues(r,
round(runif(ncell(r))* 255))))## Simulation nir
# define projection system
r.geo <- CRS("+proj=utm +zone# +south +datum=WGS84 +units=m
+no_defs") # geographical datum WGS84
proj4string(stackIm1) <- r.geo
proj4string(stackIm2) <- r.geo
#

#Create geotiff
writeRaster(stackIm1, filename="stackIm1.tif",
format="GTiff",overwrite=TRUE)
writeRaster(stackIm2, filename="stackIm2.tif",
format="GTiff",overwrite=TRUE)
#

#Subset a geotiff image 50 x 50 pixels
stackIm1.sample<-lssub("stackIm1.tif", "stackIm1.sample.tif",
centerx >> 0, centery = 0, widthx = 50, widthy = 50)
stackIm2.sample<-lssub("stackIm1.tif", "sstackIm2.sample.tif",
centerx
= 0, centery = 0, widthx = 50, widthy = 50)
#

#Calculate NDVI difference

multi.espc<-stack(c("stackIm1.sample.tif","stackIm2.sample.tif"))

band3<-raster(multi.espc,1)
band4<-raster(multi.espc,2)

ndvi<-(band4-band3)/(band4+band3)

#

--
====================================================================>> >> Alexandre dos Santos
Prote??o Florestal
Coordenador do curso T?cnico em Florestas Vice Coordenador do curso
de Engenharia Florestal IFMT - Instituto Federal de Educa??o, Ci?ncia
e Tecnologia de Mato Grosso Campus C?ceres Caixa Postal 244 Avenida
dos Ramires, s/n
Bairro: Distrito Industrial
C?ceres - MT CEP: 78.200-000
Fone: (+55) 65 8132-8112<tel:%28%2B55%29%2065%208132-8112<%28%2B55%29%2065%208132-8112>>
(TIM) (+55) 65 9686-6970<tel:%28%2B55%29%2065%209686-6970<%28%2B55%29%2065%209686-6970>>
(VIVO)
alexandre.santos at cas.ifmt.edu.br<mailto:
alexandre.santos at cas.ifmt.edu.br><mailto:alexandre.santos at cas.ifmt.ed
u.br <mai...@...br>>

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org<mailto:r-sig-geo@r-project.org><mailto:
R-sig-Geo at r-project.org<mailto:r-sig-geo@r-project.org>>
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org<mailto:r-sig-geo@r-project.org><mailto:
R-sig-Geo at r-project.org<mailto:r-sig-geo@r-project.org>>
https://stat.ethz.ch/mailman/listinfo/r-sig-geo



--
In God we trust, all others bring data.

[[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org<mailto:r-sig-geo@r-project.org>
https://stat.ethz.ch/mailman/listinfo/r-sig-geo



--
Andrew P. Vitale
Masters Student
Department of Geography
University of Nevada, Reno
(412) 915-3632<tel:%28412%29%20915-3632 <%28412%29%20915-3632>>
vitale232 at gmail.com<mailto:vit...@...com>



--
Andrew P. Vitale
Masters Student
Department of Geography
University of Nevada, Reno
(412) 915-3632
vitale232 at gmail.com<mailto:vit...@...com>

[[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****



****

** **

-- ****

*Andrew P. Vitale*****

Masters Student****

Department of Geography****

University of Nevada, Reno****

(412) 915-3632****

vitale232 at gmail.com****



--
Oscar Perpi??n Lamigueiro
Grupo de Sistemas Fotovoltaicos (IES-UPM) Dpto. Ingenier?a El?ctrica (EUITI-UPM)
URL: http://procomun.wordpress.com
Twitter: @oscarperpinan
LinkedIn: http://www.linkedin.com/in/oscarperpinan


_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
reply

Search Discussions

2 responses

  • Edzer Pebesma at Jul 3, 2013 at 6:31 pm
    Zia,


    The help of calc explains this:


         ...: Additional arguments as for ?writeRaster?


    meaning the additional agument (probs) are not passed on to fun.


    Defining the function more completely then passes probs to quantile:


    Q10 <- calc(s, fun = function(x) {quantile(x,probs = c(.1,.9),na.rm=TRUE)} )

    summary(Q10)
                   10% 90%
    Min. 128.5144 128.8897
    1st Qu. 293.2683 293.8457
    Median 371.4973 372.1300
    3rd Qu. 499.9098 500.5061
    Max. 1805.7937 1806.1727
    NA's 6097.0000 6097.0000





    On 07/03/2013 08:19 PM, Zia Uddin Ahmed wrote:
    I am trying to calculate 10- and 90-th quantiles of a raster stack . When I use following I always get quartile values from 0 to 100% with 25% interval. I am interested only 10th and 90-th quantiles of a set of raster. Any suggestion?
    Thanks
    Zia

    library(raster)
    fn <- system.file("external/test.grd", package="raster")
    r <- raster(fn)
    r2 <- raster(fn)+runif(ncell(r))
    r3 <- raster(fn)+runif(ncell(r))
    r4 <- raster(fn)+runif(ncell(r))
    s <- stack(r, r2,r3, r4)

    Q10 <- calc(s, fun = stats::quantile,probs = 0.10,na.rm=TRUE)
    summary(Q10)
    summary(Q10)
    0% 25% 50% 75% 100%
    Min. 128.4340 128.4760 128.5845 128.7812 129.0879
    1st Qu. 293.2325 293.4515 293.6626 293.9346 294.1769
    Median 371.4120 371.6039 371.6643 371.6944 371.7912
    3rd Qu. 499.8195 500.0833 500.3306 500.4950 500.5551
    Max. 1805.7800 1805.8013 1805.8857 1806.0099 1806.1505
    NA's 6097.0000 6097.0000 6097.0000 6097.0000 6097.0000
    plot(Q10)

    Q90 <- calc(s, fun = stats::quantile,probs = 0.90,na.rm=TRUE)
    summary (Q90)
    summary(Q90)
    0% 25% 50% 75% 100%
    Min. 128.4340 128.4760 128.5845 128.7812 129.0879
    1st Qu. 293.2325 293.4515 293.6626 293.9346 294.1769
    Median 371.4120 371.6039 371.6643 371.6944 371.7912
    3rd Qu. 499.8195 500.0833 500.3306 500.4950 500.5551
    Max. 1805.7800 1805.8013 1805.8857 1806.0099 1806.1505
    NA's 6097.0000 6097.0000 6097.0000 6097.0000 6097.0000

    -----Original Message-----
    From: r-sig-geo-bounces at r-project.org [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of Zia Uddin Ahmed
    Sent: Wednesday, July 03, 2013 1:18 PM
    To: Oscar Perpi??n Lamigueiro; Andrew Vitale
    Cc: r-sig-geo at r-project.org
    Subject: Re: [R-sig-Geo] Help: Individual plots all raster in a raster stack in loop

    Thank you so much.
    It works with following code:

    saveGIF(
    print(spplot(s, layout=c(1, 1))),
    height = 500, width = 350, interval = .3, outdir = getwd())

    -----Original Message-----
    From: Oscar Perpi??n Lamigueiro [mailto:oscar.perpinan at gmail.com]
    Sent: Wednesday, July 03, 2013 2:56 AM
    To: Andrew Vitale; Zia Uddin Ahmed
    Cc: r-sig-geo at r-project.org
    Subject: Re: [R-sig-Geo] Help: Individual plots all raster in a raster stack in loop

    Hello,

    Instead of using lapply, you can define layout=c(1, 1):

    library(raster)

    fn <- system.file("external/test.grd", package="raster") r <- raster(fn)
    r2 <- raster(fn)+runif(ncell(r))
    r3 <- raster(fn)+runif(ncell(r))
    r4 <- raster(fn)+runif(ncell(r))
    s <- stack(r, r2,r3, r4)

    spplot(s, layout=c(1, 1))


    On the other hand, you may be interested in this stackoverflow Q&A about animated plots:
    http://stackoverflow.com/questions/1298100/creating-a-movie-from-a-series-of-plots-in-r

    Last, there is a FAQ about the use of lattice/ggplot2 inside loops:
    http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-do-lattice_002ftrellis-graphics-not-work_003f

    Best,

    Oscar.

    Andrew Vitale <vit...@...com> writes:
    I was about to give up when I looked at your original email again.
    Looks like you just need to wrap it with print():

    lapply(1:length(s), function(x){x11(); print(spplot(s[[x]], main >> names(s[[x]])))})

    Does anyone on the list know why?

    On Tue, Jul 2, 2013 at 1:38 PM, Zia Uddin Ahmed wrote:

    Hi Andrew,****

    It works fine with plot function. But like to use ?spplot? function.
    Because like to customize plot following ways to create animated gif plot.
    ****

    However, when I use spplot function in your code, I have got plots
    with NULL values. Thanks Zia****

    ** **
    lapply(1:nlayers(s), function(x){spplot(s[[x]], main >>> names(s[[x]]));x11()})****
    [[1]]****

    NULL****

    ** **

    [[2]]****

    NULL****

    ** **

    [[3]]****

    NULL****

    ** **

    [[4]]****

    NULL****

    ** **

    # ****

    bd<-readShapePoly("BD_boundary.shp")****

    thana <- list("sp.lines", as(bd, "SpatialLines"), col="grey",
    lwd=.7,lty=3)
    ****

    at.NDVI <- seq(-3, +3, by = 0.5)****

    rgb.palette <-
    colorRampPalette(c("firebrick4","khaki1","green4"))****

    ** **

    # Plot****

    ** **

    spplot(s$May_2012, main = "",****

    sp.layout=list(thana),at=at.NDVI,****

    par.settings=list(axis.line=list(col="grey60",lwd=0.5)),****


    colorkey=list(space="right",width=1.8,at=1:13,labels=list(cex=1.2,at>>> 1:13,
    ****

    labels=c("-3.0","","","-1.5","","", "0.0", "", "", "1.5", "", "",
    "3.0"))),****

    col.regions=rgb.palette(60))****

    ** **

    # Animated GIF (This function did not work!)****

    saveGIF(****

    for (i in 1:ns){****

    print(spplot(s[], main = list(label=paste("May_",i),cex=1.5),****

    sp.layout=list(thana), at=at.NDVI,****

    par.settings=list(axis.line=list(col="grey25",lwd=0.5)),****


    colorkey=list(space="right",width=1.8,at=1:21,labels=list(cex=1.2,at>>> 1:21,
    ****

    labels=c("-3.0","","","-1.5","","", "0.0", "", "", "1.5", "", "",
    "3.0"))),****

    col.regions=rgb.palette(60)))}, ****

    height = 500, width = 350, interval = .3, outdir = getwd())****

    }****

    ** **

    *From:* Andrew Vitale [mailto:vitale232 at gmail.com]
    *Sent:* Tuesday, July 02, 2013 4:24 PM

    *To:* Zia Uddin Ahmed
    *Cc:* r-sig-geo at r-project.org
    *Subject:* Re: [R-sig-Geo] Help: Individual plots all raster in a
    raster stack in loop****

    ** **

    Zia,****

    ** **

    I like to use the plot function from the raster package.****

    ** **

    ?raster:::plot****

    ** **

    To continue your example, I'm able to plot the individual layers of a
    stack using the lapply() function. I'm not sure if there is a
    simpler way to do it, but this way works for me:****

    ** **

    library(raster)****

    library(sp)****

    library(gstat)****

    ** **

    fn <- system.file("external/test.grd", package="raster")****

    r <- raster(fn)****

    r2 <- raster(fn)+runif(ncell(r))****

    r3 <- raster(fn)+runif(ncell(r))****

    r4 <- raster(fn)+runif(ncell(r))****

    s <- stack(r, r2,r3, r4)****

    ** **

    # Plot all rasters in a raster stack****

    ** **

    lapply(1:nlayers(s), function(x){plot(s[[x]], main >>> names(s[[x]]));x11()})
    ****

    ** **

    ** **

    ####And if you want the layers all in one plot, just use****

    ** **

    plot(s)****

    ** **

    ** **

    ** **

    ** **

    ** **

    ** **

    On Tue, Jul 2, 2013 at 10:47 AM, Zia Uddin Ahmed wrote:
    ****

    I like to plot (individual plot all raster objects in raster stack in
    a loop. Help will be appreciated.
    Thanks
    Zia

    library(raster)
    library(sp)
    library(gstat)
    fn <- system.file("external/test.grd", package="raster") r <-
    raster(fn)
    r2 <- raster(fn)+runif(ncell(r))
    r3 <- raster(fn)+runif(ncell(r))
    r4 <- raster(fn)+runif(ncell(r))
    s <- stack(r, r2,r3, r4)
    names(s)
    # [1] "test.1" "test.2" "test.3" "test.4"

    # Plot all rasters in a raster stack
    print(spplot(s$test.1))
    print(spplot(s$test.2))
    print(spplot(s$test.3))
    print(spplot(s$test.4))

    From: Andrew Vitale [mailto:vitale232 at gmail.com]
    Sent: Monday, July 01, 2013 3:13 PM
    To: Zia Uddin Ahmed
    Cc: Roman Lu?trik; r-sig-geo at r-project.org
    Subject: Re: [R-sig-Geo] Error creating a standard deviation map from
    raster stack

    Did the sample code I sent not run on your system? That's strange.


    On my machine:


    library(raster)
    fn <- system.file("external/test.grd", package="raster") r <-
    raster(fn)
    r2 <- raster(fn)+runif(ncell(r))
    s <- stack(r, r2)
    nlayers(s)
    m<-mean(s)

    std<-calc(s, fun = sd)

    fun = function(x)
    { sqrt(var(x))
    }
    sd <- calc(s, fun)

    par(mfrow=c(1,2));plot(sd,main = 'sd'); plot(std, main = 'std')

    identical(sd, std)
    #[1] TRUE
    sessionInfo()
    R version 3.0.1 (2013-05-16)
    Platform: x86_64-w64-mingw32/x64 (64-bit)

    locale:
    [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United
    States.1252
    [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C [5]
    LC_TIME=English_United States.1252

    attached base packages:
    [1] stats graphics grDevices utils datasets methods base

    other attached packages:
    [1] raster_2.1-37 sp_1.0-9

    loaded via a namespace (and not attached):
    [1] grid_3.0.1 lattice_0.20-15

    On Mon, Jul 1, 2013 at 10:14 AM, Zia Uddin Ahmed <zua3 at cornell.edu<mailto:
    zua3 at cornell.edu>> wrote:
    Thanks you so much it works with following code. Zia

    fun = function(x)
    { sqrt(var(x))
    }
    sd <- calc(s, fun)
    plot(sd)

    From: Andrew Vitale
    [mailto:vitale232 at gmail.com<mailto:vitale232 at gmail.com
    ]
    Sent: Monday, July 01, 2013 12:45 PM
    To: Zia Uddin Ahmed
    Cc: Roman Lu?trik;
    r-sig-geo at r-project.org<mailto:r-sig-geo@r-project.org>

    Subject: Re: [R-sig-Geo] Error creating a standard deviation map from
    raster stack

    Check out ?calc

    library(raster)
    fn <- system.file("external/test.grd", package="raster") r <-
    raster(fn)
    r2 <- raster(fn)+runif(ncell(r))
    s <- stack(r, r2)
    nlayers(s)
    m<-mean(s)
    std<-calc(s, fun = sd)


    On Mon, Jul 1, 2013 at 9:30 AM, Zia Uddin Ahmed <zua3 at cornell.edu<mailto:
    zua3 at cornell.edu>> wrote:
    In my case, object ?s? is a set of 22 geo-tif files.
    Thanks
    zia

    From: romunov at gmail.com[mailto:
    romunov at gmail.com<mai...@...com>] On Behalf Of Roman
    Lu?trik
    Sent: Monday, July 01, 2013 12:09 PM
    To: Zia Uddin Ahmed
    Cc: r-sig-geo at r-project.org<mailto:r-sig-geo@r-project.org>
    Subject: Re: [R-sig-Geo] Error creating a standard deviation map from
    raster stack

    What's the class/structure of object `s`?

    Cheers,
    Roman
    On Mon, Jul 1, 2013 at 4:38 PM, Zia Uddin Ahmed <zua3 at cornell.edu<mailto:
    zua3 at cornell.edu><mailto:zua3 at cornell.edu<mai...@...edu>>>
    wrote:
    I am trying to calculate standard deviation map from a raster stack.
    I have got following error . Any suggestion?
    std<-sd(s)
    Error in as.double(x) :
    cannot coerce type 'S4' to vector of type 'double'

    Thanks
    Zia

    fn <- system.file("external/test.grd", package="raster") s <-
    stack(fn, fn) r <- raster(fn) s <- stack(r, fn)
    nlayers(s)
    m<-mean(s)
    std<-sd(s)
    sessionInfo()
    R version 3.0.1 (2013-05-16)
    Platform: x86_64-w64-mingw32/x64 (64-bit)

    locale:
    [1] LC_COLLATE=English_United States.1252 [2] LC_CTYPE=English_United
    States.1252 [3] LC_MONETARY=English_United States.1252 [4]
    LC_NUMERIC=C [5] LC_TIME=English_United States.1252

    attached base packages:
    [1] stats graphics grDevices utils datasets methods base

    other attached packages:
    [1] gstat_1.0-16 rgdal_0.8-10 raster_2.1-37 sp_1.0-9

    loaded via a namespace (and not attached):
    [1] grid_3.0.1 intervals_0.14.0 lattice_0.20-15 spacetime_1.0-4
    [5] xts_0.9-3 zoo_1.7-9
    -----Original Message-----
    From: r-sig-geo-bounces at r-project.org<mailto:
    r-sig-geo-bounces at r-project.org><mailto:r-sig-geo-bounces at r-project.o
    rg [mailto:
    r-sig-geo-bounces at r-project.org<mailto:r-sig-geo-bounces at r-project.or
    g
    <mailto:r-sig-geo-bounces at r-project.org<mailto:
    r-sig-geo-bounces at r-project.org>>] On Behalf Of ASANTOS
    Sent: Sunday, June 23, 2013 4:08 PM
    To: Sarah Goslee;
    r-sig-geo at r-project.org<mailto:r-sig-geo at r-project.org
    <mailto:r-sig-geo at r-project.org<mailto:r-sig-geo@r-project.org>>
    Subject: Re: [R-sig-Geo] Problem with NDVI difference with subset
    image

    Hi Dra. Goslee,

    I work in OS linux and I like very much your lssub
    function, despites your notice. I think that the problems are in
    linux system,
    because:
    landc<- stack(c("stackIm1.sample.tif","tackIm2.sample.tif"))
    Error in .local(.Object, ...) :
    `/home/asantos/Documentos/Sensoriamento remoto percevejo
    bronzeador/Contornos/tackIm2.sample.tif' does not exist in the file
    system, and is not recognised as a supported dataset name.

    Error in .rasterObjectFromFile(x, band = band, objecttype >>> "RasterLayer",
    :
    Cannot create a RasterLayer object from this file. (file does not
    exist)


    When I try to look the geotiff created
    (stackIm1.sample.tif) there are something wrong, because the pictures
    doesn't open (complete example below).

    Thanks for your attention,

    Alexandre

    require(raster)
    require(sp)
    require(rgdal)
    require(landsat)
    #
    #Create raster
    r <- raster(nc00, nr00)
    set.seed(20130622)
    stackIm1 <- stack(lapply(1, function(x) setValues(r,
    round(runif(ncell(r))* 255))))## Simulation red band
    stackIm2 <- stack(lapply(1, function(x) setValues(r,
    round(runif(ncell(r))* 255))))## Simulation nir

    # define projection system
    r.geo <- CRS("+proj=utm +zone# +south +datum=WGS84 +units=m
    +no_defs") # geographical datum WGS84
    proj4string(stackIm1) <- r.geo
    proj4string(stackIm2) <- r.geo
    #

    #Create geotiff
    writeRaster(stackIm1, filename="stackIm1.tif",
    format="GTiff",overwrite=TRUE)
    writeRaster(stackIm2, filename="stackIm2.tif",
    format="GTiff",overwrite=TRUE)
    #

    #Subset a geotiff image 50 x 50 pixels
    stackIm1.sample<-lssub("stackIm1.tif", "stackIm1.sample.tif", centerx
    = 0, centery = 0, widthx = 50, widthy = 50)
    stackIm2.sample<-lssub("stackIm1.tif", "sstackIm2.sample.tif",
    centerx = 0, centery = 0, widthx = 50, widthy = 50) #
    landc<- stack(c("stackIm1.sample.tif","tackIm2.sample.tif"))




    Em 23/06/2013 12:59, Sarah Goslee escreveu:
    Hi,

    What happens?

    Do you get an error? Where?

    What is your sessionInfo()?

    As it says in the lssub() help, this function was written for a
    particular purpose, is only known to work on linux, and may not be
    widely applicable.

    It's a "use at your own risk" kind of function, and you may well be
    better off using one of the many other methods available for
    subsetting raster data, which would save you the whole "export as
    GeoTIFF, subset, reimport" sequence.

    Sarah

    On Sun, Jun 23, 2013 at 11:36 AM, Alexandre Santos
    <alexandresantosbr at yahoo.com.br<mailto:alexandresantosbr at yahoo.com.b
    r
    <mailto:alexandresantosbr at yahoo.com.br<mailto:
    alexandresantosbr at yahoo.com.br>>> wrote:
    Dear Members,


    I'm having trouble calculating NDVI difference, first the images
    Geotiff can not view using Ubuntu 4.12 64-bit and therefore can not
    get to the NDVI, follow an example:
    require(raster)
    require(sp)
    require(rgdal)
    require(landsat)
    #
    #Create raster
    r <- raster(nc00, nr00)
    set.seed(20130622)
    stackIm1 <- stack(lapply(1, function(x) setValues(r,
    round(runif(ncell(r))* 255))))## Simulation red band
    stackIm2 <- stack(lapply(1, function(x) setValues(r,
    round(runif(ncell(r))* 255))))## Simulation nir
    # define projection system
    r.geo <- CRS("+proj=utm +zone# +south +datum=WGS84 +units=m
    +no_defs") # geographical datum WGS84
    proj4string(stackIm1) <- r.geo
    proj4string(stackIm2) <- r.geo
    #

    #Create geotiff
    writeRaster(stackIm1, filename="stackIm1.tif",
    format="GTiff",overwrite=TRUE)
    writeRaster(stackIm2, filename="stackIm2.tif",
    format="GTiff",overwrite=TRUE)
    #

    #Subset a geotiff image 50 x 50 pixels
    stackIm1.sample<-lssub("stackIm1.tif", "stackIm1.sample.tif",
    centerx >>> 0, centery = 0, widthx = 50, widthy = 50)
    stackIm2.sample<-lssub("stackIm1.tif", "sstackIm2.sample.tif",
    centerx
    = 0, centery = 0, widthx = 50, widthy = 50)
    #

    #Calculate NDVI difference

    multi.espc<-stack(c("stackIm1.sample.tif","stackIm2.sample.tif"))

    band3<-raster(multi.espc,1)
    band4<-raster(multi.espc,2)

    ndvi<-(band4-band3)/(band4+band3)

    #

    --
    ====================================================================>>> >>> Alexandre dos Santos
    Prote??o Florestal
    Coordenador do curso T?cnico em Florestas Vice Coordenador do curso
    de Engenharia Florestal IFMT - Instituto Federal de Educa??o, Ci?ncia
    e Tecnologia de Mato Grosso Campus C?ceres Caixa Postal 244 Avenida
    dos Ramires, s/n
    Bairro: Distrito Industrial
    C?ceres - MT CEP: 78.200-000
    Fone: (+55) 65 8132-8112<tel:%28%2B55%29%2065%208132-8112<%28%2B55%29%2065%208132-8112>>
    (TIM) (+55) 65 9686-6970<tel:%28%2B55%29%2065%209686-6970<%28%2B55%29%2065%209686-6970>>
    (VIVO)
    alexandre.santos at cas.ifmt.edu.br<mailto:
    alexandre.santos at cas.ifmt.edu.br><mailto:alexandre.santos at cas.ifmt.ed
    u.br <mai...@...br>>

    _______________________________________________
    R-sig-Geo mailing list
    R-sig-Geo at r-project.org<mailto:r-sig-geo@r-project.org><mailto:
    R-sig-Geo at r-project.org<mailto:r-sig-geo@r-project.org>>
    https://stat.ethz.ch/mailman/listinfo/r-sig-geo

    _______________________________________________
    R-sig-Geo mailing list
    R-sig-Geo at r-project.org<mailto:r-sig-geo@r-project.org><mailto:
    R-sig-Geo at r-project.org<mailto:r-sig-geo@r-project.org>>
    https://stat.ethz.ch/mailman/listinfo/r-sig-geo



    --
    In God we trust, all others bring data.

    [[alternative HTML version deleted]]

    _______________________________________________
    R-sig-Geo mailing list
    R-sig-Geo at r-project.org<mailto:r-sig-geo@r-project.org>
    https://stat.ethz.ch/mailman/listinfo/r-sig-geo



    --
    Andrew P. Vitale
    Masters Student
    Department of Geography
    University of Nevada, Reno
    (412) 915-3632<tel:%28412%29%20915-3632 <%28412%29%20915-3632>>
    vitale232 at gmail.com<mailto:vit...@...com>



    --
    Andrew P. Vitale
    Masters Student
    Department of Geography
    University of Nevada, Reno
    (412) 915-3632
    vitale232 at gmail.com<mailto:vit...@...com>

    [[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****



    ****

    ** **

    -- ****

    *Andrew P. Vitale*****

    Masters Student****

    Department of Geography****

    University of Nevada, Reno****

    (412) 915-3632****

    vitale232 at gmail.com****

    --
    Oscar Perpi??n Lamigueiro
    Grupo de Sistemas Fotovoltaicos (IES-UPM) Dpto. Ingenier?a El?ctrica (EUITI-UPM)
    URL: http://procomun.wordpress.com
    Twitter: @oscarperpinan
    LinkedIn: http://www.linkedin.com/in/oscarperpinan

    _______________________________________________
    R-sig-Geo mailing list
    R-sig-Geo at r-project.org
    https://stat.ethz.ch/mailman/listinfo/r-sig-geo
    _______________________________________________
    R-sig-Geo mailing list
    R-sig-Geo at r-project.org
    https://stat.ethz.ch/mailman/listinfo/r-sig-geo

    --
    Edzer Pebesma
    Institute for Geoinformatics (ifgi), University of M?nster
    Weseler Stra?e 253, 48151 M?nster, Germany. Phone: +49 251
    83 33081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de
  • Zia Uddin Ahmed at Jul 3, 2013 at 6:59 pm
    Thank you so much, I will try this.
    Zia


    -----Original Message-----
    From: r-sig-geo-bounces at r-project.org [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of Edzer Pebesma
    Sent: Wednesday, July 03, 2013 2:31 PM
    To: r-sig-geo at r-project.org
    Subject: Re: [R-sig-Geo] Quantiles (10th and 90th) of a Raster Stack?


    Zia,


    The help of calc explains this:


         ...: Additional arguments as for ?writeRaster?


    meaning the additional agument (probs) are not passed on to fun.


    Defining the function more completely then passes probs to quantile:


    Q10 <- calc(s, fun = function(x) {quantile(x,probs = c(.1,.9),na.rm=TRUE)} )

    summary(Q10)
                   10% 90%
    Min. 128.5144 128.8897
    1st Qu. 293.2683 293.8457
    Median 371.4973 372.1300
    3rd Qu. 499.9098 500.5061
    Max. 1805.7937 1806.1727
    NA's 6097.0000 6097.0000





    On 07/03/2013 08:19 PM, Zia Uddin Ahmed wrote:
    I am trying to calculate 10- and 90-th quantiles of a raster stack . When I use following I always get quartile values from 0 to 100% with 25% interval. I am interested only 10th and 90-th quantiles of a set of raster. Any suggestion?
    Thanks
    Zia

    library(raster)
    fn <- system.file("external/test.grd", package="raster") r <-
    raster(fn)
    r2 <- raster(fn)+runif(ncell(r))
    r3 <- raster(fn)+runif(ncell(r))
    r4 <- raster(fn)+runif(ncell(r))
    s <- stack(r, r2,r3, r4)

    Q10 <- calc(s, fun = stats::quantile,probs = 0.10,na.rm=TRUE)
    summary(Q10)
    summary(Q10)
    0% 25% 50% 75% 100%
    Min. 128.4340 128.4760 128.5845 128.7812 129.0879
    1st Qu. 293.2325 293.4515 293.6626 293.9346 294.1769
    Median 371.4120 371.6039 371.6643 371.6944 371.7912
    3rd Qu. 499.8195 500.0833 500.3306 500.4950 500.5551
    Max. 1805.7800 1805.8013 1805.8857 1806.0099 1806.1505
    NA's 6097.0000 6097.0000 6097.0000 6097.0000 6097.0000
    plot(Q10)

    Q90 <- calc(s, fun = stats::quantile,probs = 0.90,na.rm=TRUE) summary
    (Q90)
    summary(Q90)
    0% 25% 50% 75% 100%
    Min. 128.4340 128.4760 128.5845 128.7812 129.0879
    1st Qu. 293.2325 293.4515 293.6626 293.9346 294.1769
    Median 371.4120 371.6039 371.6643 371.6944 371.7912
    3rd Qu. 499.8195 500.0833 500.3306 500.4950 500.5551
    Max. 1805.7800 1805.8013 1805.8857 1806.0099 1806.1505
    NA's 6097.0000 6097.0000 6097.0000 6097.0000 6097.0000

    -----Original Message-----
    From: r-sig-geo-bounces at r-project.org
    [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of Zia Uddin Ahmed
    Sent: Wednesday, July 03, 2013 1:18 PM
    To: Oscar Perpi??n Lamigueiro; Andrew Vitale
    Cc: r-sig-geo at r-project.org
    Subject: Re: [R-sig-Geo] Help: Individual plots all raster in a raster
    stack in loop

    Thank you so much.
    It works with following code:

    saveGIF(
    print(spplot(s, layout=c(1, 1))),
    height = 500, width = 350, interval = .3, outdir = getwd())

    -----Original Message-----
    From: Oscar Perpi??n Lamigueiro [mailto:oscar.perpinan at gmail.com]
    Sent: Wednesday, July 03, 2013 2:56 AM
    To: Andrew Vitale; Zia Uddin Ahmed
    Cc: r-sig-geo at r-project.org
    Subject: Re: [R-sig-Geo] Help: Individual plots all raster in a raster
    stack in loop

    Hello,

    Instead of using lapply, you can define layout=c(1, 1):

    library(raster)

    fn <- system.file("external/test.grd", package="raster") r <-
    raster(fn)
    r2 <- raster(fn)+runif(ncell(r))
    r3 <- raster(fn)+runif(ncell(r))
    r4 <- raster(fn)+runif(ncell(r))
    s <- stack(r, r2,r3, r4)

    spplot(s, layout=c(1, 1))


    On the other hand, you may be interested in this stackoverflow Q&A about animated plots:
    http://stackoverflow.com/questions/1298100/creating-a-movie-from-a-ser
    ies-of-plots-in-r

    Last, there is a FAQ about the use of lattice/ggplot2 inside loops:
    http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-do-lattice_002ftrelli
    s-graphics-not-work_003f

    Best,

    Oscar.

    Andrew Vitale <vit...@...com> writes:
    I was about to give up when I looked at your original email again.
    Looks like you just need to wrap it with print():

    lapply(1:length(s), function(x){x11(); print(spplot(s[[x]], main >> names(s[[x]])))})

    Does anyone on the list know why?

    On Tue, Jul 2, 2013 at 1:38 PM, Zia Uddin Ahmed wrote:

    Hi Andrew,****

    It works fine with plot function. But like to use ?spplot? function.
    Because like to customize plot following ways to create animated gif plot.
    ****

    However, when I use spplot function in your code, I have got plots
    with NULL values. Thanks Zia****

    ** **
    lapply(1:nlayers(s), function(x){spplot(s[[x]], main >>> names(s[[x]]));x11()})****
    [[1]]****

    NULL****

    ** **

    [[2]]****

    NULL****

    ** **

    [[3]]****

    NULL****

    ** **

    [[4]]****

    NULL****

    ** **

    # ****

    bd<-readShapePoly("BD_boundary.shp")****

    thana <- list("sp.lines", as(bd, "SpatialLines"), col="grey",
    lwd=.7,lty=3)
    ****

    at.NDVI <- seq(-3, +3, by = 0.5)****

    rgb.palette <-
    colorRampPalette(c("firebrick4","khaki1","green4"))****

    ** **

    # Plot****

    ** **

    spplot(s$May_2012, main = "",****

    sp.layout=list(thana),at=at.NDVI,****

    par.settings=list(axis.line=list(col="grey60",lwd=0.5)),****


    colorkey=list(space="right",width=1.8,at=1:13,labels=list(cex=1.2,at
    1:13,
    ****

    labels=c("-3.0","","","-1.5","","", "0.0", "", "", "1.5", "", "",
    "3.0"))),****

    col.regions=rgb.palette(60))****

    ** **

    # Animated GIF (This function did not work!)****

    saveGIF(****

    for (i in 1:ns){****

    print(spplot(s[], main = list(label=paste("May_",i),cex=1.5),****

    sp.layout=list(thana), at=at.NDVI,****

    par.settings=list(axis.line=list(col="grey25",lwd=0.5)),****


    colorkey=list(space="right",width=1.8,at=1:21,labels=list(cex=1.2,at
    1:21,
    ****

    labels=c("-3.0","","","-1.5","","", "0.0", "", "", "1.5", "", "",
    "3.0"))),****

    col.regions=rgb.palette(60)))}, ****

    height = 500, width = 350, interval = .3, outdir = getwd())****

    }****

    ** **

    *From:* Andrew Vitale [mailto:vitale232 at gmail.com]
    *Sent:* Tuesday, July 02, 2013 4:24 PM

    *To:* Zia Uddin Ahmed
    *Cc:* r-sig-geo at r-project.org
    *Subject:* Re: [R-sig-Geo] Help: Individual plots all raster in a
    raster stack in loop****

    ** **

    Zia,****

    ** **

    I like to use the plot function from the raster package.****

    ** **

    ?raster:::plot****

    ** **

    To continue your example, I'm able to plot the individual layers of
    a stack using the lapply() function. I'm not sure if there is a
    simpler way to do it, but this way works for me:****

    ** **

    library(raster)****

    library(sp)****

    library(gstat)****

    ** **

    fn <- system.file("external/test.grd", package="raster")****

    r <- raster(fn)****

    r2 <- raster(fn)+runif(ncell(r))****

    r3 <- raster(fn)+runif(ncell(r))****

    r4 <- raster(fn)+runif(ncell(r))****

    s <- stack(r, r2,r3, r4)****

    ** **

    # Plot all rasters in a raster stack****

    ** **

    lapply(1:nlayers(s), function(x){plot(s[[x]], main >>> names(s[[x]]));x11()})
    ****

    ** **

    ** **

    ####And if you want the layers all in one plot, just use****

    ** **

    plot(s)****

    ** **

    ** **

    ** **

    ** **

    ** **

    ** **

    On Tue, Jul 2, 2013 at 10:47 AM, Zia Uddin Ahmed wrote:
    ****

    I like to plot (individual plot all raster objects in raster stack
    in a loop. Help will be appreciated.
    Thanks
    Zia

    library(raster)
    library(sp)
    library(gstat)
    fn <- system.file("external/test.grd", package="raster") r <-
    raster(fn)
    r2 <- raster(fn)+runif(ncell(r))
    r3 <- raster(fn)+runif(ncell(r))
    r4 <- raster(fn)+runif(ncell(r))
    s <- stack(r, r2,r3, r4)
    names(s)
    # [1] "test.1" "test.2" "test.3" "test.4"

    # Plot all rasters in a raster stack
    print(spplot(s$test.1))
    print(spplot(s$test.2))
    print(spplot(s$test.3))
    print(spplot(s$test.4))

    From: Andrew Vitale [mailto:vitale232 at gmail.com]
    Sent: Monday, July 01, 2013 3:13 PM
    To: Zia Uddin Ahmed
    Cc: Roman Lu?trik; r-sig-geo at r-project.org
    Subject: Re: [R-sig-Geo] Error creating a standard deviation map
    from raster stack

    Did the sample code I sent not run on your system? That's strange.


    On my machine:


    library(raster)
    fn <- system.file("external/test.grd", package="raster") r <-
    raster(fn)
    r2 <- raster(fn)+runif(ncell(r))
    s <- stack(r, r2)
    nlayers(s)
    m<-mean(s)

    std<-calc(s, fun = sd)

    fun = function(x)
    { sqrt(var(x))
    }
    sd <- calc(s, fun)

    par(mfrow=c(1,2));plot(sd,main = 'sd'); plot(std, main = 'std')

    identical(sd, std)
    #[1] TRUE
    sessionInfo()
    R version 3.0.1 (2013-05-16)
    Platform: x86_64-w64-mingw32/x64 (64-bit)

    locale:
    [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United
    States.1252
    [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C [5]
    LC_TIME=English_United States.1252

    attached base packages:
    [1] stats graphics grDevices utils datasets methods base

    other attached packages:
    [1] raster_2.1-37 sp_1.0-9

    loaded via a namespace (and not attached):
    [1] grid_3.0.1 lattice_0.20-15

    On Mon, Jul 1, 2013 at 10:14 AM, Zia Uddin Ahmed <zua3 at cornell.edu<mailto:
    zua3 at cornell.edu>> wrote:
    Thanks you so much it works with following code. Zia

    fun = function(x)
    { sqrt(var(x))
    }
    sd <- calc(s, fun)
    plot(sd)

    From: Andrew Vitale
    [mailto:vitale232 at gmail.com<mailto:vitale232 at gmail.com
    ]
    Sent: Monday, July 01, 2013 12:45 PM
    To: Zia Uddin Ahmed
    Cc: Roman Lu?trik;
    r-sig-geo at r-project.org<mailto:r-sig-geo@r-project.org>

    Subject: Re: [R-sig-Geo] Error creating a standard deviation map
    from raster stack

    Check out ?calc

    library(raster)
    fn <- system.file("external/test.grd", package="raster") r <-
    raster(fn)
    r2 <- raster(fn)+runif(ncell(r))
    s <- stack(r, r2)
    nlayers(s)
    m<-mean(s)
    std<-calc(s, fun = sd)


    On Mon, Jul 1, 2013 at 9:30 AM, Zia Uddin Ahmed <zua3 at cornell.edu<mailto:
    zua3 at cornell.edu>> wrote:
    In my case, object ?s? is a set of 22 geo-tif files.
    Thanks
    zia

    From: romunov at gmail.com[mailto:
    romunov at gmail.com<mai...@...com>] On Behalf Of Roman
    Lu?trik
    Sent: Monday, July 01, 2013 12:09 PM
    To: Zia Uddin Ahmed
    Cc: r-sig-geo at r-project.org<mailto:r-sig-geo@r-project.org>
    Subject: Re: [R-sig-Geo] Error creating a standard deviation map
    from raster stack

    What's the class/structure of object `s`?

    Cheers,
    Roman
    On Mon, Jul 1, 2013 at 4:38 PM, Zia Uddin Ahmed <zua3 at cornell.edu<mailto:
    zua3 at cornell.edu><mailto:zua3 at cornell.edu<mai...@...edu>>>
    wrote:
    I am trying to calculate standard deviation map from a raster stack.
    I have got following error . Any suggestion?
    std<-sd(s)
    Error in as.double(x) :
    cannot coerce type 'S4' to vector of type 'double'

    Thanks
    Zia

    fn <- system.file("external/test.grd", package="raster") s <-
    stack(fn, fn) r <- raster(fn) s <- stack(r, fn)
    nlayers(s)
    m<-mean(s)
    std<-sd(s)
    sessionInfo()
    R version 3.0.1 (2013-05-16)
    Platform: x86_64-w64-mingw32/x64 (64-bit)

    locale:
    [1] LC_COLLATE=English_United States.1252 [2]
    LC_CTYPE=English_United
    States.1252 [3] LC_MONETARY=English_United States.1252 [4]
    LC_NUMERIC=C [5] LC_TIME=English_United States.1252

    attached base packages:
    [1] stats graphics grDevices utils datasets methods base

    other attached packages:
    [1] gstat_1.0-16 rgdal_0.8-10 raster_2.1-37 sp_1.0-9

    loaded via a namespace (and not attached):
    [1] grid_3.0.1 intervals_0.14.0 lattice_0.20-15 spacetime_1.0-4
    [5] xts_0.9-3 zoo_1.7-9
    -----Original Message-----
    From: r-sig-geo-bounces at r-project.org<mailto:
    r-sig-geo-bounces at r-project.org><mailto:r-sig-geo-bounces at r-project.
    o rg [mailto:
    r-sig-geo-bounces at r-project.org<mailto:r-sig-geo-bounces at r-project.o
    r
    g
    <mailto:r-sig-geo-bounces at r-project.org<mailto:
    r-sig-geo-bounces at r-project.org>>] On Behalf Of ASANTOS
    Sent: Sunday, June 23, 2013 4:08 PM
    To: Sarah Goslee;
    r-sig-geo at r-project.org<mailto:r-sig-geo at r-project.org
    <mailto:r-sig-geo at r-project.org<mailto:r-sig-geo@r-project.org>>
    Subject: Re: [R-sig-Geo] Problem with NDVI difference with subset
    image

    Hi Dra. Goslee,

    I work in OS linux and I like very much your lssub
    function, despites your notice. I think that the problems are in
    linux system,
    because:
    landc<- stack(c("stackIm1.sample.tif","tackIm2.sample.tif"))
    Error in .local(.Object, ...) :
    `/home/asantos/Documentos/Sensoriamento remoto percevejo
    bronzeador/Contornos/tackIm2.sample.tif' does not exist in the file
    system, and is not recognised as a supported dataset name.

    Error in .rasterObjectFromFile(x, band = band, objecttype >>> "RasterLayer",
    :
    Cannot create a RasterLayer object from this file. (file does not
    exist)


    When I try to look the geotiff created
    (stackIm1.sample.tif) there are something wrong, because the
    pictures doesn't open (complete example below).

    Thanks for your attention,

    Alexandre

    require(raster)
    require(sp)
    require(rgdal)
    require(landsat)
    #
    #Create raster
    r <- raster(nc00, nr00)
    set.seed(20130622)
    stackIm1 <- stack(lapply(1, function(x) setValues(r,
    round(runif(ncell(r))* 255))))## Simulation red band
    stackIm2 <- stack(lapply(1, function(x) setValues(r,
    round(runif(ncell(r))* 255))))## Simulation nir

    # define projection system
    r.geo <- CRS("+proj=utm +zone# +south +datum=WGS84 +units=m
    +no_defs") # geographical datum WGS84
    proj4string(stackIm1) <- r.geo
    proj4string(stackIm2) <- r.geo
    #

    #Create geotiff
    writeRaster(stackIm1, filename="stackIm1.tif",
    format="GTiff",overwrite=TRUE)
    writeRaster(stackIm2, filename="stackIm2.tif",
    format="GTiff",overwrite=TRUE)
    #

    #Subset a geotiff image 50 x 50 pixels
    stackIm1.sample<-lssub("stackIm1.tif", "stackIm1.sample.tif",
    centerx = 0, centery = 0, widthx = 50, widthy = 50)
    stackIm2.sample<-lssub("stackIm1.tif", "sstackIm2.sample.tif",
    centerx = 0, centery = 0, widthx = 50, widthy = 50) #
    landc<- stack(c("stackIm1.sample.tif","tackIm2.sample.tif"))




    Em 23/06/2013 12:59, Sarah Goslee escreveu:
    Hi,

    What happens?

    Do you get an error? Where?

    What is your sessionInfo()?

    As it says in the lssub() help, this function was written for a
    particular purpose, is only known to work on linux, and may not be
    widely applicable.

    It's a "use at your own risk" kind of function, and you may well be
    better off using one of the many other methods available for
    subsetting raster data, which would save you the whole "export as
    GeoTIFF, subset, reimport" sequence.

    Sarah

    On Sun, Jun 23, 2013 at 11:36 AM, Alexandre Santos
    <alexandresantosbr at yahoo.com.br<mailto:alexandresantosbr at yahoo.com.
    b
    r
    <mailto:alexandresantosbr at yahoo.com.br<mailto:
    alexandresantosbr at yahoo.com.br>>> wrote:
    Dear Members,


    I'm having trouble calculating NDVI difference, first the images
    Geotiff can not view using Ubuntu 4.12 64-bit and therefore can not
    get to the NDVI, follow an example:
    require(raster)
    require(sp)
    require(rgdal)
    require(landsat)
    #
    #Create raster
    r <- raster(nc00, nr00)
    set.seed(20130622)
    stackIm1 <- stack(lapply(1, function(x) setValues(r,
    round(runif(ncell(r))* 255))))## Simulation red band
    stackIm2 <- stack(lapply(1, function(x) setValues(r,
    round(runif(ncell(r))* 255))))## Simulation nir
    # define projection system
    r.geo <- CRS("+proj=utm +zone# +south +datum=WGS84 +units=m
    +no_defs") # geographical datum WGS84
    proj4string(stackIm1) <- r.geo
    proj4string(stackIm2) <- r.geo
    #

    #Create geotiff
    writeRaster(stackIm1, filename="stackIm1.tif",
    format="GTiff",overwrite=TRUE)
    writeRaster(stackIm2, filename="stackIm2.tif",
    format="GTiff",overwrite=TRUE)
    #

    #Subset a geotiff image 50 x 50 pixels
    stackIm1.sample<-lssub("stackIm1.tif", "stackIm1.sample.tif",
    centerx >>> 0, centery = 0, widthx = 50, widthy = 50)
    stackIm2.sample<-lssub("stackIm1.tif", "sstackIm2.sample.tif",
    centerx
    = 0, centery = 0, widthx = 50, widthy = 50)
    #

    #Calculate NDVI difference

    multi.espc<-stack(c("stackIm1.sample.tif","stackIm2.sample.tif"))

    band3<-raster(multi.espc,1)
    band4<-raster(multi.espc,2)

    ndvi<-(band4-band3)/(band4+band3)

    #

    --
    ===================================================================>>> >>> >>> Alexandre dos Santos
    Prote??o Florestal
    Coordenador do curso T?cnico em Florestas Vice Coordenador do curso
    de Engenharia Florestal IFMT - Instituto Federal de Educa??o,
    Ci?ncia e Tecnologia de Mato Grosso Campus C?ceres Caixa Postal 244
    Avenida dos Ramires, s/n
    Bairro: Distrito Industrial
    C?ceres - MT CEP: 78.200-000
    Fone: (+55) 65 8132-8112<tel:%28%2B55%29%2065%208132-8112<%28%2B55%29%2065%208132-8112>>
    (TIM) (+55) 65 9686-6970<tel:%28%2B55%29%2065%209686-6970<%28%2B55%29%2065%209686-6970>>
    (VIVO)
    alexandre.santos at cas.ifmt.edu.br<mailto:
    alexandre.santos at cas.ifmt.edu.br><mailto:alexandre.santos at cas.ifmt.e
    d u.br <mai...@...br>>

    _______________________________________________
    R-sig-Geo mailing list
    R-sig-Geo at r-project.org<mailto:r-sig-geo@r-project.org><mailto:
    R-sig-Geo at r-project.org<mailto:r-sig-geo@r-project.org>>
    https://stat.ethz.ch/mailman/listinfo/r-sig-geo

    _______________________________________________
    R-sig-Geo mailing list
    R-sig-Geo at r-project.org<mailto:r-sig-geo@r-project.org><mailto:
    R-sig-Geo at r-project.org<mailto:r-sig-geo@r-project.org>>
    https://stat.ethz.ch/mailman/listinfo/r-sig-geo



    --
    In God we trust, all others bring data.

    [[alternative HTML version deleted]]

    _______________________________________________
    R-sig-Geo mailing list
    R-sig-Geo at r-project.org<mailto:r-sig-geo@r-project.org>
    https://stat.ethz.ch/mailman/listinfo/r-sig-geo



    --
    Andrew P. Vitale
    Masters Student
    Department of Geography
    University of Nevada, Reno
    (412) 915-3632<tel:%28412%29%20915-3632 <%28412%29%20915-3632>>
    vitale232 at gmail.com<mailto:vit...@...com>



    --
    Andrew P. Vitale
    Masters Student
    Department of Geography
    University of Nevada, Reno
    (412) 915-3632
    vitale232 at gmail.com<mailto:vit...@...com>

    [[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****



    ****

    ** **

    -- ****

    *Andrew P. Vitale*****

    Masters Student****

    Department of Geography****

    University of Nevada, Reno****

    (412) 915-3632****

    vitale232 at gmail.com****

    --
    Oscar Perpi??n Lamigueiro
    Grupo de Sistemas Fotovoltaicos (IES-UPM) Dpto. Ingenier?a El?ctrica
    (EUITI-UPM)
    URL: http://procomun.wordpress.com
    Twitter: @oscarperpinan
    LinkedIn: http://www.linkedin.com/in/oscarperpinan

    _______________________________________________
    R-sig-Geo mailing list
    R-sig-Geo at r-project.org
    https://stat.ethz.ch/mailman/listinfo/r-sig-geo
    _______________________________________________
    R-sig-Geo mailing list
    R-sig-Geo at r-project.org
    https://stat.ethz.ch/mailman/listinfo/r-sig-geo

    --
    Edzer Pebesma
    Institute for Geoinformatics (ifgi), University of M?nster Weseler Stra?e 253, 48151 M?nster, Germany. Phone: +49 251
    83 33081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de


    _______________________________________________
    R-sig-Geo mailing list
    R-sig-Geo at r-project.org
    https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Related Discussions

Discussion Navigation
viewthread | post

2 users in discussion

Zia Uddin Ahmed: 2 posts Edzer Pebesma: 1 post