[R-sig-Geo] polygon/raster overlay error
Tomislav Hengl
hengl at spatial-analyst.net
Mon Apr 19 09:28:15 CEST 2010
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
>
More information about the R-sig-Geo
mailing list