[R-sig-Geo] points-in-polygons calculation

[Anne Ghisla]

Il giorno mar, 12/02/2008 alle 18.57 +0100, Roger Bivand ha scritto:

> Well, if you showed the code used for doing this, it might be easier to 
> advise. 
you're right. the code is in attachment. It's still a beta version, and
needs much improvement.
the attributes that identify an home range are year, period and ID
together. that's why the code contains three nested for cycles.

> But it does strike me that it might be simpler just to read the points 
> with rgdal, and use the overlay() method in sp to find which points belong 
> to which polygons, all in one go. Using the output membership vector, you 
> could simply run by() or aggregate on as(pts, "data.frame") to get your 
> summaries per home range. The only reason for caution might be if the home 
> range polygons overlap, in which case judicious use of lapply() on the 
> "polygons" slot of the SpatialPolygons* might be needed, because some 
> points might fall into multiple polygons.
in fact most polygons overlap, therefore I found easier to work on one
polygon at time with PBSmapping functions.

> The IDs of the Polygons in the SpatialPolygons* object SP are in:
> 
> sapply(slot(SP, "polygons"), function(x) slot(x, "ID"))

thank you for hints and explanations!

regards,
Anne

-------------- next part --------------

rm(list=ls())

require(rgdal) 
require(PBSmapping) 

polygons <- readOGR('path2shapefiledsn', 'shapefileName') # import con rgdal
points <- importShapefile('path2shapefileNameWithoutExt')  # import con pbsmapping 

[...]

# correspondence table. the energy values to be picked up for each homerange depend on the its time frame
rules <- read.table(\"$file\", header=TRUE, sep=\"\t\", dec=\".\") #per esempio
rules <- data.frame(
	year = c('2000','2002','2002','2004','2004','2006','2006'),
    seas = c('spr','spr','aut','spr','aut','spr','aut'),
	energy = c('EN_F1999','EN_F2001','EN_F2002','EN_F2003','EN_F2004','EN_F2005','EN_F2006')
)

for (y in as.character(unique(polygons$YEAR))) {
	for (p in as.character(unique(polygons$PERIOD))) {
		for (i in as.character(unique(polygons$ID))) {
			try(
				if (nrow(polygons[polygons$YEAR == y & polygons$PERIOD == p & polygons$ID == i,]@data) == 1) {
					poly.tmp <- polygons[polygons$YEAR == y & polygons$PERIOD == p & polygons$ID == i,]
					poly.PBS.tmp <- SpatialPolygons2PolySet(poly.tmp)
					intersection <- findPolys(points, poly.PBS.tmp)
					fld2 <- as.character(rules[rules$year == y & rules$seas == p,]$energy)
					energyHR <- mean(points[points$EID %in% intersection$EID,fld2],na.rm = TRUE)
					polygons at data[polygons$YEAR == y & polygons$MONTH == m & polygons$ID == i,"EN_HR"] <- energyHR
				}
				, silent = TRUE)
		}
	}
}
writeOGR(polygons, 'dsn', 'namefile', driver='ESRI Shapefile')

-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: Questa ? una parte del messaggio	firmata digitalmente
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20080213/203e0697/attachment.bin>