Как пространственный полигон% над% полигоном работает при агрегировании значений в r?

12

Я работаю над проектом по экологической эпидемиологии, где у меня есть точечные воздействия (~ 2000 операций по свиноводству - МГО). Эти МГО распыляют на близлежащие поля, но капли воды и запах фекалий могут преодолевать мили. Таким образом, эти точечные экспозиции получают 3-миллиметровые буферы, и я хочу знать количество экспозиций IHO (различных видов - сумма количества навоза, количество свиней, что угодно; самое простое, просто количество перекрывающихся буферов экспозиции) на блоки переписи NC (~ 200 000). Блоки переписи исключений (синие) - это (1) все, что находится в топ-5 самых густонаселенных городов, и (2) округа, которые не граничат с округом с IHO в нем (примечание: это было сделано с помощью функции gRelate и кодов DE-9IM - очень гладко!) См. Изображение ниже для визуального

введите описание изображения здесь

Последний шаг - агрегирование буферизованного представления рисков для каждого блока переписи. Вот где я в тупике.

До сих пор я хорошо проводил время с функциями% over% в пакете sp, но из овернета я понимаю, что poly-poly и poly-line over реализованы в rgeos. Виньетка покрывает только poly-poly и самоссылающиеся poly, а не агрегацию, поэтому я немного запутался в том, что мои варианты для poly-poly с агрегацией функций, такие как sum или mean.

В качестве тестового примера рассмотрим приведенный ниже несколько подробный фрагмент, работающий с файлом границ стран мира. Это должно быть в состоянии быть скопировано и выполнено как есть, так как я использую случайное начальное число для точек и так как я загружаю и разархивирую файл мира в коде.

Сначала мы создаем 100 точек, затем используем функцию over с аргументом fn для добавления элемента во фрейм данных. Здесь много точек, но взгляните на Австралию: 3 балла, номер 3 в качестве метки. Все идет нормально.

введите описание изображения здесь

Теперь мы трансформируем геометрию, чтобы мы могли создавать буферы, преобразовывать обратно и отображать эти буферы. (Включено в предыдущую карту, так как я ограничен двумя ссылками.) Мы хотим знать, сколько буферов перекрывают каждую страну - в случае Австралии, на глаз, это 4. Я не могу на всю жизнь понять, что происходит хотя, чтобы получить это с помощью функции over. Смотрите мой беспорядок попыток в последних строках кода.

РЕДАКТИРОВАТЬ: Обратите внимание, что комментатор на r-sis-geo упомянул агрегатную функцию - также упоминается в вопросе обмена стека 63577 - так что обход этой функции / поток может быть через эту функцию, но я не понимаю, почему мне нужно идти агрегировать для полиополии, когда кажется, что сверх имеет такую ​​функциональность для других пространственных объектов.

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?
Майк Долан Флисс
источник
Цени перенаправление - я должен удалить отсюда и сделать репост там? Какой лучший ход? Благодарю.
Майк Долан Флисс

Ответы:

5

Спасибо за ясный вопрос и воспроизводимый пример.

Ваше понимание верно, и это сводится к ошибке в rgeos :: over, которая была исправлена месяц назад, но еще не превратилась в релиз CRAN. Ниже приведен обходной вариант, если вас интересует только количество пересечений:

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

Я использую NROWздесь вместо того, lengthчтобы он работал не с rgeos (0.3-8, от CRAN), а с исправленными (0.3-10, от r-forge). Более раннее предложение об использовании

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

также подсчитывает количество пересечений, но только с установленной версией rgeos. Его преимущество, помимо более интуитивного имени, состоит в том, что он напрямую возвращает Spatialобъект с геометрией world.map.

Чтобы заставить работать rgeos 0.3-8, добавьте

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

в ваш сценарий, прежде чем использовать over.

Эдзер Пебесма
источник
Очень полезно, спасибо. Я особенно хочу отметить ваше предложение решения, которое работает до и после исправления. Не могли бы вы остановиться на следующем: (1) В чем заключается ошибка, с которой я сталкиваюсь здесь - rgeos :: over возвращает пространственную полигональную географию, а не пространственный фрейм данных? Разве некоторые функции не возвращают фреймы данных ...? (2) Как обычно предполагается работать с совокупностью и более? Я немного запутался относительно их предполагаемых различий и вариантов использования. Очень ценю ваше взвешивание, спасибо. И sidenote: какие-либо предложения для понимания цикла выпуска CRAN?
Майк Долан Флисс
Кроме того, что касается первоначального вопроса: мне нужно посчитать количество экспозиций, но мне действительно нужно их суммировать - такие вещи, как количество свиней в каждой экспозиции. Подсчет перекрытий - это начало ... но звучит так, как будто мне нужно найти новые rgeos, да? Нет способа сделать это функциональное агрегирование (не просто подсчет) без него?
Майк Долан Флисс
(1) rgeos :: over для подписи SpatialPolygons,SpatialPolygonsDataFrameдолжен возвращать a data.frame, но возвращает индексный вектор, идентичный тому, когда yбыл бы SpatialPolygons. sp::aggregateделает то, что вы делаете, более удобным для пользователя способом, возвращая Spatialобъект вместо data.frame. Пакеты CRAN обслуживаются волонтерами.
Эдзер Пебесма
ОК, спасибо Эдзер. Похоже, что агрегат опирается на rgeos, поэтому для того, чтобы вывести эту функциональность перед циклом выпуска CRAN (где бы он ни был), мне нужно выяснить, как загрузить новейшие rgeo и как с этим работать. Спасибо. И спасибо за всю вашу работу над пакетом !!
Майк Долан Флисс
Кроме того, Эдзер, большое спасибо за заметку о R-sis-geo. Я не был уверен, где лучшее место для публикации, поэтому я рад, что тема теперь указывает здесь.
Майк Долан Флисс
1

Тем временем я подхватил быстрый (и плохо закодированный) заменитель, который создает нужный мне фрейм данных, поскольку на мой вопрос не совсем ответили вышеупомянутое решение только для подсчета или «отработать новые rgeos», которые я Я не достаточно опытен, чтобы понять, как это сделать.

Эта функция явно (1) неполна (обратите внимание, как я игнорирую аргумент fn) и (2) неэффективна, так как я иду на нее без мощных манипуляций с массивом R / sapply ... (ясно, что я прихожу из других языков без эта сила) но, честно говоря, я все еще не понимаю, что возвращает структура функции over (список списков ...? И пустые списки, если нет данных?). Что бы это ни стоило (редактирование приветствуется), эта функция успешно выполняет ту работу, которая мне нужна, и имитирует действие других над функциями.

Редактирование приветствуется:

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)
}
Майк Долан Флисс
источник
1
Я добавил альтернативу, чтобы исправить rgeos 0.3-8, к моему ответу.
Эдзер Пебесма