Определение того, находится ли растр внутри объекта SpatialPolygons, без него или пересекает его

У меня есть много растров, для которых я хотел бы проверить, полностью ли они содержатся в пространственном многоугольнике, полностью без пространственного многоугольника или пересекаются с пространственным многоугольником (это может означать, что многоугольник полностью находится в растре или многоугольник и растровое перекрытие). Я делаю эту проверку, чтобы сэкономить время, когда это возможно.

Вот пример:

                                      # create 3 example rasters
r <- raster()
r[] <- rnorm(n = ncell(r))
e1 <- extent(c(45,55,45,50))
r1 <- crop(r,e1)
e2 <- extent(c(20,25,25,30))
r2 <- crop(r,e2)
e3 <- extent(c(38,55,57,65))
r3 <- crop(r,e3)


                                        #create SpatialPolygons

x <- c(40,60)
y <- c(40,60)
m <- expand.grid(x,y)
m <- m[c(1,2,4,3),]
p1 <- Polygon(m)
p1 <- Polygons(list(p1),1)

x <- c(10,15)
y <- c(10,15)
m <- expand.grid(x,y)
m <- m[c(1,2,4,3),]
p2 <- Polygon(m)
p2 <- Polygons(list(p2),2)

x <- c(30,45)
y <- c(70,80)
m <- expand.grid(x,y)
m <- m[c(1,2,4,3),]
p3 <- Polygon(m)
p3 <- Polygons(list(p3),3)


poly <- SpatialPolygons(list(p1,p2,p3))

построение этих:

три растра и три многоугольника

Я буду читать каждый растр отдельно и проверять, находится ли он внутри, вне или пересекается с пространственными полигонами.

Как вы думаете, какой способ сделать это в R будет наиболее эффективным? У меня есть тысячи растров размером 4 Мб, которые я планирую замаскировать параллельно, и мне хотелось бы, чтобы эта проверка немного ускорила процесс.

Обратите внимание, есть еще этот вопрос: https://gis.stackexchange.com/questions/34535/detect-whether-there-is-a-spatial-polygon-in-a-spatial-extent

Однако я не думаю, что это дает то, что я ищу. Например, все растры находятся в пределах экстента пространственных многоугольников, но не все находятся внутри пространственных многоугольников.

Функции, подобные тем, что есть в rgeos (gIntersects, gContains), вероятно, были бы удобны. Я не уверен, являются ли они наиболее эффективными или как мне преобразовать растр (или его размер) в объект sp.

Благодарность!


person Tedward    schedule 08.01.2016    source источник


Ответы (2)


Вы также можете использовать для этого gRelate. Он возвращает код DE-9IM, который описывает отношения между внутренними, граничными и внешними компонентами. двух геометрий.

library(rgeos)
x <- sapply(rlist, function(x) 
  gsub('[^F]', 'T', gRelate(as(extent(x), 'SpatialPolygons'), poly)))

Затем вы можете сравнить строки с интересующими вас отношениями. Например, мы можем определить within, disjoint и overlaps следующим образом (но обратите внимание, что некоторые другие пересечения являются необязательными для данных отношений - «внутри» - это определен GEOS как T*F**F***," непересекающийся "как FF*FF**** и" перекрывается "как T*T***T**):

pat <- c(TFFTFFTTT='within', FFTFFTTTT='disjoint', TTTTTTTTT='overlaps')
pat[x]

##  TFFTFFTTT  FFTFFTTTT  TTTTTTTTT 
##   "within" "disjoint"  "overlaps" 

Кажется, что он немного быстрее, чем подход _10 _ / _ 11_, но @ Tedward post более понятен и больше соответствует определениям GEOS (хотя возможность создания конкретных определений взаимосвязей может быть желательной).


Элементы строк DE-9IM представляют по порядку:

  1. Пересекает ли внутренность геометрии A внутренность геометрии B?
  2. Пересекает ли граница A внутренность B?
  3. Пересекает ли внешность A внутренность B?
  4. Пересекает ли внутренность геометрии A границу геометрии B?
  5. Пересекает ли граница A границу B?
  6. Пересекает ли внешность A границу B?
  7. Пересекает ли внутренность геометрии A внешность геометрии B?
  8. Пересекает ли граница A внешность B?
  9. Пересекается ли внешность A с внешней стороной B?
person jbaums    schedule 09.01.2016
comment
Спасибо за объяснение gRelate. Гибкость очень полезна, и о ней полезно знать. - person Tedward; 09.01.2016

Вот что я сделал для решения проблемы:

library(rgeos)
rlist <- list(r1,r2,r3)

lapply(rlist, function(raster) {
  ei <- as(extent(raster), "SpatialPolygons")
  if (gContainsProperly(poly, ei)) {
    print ("fully within")
  } else if (gIntersects(poly, ei)) {
    print ("intersects")
  } else {
    print ("fully without")
  }
})

Пожалуйста, дайте мне знать, если вы знаете более эффективное решение.

person Tedward    schedule 08.01.2016
comment
Вы сэкономите немного времени, удалив вызовы print и назначение proj4string. - person jbaums; 09.01.2016
comment
Спасибо, звонки в печать - просто заполнители. Я собираюсь вызвать некоторые другие функции. Хороший звонок по заданию proj4string. - person Tedward; 09.01.2016
comment
(Я имел в виду, что вы можете просто иметь, например, 'fully within', а не print('fully within'), но если они просто заполнители, не беспокойтесь об этом) - person jbaums; 09.01.2016
comment
Попался, это приятно знать. Большое спасибо! - person Tedward; 09.01.2016