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.