[R-sig-Geo] Changing pixel values in a stack file based on regular patterns

[Robert J. Hijmans]

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
>