jointure spatiale sur deux fonctionnalités simples {sf} avec plus de 1 mil. entrées aussi vite que possible

Question

J'espère que ce n'est pas trop trivial mais je ne peux vraiment pas trouver une réponse et je suis trop nouveau sur le sujet pour trouver des alternatives moi-même. Voici donc le problème:

J'ai deux fichiers de formes x et y qui représentent différents niveaux de traitement d'une image satellite Sentinel2. x contient environ 1.300.000 polygones/Segments couvrant complètement l'image s'étendent sans autres informations vitales.
y contient environ 500 polygones représentant la zone sans nuage de l'image (couvrant également la majeure partie de l'image à l'exception de quelques "trous de nuages") ainsi que des informations sur l'image utilisée dans 4 colonnes (Capt .. Si vous essayez d'ajouter l'information d'image à x à des endroits, x est couvert par y. assez simple? Je ne peux pas trouver un moyen de le faire hapen sans prendre des jours. J'ai lu x dans comme une caractéristique simple {sf}, comme le lire avec shapefile/readOGR prend des âges. J'ai essayé différentes choses avec y

lorsque j'essaie de fusionner (x, y) je ne peux prendre qu'une seule sf car la fusion ne supporte pas deux sf. Fusionner x (comme sf) et y (comme shp) me donne l'erreur "impossible d'allouer un vecteur de taille 13.0 Gb"

alors j'ai essayé sf :: st_join (x, y), qui supporte les deux variables à être sf mais n'a pas fini pour 28 heures maintenant

sf :: st_intersect (x, y) a pris environ 9 minutes pour un sous-ensemble de 10.000 segments, ce qui peut ne pas être beaucoup plus rapide pour l'ensemble de la pièce.

est-ce que x peut être résolu en quelques parties plus petites ou y a-t-il une autre solution simple? Puis-je faire quelque chose avec mon espace de travail pour faire fonctionner la fusion ou n'y a-t-il simplement aucun raccourci pour rejoindre cette quantité de polygones?

Merci beaucoup d'avance et j'espère que ma description n'est pas trop floue!

mon petit poste de travail:

win 7 64 bits 8 Go de RAM Intel i7-4790 @ 3,6 GHz

Cheers, Matthias

Answer 1

Je sont souvent confrontés à ce genre de problèmes et comme @ manotheshark2 affirme, je préfère travailler dans une boucle qui englobe ma couche vectorielle. Voici mon conseil:

Chargez vos données

library(raster) 
library(rgdal) 
x <- readOGR('C:/', 'sentinelCovers') 
y <- readOGR('C:/', 'cloudHoles') 

Assigner un ID y pour identifier les x polygones intersecte les polygones y et créer la colonne dans le tableau x

x$xyID <- NA # Answer col 
y$yID <- 1:nrow(y@data) # ID col 

Exécuter une boucle subseting x

for (posX in 1:nrow(x@data)){ 
    pol.x <- x[posX, ] 
    intX <- raster::intersect(pol.x, y) 
    # x$xyID[posX] <- intX@data$yID ## Run this if there's unique y polygons 
    # x$xyID[posX] <- paste0(intX@data$yID, collapse = ',') ## Run this if there's multiple y polygons 
} 

Vous pouvez vérifier si est préférable d'exécuter la boucle sur la couche de xoy

x$xyID <- NA # Answer col 
x$xID <- 1:nrow(x@data) # ID Col 

for (posY in 1:nrow(y@data)){ 
    pol.y <- y[posY, ] 
    intY <- tryCatch(raster::intersect(pol.y, x), finally = 'NULL') 
    if (is.null(intY)) next 
    x$xyID[x@data$xID %in% intY@data$xID] <- pol.y$yID 
}