[R] IFELSE across large array?

[Sander Oom]

Dear all,

Thanks to the help from the R mailing list we now have a script performing 
the required tasks, i.e. applying a mask and filter to a 3D array. See 
script below. Any further suggestions to simplify the code are very 
welcome. We are currently debugging the larger script of which this section 
is a small part.

The script could serve as a framework for many typical GIS neighborhood 
analyses, in this case the 'block majority'. The 'var.range' variable is 
derived from variography earlier in the main script.

Yours,

Sander and Alessandro

##########INPUT SECTION1##############################
numRow<-7
numCol<-8
numReal<-2
var.range<-  2.1
cellsize<- 0.25

################################################
##CALCULATE WINDOW SIZES BASED ON INPUT ABOVE#####
#window size ...minimum is 1

numCells<-round((var.range*cellsize)/2,0)
if(numCells < 1) {numCells<-1}

###################################################
library(abind)

#Sim: 0=Absent; 1=Present; 10=NA
vectSim <- c(1,1,0,0,1,1,10,1,0,10,1,1,1,0,1,1,1,1,1,0,0,1,0,0,10,0,0,
1,1,1,0,1,1,0,0,0,1,1,1,0,1,1,0,1,0,1,10,0,1,0,0,1,0,0,10,0,0,0,1,1,0,
0,1,0,1,0,0,1,0,0,1,1,1,1,1,10,1,10,1,0,1,0,1,0,0,1,1,0,10,10,1,1,0,1,
0,1,1,1,0,0,0,1,0,0,0,0,0,0,1,1,0,0)
#Mask: 10=NA; 1=DATA
vectMask <- c(10,1,10,1,1,10,1,10,1,10,10,1,10,1,10,1,1,10,1,10,1,10,10,1,10,1,
10,1,1,10,1,10,1,10,10,1,10,1,10,1,1,10,1,10,1,10,10,1,10,1,10,1,1,10,1,10)
length(vectSim)
length(vectMask)
Sim <- array(vectSim, c(numRow,numCol,numReal))
Sim[Sim==10]<-NA
Sim
Mask <- array(data=vectMask, c(numRow,numCol,numReal))
Mask

##expand array
junkcol<-array(rep(NA,numReal*numRow),c(numRow,numCells,numReal))
#add columns borders

#cols
Sim<-abind(junkcol,Sim,junkcol, along=2)
dim(Sim)[2]
junkrow<-array(rep(NA,dim(Sim)[2]),c(numCells,dim(Sim)[2],numReal))
Sim<-abind(junkrow,Sim,junkrow,along=1)

##DEFINE PORTION OF ARRAY ON WHICH MOVING WINDOW RUNS
##IT AVOIDS THE na BORDER

#maximum row and col indexes to consider are equal
# to num num rows and num cols in the original matrix

minr<-1+numCells
maxr<-dim(Sim)[1]-numCells
minc<-1+numCells
maxc<-dim(Sim)[2]-numCells

clean<-function(a,nr,r,c) {
   rmin<-row(a)[r-numCells]
   rmax<-row(a)[r+numCells]
   cmin<-col(a)[nr*(c-numCells)]
   cmax<-col(a)[nr*(c+numCells)]

   junk<-a[rmin:rmax,cmin:cmax]
   junk[numCells,numCells]<-NA
   junk<-as.vector(junk)
   sample(na.omit(junk),1)
}

for (n in 1:numReal)  {
   realiz<-Sim[1:dim(Sim)[1],1:dim(Sim)[2],n]
   rz<-realiz
   for (r in minr:maxr)  {
     for (c in minc:maxc)  {
       if(is.na(realiz[r,c]) ) {
         realiz[r,c]<-clean(a=realiz,nr=numRow,r=r,c=c)
         Sim[1:dim(Sim)[1],1:dim(Sim)[2],n]<-realiz
       }#end if
     }
   }
}##end overall loop

# Return to original array dimensions, removing extra NA
minr <- (numCells+1)
maxr <- (dim(Sim)[1]-numCells)
minc <- (numCells+1)
maxc <- (dim(Sim)[2]-numCells)

expSim <- Sim
expSim
Sim <- expSim[minr:maxr, minc:maxc, 1:dim(Sim)[3]]
Sim

# Clip simulation results using the mask
Mask

Sim[Sim==10] <- NA
Mask[Mask==10] <- NA
A <- Sim + Mask -1
A

#################################




At 01:55 2004/11/24, Ray Brownrigg wrote:
> > From: "Liaw, Andy" <andy_liaw at merck.com>
> > Date: Tue, 23 Nov 2004 12:28:48 -0500
> >
> > I'll give it half a crack:
> >
> > Steps a through c can be done via nested ifelse(), treating A and M as
> > vectors (as they really are).  Step d is the hard one.  I'd write a simple
> > Fortran code and use .Fortran() for that.
> >
> > I don't see how any of the *apply() functions can help here, as your
> > operations are element-wise, not dimension-wise.
> >
> > Andy
> >
>The original message mentions that the value 10 is actually "NODATA",
>and if one recodes 10 as NA, steps a) to c) become trivial, namely:
>
>A[A == 10] <- NA
>M[M == 10] <- NA
>return(A + M - 1)
>
>Then if step d) is performed first (i.e. appropriate values in A are
>replaced by the 'most common neighbour' [perhaps using
>round(mean(.., na.rm=T))] this still works, but would have to be repeated
>for each replication (the third dimension).
>
>Ray Brownrigg
>
> > > From: Sander Oom
> > >
> > > Dear all,
> > >
> > > As our previous email did not get any response, we try again with a
> > > reformulated question!
> > >
> > > We are trying to do something which needs an efficient loop
> > > over a huge
> > > array, possibly functions such as apply and related (tapply,
> > > lapply...?), but can't really understand syntax and examples in
> > > practice...i.e. cant' make it work.
> > >
> > > to be more specific:
> > > we are trying to apply a mask to a 3D array.
> > > By this I mean that when "overlaying" [i.e. comparing element
> > > by element]
> > > the mask on to the array the mask should change array
> > > elements according to
> > > the values of both array and mask elements
> > >
> > > the mask has two values: 1 and 10.
> > >
> > > the array elements have 3 values: 0, 1,  or 10
> > >
> > > sticking for the moment to the single 2d array case
> > >
> > > for example:
> > > [A= array]  10    0 10 1  10  0
> > >                   1   10   1 0   0 10
> > >
> > > [ M=mask]     1  10  10 1   1  1
> > >                  10    1   1  1 10 10
> > >
> > > I would like the array elements to:
> > >
> > > a) IF A(ij) !=10 and  Mij = 1
> > >               leave A(ij) unchanged
> > >
> > > b)  IF   A(ij) != 10 and M(ij) =10
> > >                 change A(ij) to M(ij) i.e mask value (10)
> > >
> > > c)IF A(ij) = 10 and M(ij) = 10
> > >                leave (Aij) unchanged
> > >
> > > d) IF A(ij) = 10 and M(ij) !=10
> > >             replace A(ij) with the majority value in the
> > > 8-neighborhood
> > >
> > > (or whatever if it is an edge element) BUT ignoring 10s in this
> > > neighborhood (i.e. with either 1 or 0, whichever is in majority)
> > >
> > > because the array is 3d I would like to repeat the thing with
> > > all the k
> > > elements (2d sub-arrays) of the array in question, using the
> > > same mask for
> > > al k elements
> > >
> > > Would you be able to suggest a strategy to do this?
> > >
> > > thanks very much
> > >
> > > Alessandro and Sander.
>
>______________________________________________
>R-help at stat.math.ethz.ch mailing list
>https://stat.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html