[R-sig-Geo] Help: Individual plots all raster in a raster stack in loop
Zia Uddin Ahmed
zua3 at cornell.edu
Wed Jul 3 19:18:09 CEST 2013
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 <vitale232 at gmail.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 <zua3 at cornell.edu> 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 <zua3 at cornell.edu> 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 at 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> [mailto:
>> romunov at gmail.com<mailto:romunov at gmail.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 at 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<mailto:zua3 at cornell.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.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 at 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(nc=1000, nr=1000)
>> 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=23 +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(nc=1000, nr=1000)
>> >> 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=23 +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 <mailto:alexandre.santos at cas.ifmt.edu.br>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> 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 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<mailto:R-sig-Geo at r-project.org><mailto:
>> R-sig-Geo at r-project.org<mailto:R-sig-Geo at 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 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<tel:%28412%29%20915-3632 <%28412%29%20915-3632>>
>> vitale232 at gmail.com<mailto:vitale232 at gmail.com>
>>
>>
>>
>> --
>> Andrew P. Vitale
>> Masters Student
>> Department of Geography
>> University of Nevada, Reno
>> (412) 915-3632
>> vitale232 at gmail.com<mailto:vitale232 at gmail.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
More information about the R-sig-Geo
mailing list