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
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,[EMAIL PROTECTED]) == 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)
[EMAIL PROTECTED] == y & polygons$MONTH
== m & polygons$ID == i,"EN_HR"] <- energyHR
}
, silent = TRUE)
}
}
}
writeOGR(polygons, 'dsn', 'namefile', driver='ESRI Shapefile')
signature.asc
Description: Questa รจ una parte del messaggio firmata digitalmente
_______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
