points-in-polygons calculation

Hi all,

during a spatial analysis I need to compute a point-in-polygon
calculation, and I chose findPolys() in PBSmapping package.

For the moment the script is built as follows:
- import the polygon shapefile with rgdal
- import the point shapefile with PBSmapping
- nested for cycles that
 - isolate each polygon (or region) through a superkey on the
attributes in data slot. The selected element is a homerange. It's then
converted in a PolySet, for next step
 - finds which points fall inside the homerange with findPolys()
 - picks up a defined value among the attributes of points and computes
mean of all values of the selected points, then writes it in a new
column of the attribute table

Well, if you showed the code used for doing this, it might be easier to 
advise. Are you using the maptools function SpatialPolygons2PolySet()? You 
can then use the PID and SID values as you would otherwise in PBSmapping.

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.

The IDs of the Polygons in the SpatialPolygons* object SP are in:

sapply(slot(SP, "polygons"), function(x) slot(x, "ID"))

Hope this helps,

Roger

Can the for cycle be based on a unique key like PolygonsID? Among the
methods for sp classes I didn't find it. I found functions like
getPolygonsIDSlot but they have not to be used directly.

Any suggestions?

thank you in advance

Anne Ghisla
Universit? degli Studi dell'Insubria, Varese

Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no

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>