Error in if (dxsd > tolerance) { : missing value where TRUE/FALSE needed

Dear R-sig-Geo'ers

I am currently trying to extract some NDVI values from a group of NDVI images for specific shapefile polygons. I have 17 individual polygons and 23 NDVI images. I have written the script below to select individual polygons and to then extract the average NDVI within each polygon and store this in a list. Once all 23 images have been processed the script moves to the next polygon and so on and so on, until all polygons have been average for all 23 images. The code essentially extracts the NDVI data as a time series for each polygon.

Unfortunately, as everyone knows, Landsat 7 ETM has a broken scan line corrector meaning that all post 2003 data have large areas of missing data. In my pre-processing I convert these areas to nan using gdal_merge.py (I have to stich and subset two scenes). The code below works well until a large portion of a polygon lies on the nan area of the image (some polygons that only have a small area in the nan process fine). When this happens I get the following error message

Error in if (dxsd > tolerance) { : missing value where TRUE/FALSE needed

The error is related to the following command

x.clip <- temp[!is.na(overlay(temp, brede_add)),]

I have no idea how to solve this error, as a backup plan I am using r.fillnulls in grass to get rid of the stripes but this takes a very long time and is not ideal.

Any additional help or pointers would be greatly appreciated

Many thanks and kind regards

****************************************************************************************************************************************************************************************************************************************************

library(maptools)
library(rgdal)
library(spatstat)

input <- "/media/win-drive/wes2006/Projects/Evapotranspiration/Data/Landsat_NDVI"
maps <- data.frame((list.files(input, pattern=".img",full.names=FALSE)),(list.files(input, pattern=".img", full.names=TRUE)))
names(maps) <- c("File", "Path")
fn=data.frame(maps$File)

rastx<-284055.000
rasty<-6302475.000

Lst <- list()

brede <- readShapePoly("BredeCatchmentSites.shp", IDvar="COUNT", proj4string=CRS("+proj=utm +zone=34 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))

x=0
while (x <= 17)
{
? ? ? ?x<-x+1
? ? ? ? ? ? ? ?nfiles <- length(maps$Path)
? ? ? ? ? ? ? ?for (n in 1:nfiles)
? ? ? ? ? ? ? ?{
? ? ? ? ? ? ? ? ? brede_add <- brede[brede$COUNT[x], ]
? ? ? ? ? ? ? ? ? box <- as.list(as.numeric(bbox(brede_add)))
? ? ? ? ? ? ? ? ? xstart <- (as.numeric(box[1])-rastx)/30
? ? ? ? ? ? ? ? ? ystart <- abs((as.numeric(box[4])-rasty)/30)

? ? ? ? ? ? ? ? ? row <- round((as.numeric(box[3])-as.numeric(box[1]))/30)
? ? ? ? ? ? ? ? ? col <- round((as.numeric(box[4])-as.numeric(box[2]))/30)

? ? ? ? ? ? ? ? ? temp <- readGDAL(as.character(maps$Path[n]),silent=TRUE, offset=c(ystart,xstart), region.dim=c(col,row))
? ? ? ? ? ? ? ? ? fullgrid(temp)=FALSE

? ? ? ? ? ? ? ? ? x.clip <- temp[!is.na(overlay(temp, brede_add)),] ? <<<<<<<<<<<-------------HERE IS THE PROBLEM-------------

? ? ? ? ? ? ? ? ? a <- mean(as.data.frame(x.clip$band1))

? ? ? ? ? ? ? ? ? Lst[n] <- list(as.numeric(a))
? ? ? ? ? ? ? ? ? print(n)
? ? ? ? ? ? ? ? ? x.clip<-0
? ? ? ? ? ? ? ? }

? ? ? ?data <- as.numeric(Lst)
? ? ? ?data2 <- cbind(data,fn)
? ? ? ?write.table(data, paste(brede$NBALID[x]))
? ? ? ?print(paste(brede$NBALID[x]))
}



Wesley Roberts MSc.
Researcher: Earth Observation (Ecosystems)
Natural Resources and the Environment
CSIR
Tel: +27 (21) 888-2490
Fax: +27 (21) 888-2693

"To know the road ahead, ask those coming back."
- Chinese proverb


--
This message is subject to the CSIR's copyright terms and conditions, e-mail legal notice, and implemented Open Document Format (ODF) standard.
The full disclaimer details can be found at http://www.csir.co.za/disclaimer.html.


and is believed to be clean. ?MailScanner thanks Transtec Computers for their support.