[R-sig-Geo] Changing pixel values in a stack file based on regular patterns
Robert J. Hijmans
r.hijmans at gmail.com
Sat Dec 18 20:04:39 CET 2010
Dear Paulo,
I think your example works for me:
#function (very basic) -- sure can be improved...
f = function(x)
{
if(any(x%in%1) & any(x%in%3))#if land uses 1 and 3 are present in
that pixel over time...
{
r.0 = grepl("3",x)#find the postion of LU3
r.M1 = grepl("1",x[which(r.0)-1])#pixel yr-1 of LU3 with LU1
r.P1 = grepl("1",x[which(r.0)+1])#pixel yr+1 of LU3 with LU1
if(any(any(r.0)&any(r.M1&r.P1)))#if there is a LU3 somewhere
and LU1 before and after that
{
x[r.0]=1 #replace the LU3 by LU1
}
}
return(x)
}
library(raster)
r <- r2 <- r3 <- r4 <- raster(ncol=5, nrow=5)
set.seed(23)
r[] <- as.numeric(sample(rep(1:3,1000),25))
r2[] <- as.numeric(sample(rep(1:3,1000),25))
r3[] <- as.numeric(sample(rep(1:3,1000),25))
r4[] <- as.numeric(sample(rep(1:3,1000),25))
rr = stack(r,r2,r3,r4)
ee = calc(rr, fun=f)
rv <- getValues(rr)
ev <- getValues(ee)
# which rows have different values?
dif <- rv != ev
rows <- which(apply(dif, 1, sum) > 0)
> rv[rows,]
layer1 layer2 layer3 layer4
[1,] 1 3 1 2
[2,] 1 3 1 3
[3,] 3 1 1 3
[4,] 1 3 1 1
> ev[rows,]
[1,] 1 1 1 2
[2,] 1 1 1 1
[3,] 1 1 1 1
[4,] 1 1 1 1
If this does not work for you, perhaps you need to update the raster package.
Here is an alternative formulation of your function that might be useful
f = function(x) {
# changes 131 to 111, but does not touch 311 or 113
a <- which(x == 3)
a <- subset(a, a > 1 & a < length(x))
b <- which(x[a-1] ==1 & x[a+1] == 1)
x[a[b]] <- 1
return(x)
}
Robert
On Sat, Dec 18, 2010 at 9:39 AM, Paulo Brando <paulobrando at gmail.com> wrote:
> Dear All,
>
> I'm trying to improve a land-use classification based on MODIS product from
> 2000-2010 (one map per year). Each pixel was classified into 0 (no class), 1
> (LU1), 2 (LU2), or 3 (LU3). The overall result of this classification is
> quite good, but there are some cases that do not make much sense. For
> example, some pixels have LU1 in year one, LU3 in year 2, and LU1 again in
> year 3, followed by LU1 during the remaining of the time-series.
>
> I wrote a simple function to deal with this problem, but haven't
> figured out why it is not working with calc.
>
> #function (very basic) -- sure can be improved...
> f = function(x)
> {
> if(any(x%in%1) & any(x%in%3))#if land uses 1 and 3 are present in that
> pixel over time...
> {
> r.0 = grepl("3",x)#find the postion of LU3
> r.M1 = grepl("1",x[which(r.0)-1])#pixel yr-1 of LU3 with LU1
> r.P1 = grepl("1",x[which(r.0)+1])#pixel yr+1 of LU3 with LU1
>
> if(any(any(r.0)&any(r.M1&r.P1)))#if there is a LU3 somewhere and LU1
> before and after that -- problematic when there are two 131 patterns...
> {
> x[r.0]=1#replace the LU3 by LU1
> return(x)}
> else{return(x)}
> }
> else{return(x)}
> }
>
> #testing the function:
> x = c(2,2,2,0,1,1,3,1,1,1)
> f(x)
> #[1] 2 2 2 0 1 1 1 1 1 1
>
> #it seems to work (at least when there is only one series of 1-3-1), but
> when I include it in the calc #function, something unexpected happens.
>
> #Example:
> library(raster)
>
> #creating rasters...
> r <- raster(ncol=5, nrow=5)
> r2 <- raster(ncol=5, nrow=5)
> r3 <- raster(ncol=5, nrow=5)
> r4 <- raster(ncol=5, nrow=5)
>
> set.seed(23)
> r[] <- as.numeric(sample(rep(1:3,1000),25))
> r2[] <- as.numeric(sample(rep(1:3,1000),25))
> r3[] <- as.numeric(sample(rep(1:3,1000),25))
> r4[] <- as.numeric(sample(rep(1:3,1000),25))
>
> rr = stack(r,r2,r3,r4)
>
> #############
> f = function(x)
> {
> if(any(x%in%1) & any(x%in%3))#if land uses 1 and 3 are present in that
> pixel over time...
> {
> r.0 = grepl("3",x)#find the postion of LU3
> r.M1 = grepl("1",x[which(r.0)-1])#pixel yr-1 of LU3 with LU1
> r.P1 = grepl("1",x[which(r.0)+1])#pixel yr+1 of LU3 with LU1
>
> if(any(any(r.0)&any(r.M1&r.P1)))#if there is a LU3 somewhere and LU1
> before and after that
> {
> x[r.0]=1#replace the LU3 by LU1
> return(x)}
> else{return(x)}
> }
> else{return(x)}
> }
>
> x = c(2,2,2,0,1,1,3,1,1,1)#example of land use over time for a single pixel
> (each value represents a yr)
>
> f(x)
>
> ############################
> ee = calc(rr, fun=f)
>
> ss = stack(rr,ee)
> plot(ss)
>
> #it removes all LU2?
>
> Thanks, Paulo
>
> [[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