Many thanks Edzer, gridded is now working a charm on my data.
Regards, Wesley >>> Edzer Pebesma <edzer.pebe...@uni-muenster.de> 03/18/09 9:58 PM >>> A quick solution could be to convert it into a SpatialPointsDataFrame object (gridded(x) = FALSE), before doing the selection that broke your script. If it still breaks, you might check whether the selection is empty. -- Edzer Roger Bivand schrieb: > 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 <wrobe...@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@stat.math.ethz.ch >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>> >> >> _______________________________________________ >> R-sig-Geo mailing list >> R-sig-Geo@stat.math.ethz.ch >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > > ------------------------------------------------------------------------ > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@stat.math.ethz.ch > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- 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@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo