Comment fonctionne le polygone spatial% sur% polygone lors de l'agrégation des valeurs dans r?

12

Je travaille sur un projet d'épidémiologie environnementale où j'ai des expositions ponctuelles (~ 2 000 porcs industriels - OHI). Ces OHI pulvérisent sur les champs voisins, mais les gouttelettes d'eau et les odeurs des matières fécales peuvent parcourir des kilomètres. Donc, ces expositions ponctuelles obtiennent des tampons de 3 mi, et je veux connaître le nombre d'expositions de l'OHI (de différents types - somme de la quantité de fumier, nombre de porcs, peu importe; le plus simple, juste le nombre de tampons d'exposition qui se chevauchent) par blocs de recensement NC (~ 200 000). Les blocs de recensement d'exclusion (en bleu) sont (1) n'importe quoi dans les 5 villes les plus peuplées et (2) les comtés qui ne bordent pas un comté avec une OHI (note: cela a été fait avec la fonction gRelate et les codes DE-9IM - très lisse!). Voir l'image ci-dessous pour un visuel

entrez la description de l'image ici

La dernière étape consiste à agréger la représentation de l'exposition tamponnée à chaque bloc de recensement. Voici où je suis perplexe.

Jusqu'à présent, j'ai passé de bons moments avec les fonctions% over% dans le package sp, mais comprenez d'après la vignette over que poly-poly et poly-line over sont implémentés dans les rgeos. La vignette ne couvre que le poly-ligne et le poly auto-référencé, et non l'agrégation, donc je suis un peu confus sur mes options pour le poly-poly avec l'agrégation de fonctions, comme la somme ou la moyenne.

Pour un cas de test, considérez l'extrait ci-dessous, quelque peu verbeux, qui fonctionne avec le fichier des frontières des pays du monde. Cela devrait pouvoir être copié et exécuté tel quel, puisque j'utilise une graine aléatoire pour les points et que je télécharge et décompresse le fichier mondial en code.

Tout d'abord, nous créons 100 points, puis utilisons la fonction over avec l'argument fn pour additionner l'élément dans le bloc de données. Il y a beaucoup de points ici, mais jetez un œil à l'Australie: 3 points, le numéro 3 comme étiquette. Jusqu'ici tout va bien.

entrez la description de l'image ici

Maintenant, nous transformons les géométries pour pouvoir créer des tampons, les retransformer et mapper ces tampons. (Inclus sur la carte précédente, car je suis limité à deux liens.) Nous voulons savoir combien de tampons chevauchent chaque pays - dans le cas de l'Australie, à l'œil, c'est 4. Je ne peux pas pour la vie de moi comprendre ce qui se passe mais pour obtenir cela avec la fonction over. Voir mon désordre d'une tentative dans les dernières lignes de code.

EDIT: Notez qu'un commentateur sur r-sis-geo a mentionné la fonction d'agrégat - également référencée sur la question d'échange de pile 63577 - donc un contournement / flux pourrait passer par cette fonction, mais je ne comprends pas pourquoi j'aurais besoin d'aller s'agréger pour le polypole quand fini semble avoir cette fonctionnalité pour d'autres objets spatiaux.

require(maptools)
require(sp)
require(rgdal)
require(rgeos)

download.file("http://thematicmapping.org/downloads/TM_WORLD_BORDERS_SIMPL-0.3.zip", destfile="world.zip")
unzip("world.zip")
world.map = readOGR(dsn=".", "TM_WORLD_BORDERS_SIMPL-0.3", stringsAsFactors = F)
orig.world.map = world.map #hold the object, since I'm going to mess with it.

