[R-sig-Geo] polygon/raster overlay error

[Tomislav Hengl]

Evann Smith wrote:
> Hello,
> 
> I am attempting to calculate the mean value in a raster (3432 x 8640) for
> each of 689 polygons (code at end of email). To that end, I have a few
> questions:
> 
> First, the for loop runs fine for a while (2 hours, maybe) before
> terminating itself and returning this error:
> 
> Error in unlist(lapply(rings, function(x, pts) pointsInPolygon(pts, x,  :
>   negative length vectors are not allowed
> 
> Based on the output matrix, it seems to have gotten hung up on polygon 24.
> Any ideas why this might be and how I can fix it? In reading the shapefile I
> both delete null objects and force rings.
> 
> Second, this process (if it works) will take forever. I've tried an approach
> where I covert the polygons to a raster, but this seems to result in data
> loss and produce a lot of null values. I've also read on this list a
> discussion about just using those points in the bounding box. However, it
> was entirely unclear to me how to implement that recommendation.

Maybe you should try using the SAGA module "Grid Statistics for 
Polygons". SAGA (implemented in C++) is pretty efficient in processing 
large data. I would be interested to learn how much time did it take on 
the end with your specific dataset.

 > library(RSAGA)
 > rsaga.get.usage("shapes_grid", 2)
SAGA CMD 2.0.4
library path:   C:/Progra~1/saga_vc/modules
library name:   shapes_grid
module name :   Grid Statistics for Polygons
Usage: 2 -GRID <str> -POLY <str> -RESULT <str> [-QUANTILES] 
[-QUANTILE_STEP <num>]
   -GRID:<str>           Grid
         Grid (input)
   -POLY:<str>           Polygons
         Shapes (input)
   -RESULT:<str>         Result
         Shapes (output)
   -QUANTILES            Quantiles
         Boolean
   -QUANTILE_STEP:<num>  Quantile Step
         Choice
         Available Choices:
         [0] median
         [1] quartiles
         [2] deciles


T. Hengl
http://home.medewerker.uva.nl/t.hengl/

> 
> Any help would be greatly appreciated.
> 
> Cheers,
> Evann
> 
> 
> R Code:
> # Projection
> prj <- CRS("+proj=longlat +ellps=WGS84")
> 
> # Shapefile
> polys <- readShapePoly("GeoEPR/geoepr-20100212-1248.shp", proj4string=prj,
> verbose=TRUE, delete_null_obj=TRUE, force_ring=TRUE)
> summary(polys)
> 
> # Raster
> popimg <- readGDAL("PopDens/glds05ag.bil")
> proj4string(popimg) <- prj
> summary(popimg)
> 
> # Values
> y <- slot(popimg,"data")
> head(y)
> 
> # Total polygons
> n <- length(slot(polys,"polygons"))
> n
> 
> # Empty matrix
> output1 <- as.data.frame(matrix(NA, nrow=n, ncol=4))
> names(output1) <- c("cowgroup", "start", "end", "popdens")
> head(output1)
> 
> # Beastly for loop
> for(i in 1:n) {
>    output1[i,1] <- slot(polys, "data")[i,1]
>    output1[i,2] <- slot(polys, "data")[i,2]
>    output1[i,3] <- slot(polys, "data")[i,3]
>    overlap <- overlay(popimg,polys[i,])
>    output1[i,4] <- mean(y[overlap==1 & !is.na(overlap),1], na.rm=TRUE)
>    }
> 
> ============================================
> Evann Smith
> Doctoral Student
> Department of Government, Harvard University
> ============================================
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>