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" <[EMAIL PROTECTED]> > 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.
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