#Let's create 500 random lat/long points with a single value in the data frame: the number 1
set.seed(1)
n=100
lat.v = runif(n, -90, 90)
lon.v = runif(n, -180, 180)
coords.df = data.frame(lon.v, lat.v)
val.v = data.frame(rep(1,n))
names(val.v) = c("val")
names(coords.df) = c("lon", "lat")
points.spdf = SpatialPointsDataFrame(coords=coords.df, proj4string=CRS("+proj=longlat +datum=WGS84"), data=val.v)
points.spdf = spTransform(points.spdf, CRS(proj4string(world.map)))
plot(world.map, main="World map and points") #replot the map
plot(points.spdf, col="red", pch=20, cex=1, add=T) #...and add points.

#Let's use over with the point data
join.df = over(geometry(world.map), points.spdf,  fn=sum)
plot(world.map, main="World with sum of points, 750mi buffers") #Note - happens to be the count of points, but only b/c val=1.
plot(points.spdf, col="red", pch=20, cex=1, add=T) #...and add points.
world.map@data = data.frame(c(world.map@data, join.df))
#world.map@data = data.frame(c(world.map@data, over(world.map, points.spdf, fun="sum")))
invisible(text(getSpPPolygonsLabptSlots(world.map), labels=as.character(world.map$val), cex=1))
#Note I don't love making labels like above, and am open to better ways... plus I think it's deprecated/ing

#Now buffer...
pointbuff.spdf = gBuffer(spTransform(points.spdf, CRS("+init=EPSG:3358")), width=c(750*1609.344), byid=T)
pointbuff.spdf = spTransform(pointbuff.spdf, world.map@proj4string)
plot(pointbuff.spdf, col=NA, border="pink", add=T)



#Now over with the buffer (poly %over% poly).  How do I do this?
world.map = orig.world.map
join.df = data.frame(unname(over(geometry(world.map), pointbuff.spdf, fn=sum, returnList = F)) ) #Seems I need to unname this...?
names(join.df) = c("val")
world.map@data = data.frame(c(world.map@data, join.df)) #If I don't mess with the join.df, world.map's df is a mess..
plot(world.map, main="World map, points, buffers...and a mess of wrong counts") #replot the map
plot(points.spdf, col="red", pch=20, cex=1, add=T) #...and add points.
plot(pointbuff.spdf, col=NA, border="pink", add=T)
invisible(text(getSpPPolygonsLabptSlots(world.map), labels=as.character(world.map$val), cex=1)) 
#^ But if I do strip it of labels, it seems to be misassigning the results?
# Australia should now show 4 instead of 3.  I'm obviously super confused, probably about the structure of over poly-poly returns.  Help?
Mike Dolan Fliss
la source
Appréciez la redirection - dois-je supprimer d'ici et republier là-bas? Quel est le meilleur coup? Merci.
Mike Dolan Fliss

Réponses:

5

Merci pour la question claire et l'exemple reproductible.

Votre compréhension est correcte, et cela se résume à un bogue dans rgeos :: over, qui a été corrigé il y a un mois mais qui n'est pas encore entré dans une version CRAN. Ce qui suit est une solution de contournement si vous êtes uniquement intéressé par le nombre d'intersections:

world.map$val = sapply(over(geometry(world.map), pointbuff.spdf, returnList = TRUE), NROW)

J'utilise NROWici au lieu de lengthpour que cela fonctionne avec les mauvais rgeos (0.3-8, de CRAN) ainsi que les corrigés (0.3-10, de r-forge). La suggestion précédente d'utiliser

a = aggregate(pointbuff.spdf, world.map, sum)

compte également le nombre d'intersections, mais uniquement avec la version de rgeos fixe installée. Son avantage, outre un nom plus intuitif, est qu'il renvoie directement un Spatialobjet, avec la géométrie de world.map.

Pour faire fonctionner les rgeos 0.3-8, ajoutez

setMethod("over",
    signature(x = "SpatialPolygons", y = "SpatialPolygonsDataFrame"),
        rgeos:::overGeomGeomDF)

à votre script, avant de l'utiliser over.

