Assigner un polygone à un point de données dans un cadre de données R

Question

J'ai deux cadres de données :

• points contient une série de points avec x, y coordonnées.
• poly contient les coordonnées de deux polygones (j'en ai plus de 100 en réalité, mais je reste simple ici).

Je veux pouvoir ajouter à la trame de données points une colonne supplémentaire appelée Area qui contient le nom du polygone dans lequel se trouve le point.

poly <- data.frame(
pol= c("P1", "P1","P1","P1","P1","P2","P2","P2","P2", "P2"),
x=c(4360, 7273, 7759, 4440, 4360, 8720,11959, 11440,8200, 8720),
y=c(1009, 9900,28559,28430,1009,9870,9740,28500,28040,9870))

points <- data.frame(
       object = c("P1", "P1","P1","P2","P2","P2"),
       timestamp= c(1485670023468,1485670023970, 1485670024565, 1485670025756,1485670045062, 1485670047366),
       x=c(6000, 6000, 6050, 10000, 10300, 8000),
       y=c(10000, 20000,2000,5000,20000,2000))

plot(poly$x, poly$y, type = 'l')
text(points$x, points$y, labels=points$object )

Donc, essentiellement dans cet exemple, les 2 premières lignes devraient avoir Area= "P1" tandis que le dernier point doit être vide car il n'est contenu dans aucun polygone.

J'ai essayé d'utiliser la fonction in.out mais je n'ai pas été capable de construire mon cadre de données comme je l'ai décrit.

Toute aide est très appréciée !

Answer 1

Bien que cela utilise un for boucle, il est pratiquement assez rapide.

library(mgcv)

x <- split(poly$x, poly$pol)
y <- split(poly$y, poly$pol)

todo <- 1:nrow(points)
Area <- rep.int("", nrow(points))
pol <- names(x)

# loop through polygons
for (i in 1:length(x)) {
  # the vertices of i-th polygon
  bnd <- cbind(x[[i]], y[[i]])
  # points to allocate
  xy <- with(points, cbind(x[todo], y[todo]))
  inbnd <- in.out(bnd, xy)
  # allocation
  Area[todo[inbnd]] <- pol[i]
  # update 'todo'
  todo <- todo[!inbnd]
  }

points$Area <- Area

Deux raisons pour son efficacité :

• for La boucle passe par les polygones, pas par les points. Ainsi, si vous avez 100 polygones et 100 000 points à allouer, la boucle n'a que 100 itérations et non 100 000. À l'intérieur de chaque itération, la puissance de vectorisation de la fonction C in.out est exploité ;
• Il fonctionne de manière progressive. Une fois qu'un point a été attribué, il sera exclu de l'attribution ultérieurement. todo contrôle les points à allouer dans la boucle. Au fur et à mesure, l'ensemble de travail se réduit.

Answer 2

Vous pouvez utiliser la fonction point.in.polygon du paquet sp :

points$Area = apply(points, 1, function(p)ifelse(point.in.polygon(p[3], p[4], poly$x[which(poly$pol==p[1])], poly$y[which(poly$pol==p[1])]), p[1], NA))

vous donne

  object   timestamp     x     y Area
1     P1 1.48567e+12  6000 10000   P1
2     P1 1.48567e+12  6000 20000   P1
3     P1 1.48567e+12  6050  2000 <NA>
4     P2 1.48567e+12 10000  5000 <NA>
5     P2 1.48567e+12 10300 20000   P2
6     P2 1.48567e+12  8000  2000 <NA>