[R-sig-Geo] Error in if (dxsd > tolerance) { : missing value where TRUE/FALSE needed

[Roger Bivand]

On Wed, 18 Mar 2009, David Wahlund wrote:

> Hi
> I got the same error trying to do something similar. For me the
> problem was with the polygon data. One of the shapes was corrupted. My
> tip is to try and identify what polygon renders the error and confirm
> that it is OK.

Yes, the underlying problem is that the "[" method for a SpatialGrid is 
trying to re-assemble an object from the remaining data, and the heuristic 
being used to "guess" the step of the grid when many rows and columns are 
missing is failing. If you manage to find out which it is, you might 
consider sending me a copy of the objects (save()), so that I can have a 
look at the heuristic. In general, you can drill down to the failing 
function by using traceback() afterwards - here it is points2grid() that 
has failed.

Roger

>
> -David Wahlund
>
> On Wed, Mar 18, 2009 at 16:01, Wesley Roberts <wroberts at csir.co.za> wrote:
>> 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.
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no