Edzer Pebesma
la source
Très utile, merci. Je tiens particulièrement à célébrer votre offre d'une solution qui fonctionne avant et après la correction. Pourriez-vous élaborer sur: (1) Quel est le bogue que je rencontre ici-rgeos :: over renvoie une géographie de polygone spatial, pas une trame de données poly spatiale? Certaines fonctions ne renvoient-elles pas simplement des trames de données ...? (2) Comment cela est-il généralement censé fonctionner avec des agrégats et plus? Je suis un peu confus quant à leurs différences et à leurs cas d'utilisation. J'apprécie vraiment votre pesée, merci. Et sidenote: des suggestions pour comprendre le cycle de sortie de CRAN?
Mike Dolan Fliss
Aussi, en ce qui concerne la question d'origine: j'ai besoin de compter le nombre d'expositions, mais j'ai aussi vraiment besoin de les additionner - des choses comme le nombre de porcs dans chaque exposition. Compter les chevauchements est un début ... mais il me semble que la solution dont j'ai besoin est de récupérer les derniers rgeos, oui? Pas moyen de faire cette agrégation fonctionnelle (pas seulement de compter) sans elle?
Mike Dolan Fliss
(1) rgeos :: over pour la signature SpatialPolygons,SpatialPolygonsDataFramedevrait retourner un data.frame, mais retourne un vecteur d'index identique à ce yqu'il aurait été SpatialPolygons. sp::aggregatefait ce que vous faites de plus d'une manière plus conviviale, en renvoyant l' Spatialobjet au lieu de data.frame. Les packages CRAN sont gérés par des bénévoles.
Edzer Pebesma
D'accord, merci Edzer. Il semble que l'agrégat repose sur des rgeos, donc pour obtenir cette fonctionnalité avant le cycle de sortie de CRAN (quoi qu'il en soit), je vais avoir besoin de savoir comment télécharger les rgeos les plus récents et travailler sur cela. Je vous remercie. Et merci pour tout votre travail sur le package !!
Mike Dolan Fliss
Aussi, Edzer, merci beaucoup pour la note sur R-sis-geo. Je ne savais pas où était le meilleur endroit pour publier, donc je suis heureux que le fil pointe maintenant ici.
Mike Dolan Fliss
1

J'ai concocté un sur-remplaçant rapide (et mal codé) en attendant qui crée la trame de données dont j'ai besoin, car ma question n'est pas tout à fait répondue par la solution de comptage uniquement ci-dessus ou "travailler sur les nouveaux rgeos", que je ne suis pas assez habile pour comprendre comment faire.

Cette fonction est clairement (1) incomplète (remarquez comment j'ignore l'argument fn) et (2) inefficace, car j'y arrive sans les puissantes manipulations de tableau de R / sapply ... (clairement je viens d'autres langues sans ce pouvoir) mais honnêtement, je ne comprends toujours pas ce que la structure de la fonction over renvoie (liste de listes ...? Et listes blanches si NA?). Pour ce qu'elle vaut (modifications bienvenues), cette fonction fait le travail que j'ai besoin de faire, avec succès, et imite l'action des autres sur les fonctions.

Modifications bienvenues:

overhelper <- function(pol, pol.df, fn=sum, verbose=F){
   if(verbose) {cat("Building over geometry...\n"); t=Sys.time(); t}
   geolist = over(geometry(pol), pol.df, returnList = T)
   if(verbose) {cat("Geometry done. Aggregating df. \n"); Sys.time()-t;t=Sys.time();t;}
   results = data.frame(matrix(0,nrow=length(pol), ncol=ncol(pol.df)))
   names(results) = names(pol.df)
   end = length(geolist)

   for (i in 1:end){
     if(verbose) cat(i, "...")
     results[i,] = sapply(pol.df@data[unlist(geolist[i]),], fn)
   }
   if(verbose) cat("Aggregation done! (", Sys.time()-t, ") \n Returning result vector.")
   return (results)
}
Mike Dolan Fliss
la source
1
J'ai ajouté une alternative pour obtenir des rgeos 0.3-8 fixes, à ma réponse.
Edzer Pebesma