Подсчет точек внутри полигонов по коэффициенту в R

Question

У меня есть два набора данных координат xy. Первая имеет координаты xy плюс столбец тегов с моими уровнями факторов. Я назвал data.frameqq и это выглядит следующим образом:

structure(list(x = c(5109, 5128, 5137, 5185, 5258, 5324, 5387, 
5343, 5331, 5347, 5300, 5180, 4109, 4082, 4091, 4139, 4212, 4279, 
4291, 4297, 4285, 4301, 4254, 4181), y = c(1692, 1881, 2070, 
2119, 2144, 2065, 1987, 1813, 1705, 1649, 1631, 1654, 1847, 2015, 
2204, 2253, 2278, 2282, 2166, 1947, 1839, 1783, 1765, 1783), 
    tag = c("MPN_right", "MPN_right", "MPN_right", "MPN_right", 
    "MPN_right", "MPN_right", "MPN_right", "MPN_right", "MPN_right", 
    "MPN_right", "MPN_right", "MPN_right", "MPN_left", "MPN_left", 
    "MPN_left", "MPN_left", "MPN_left", "MPN_left", "MPN_left", 
    "MPN_left", "MPN_left", "MPN_left", "MPN_left", "MPN_left" 
    )), .Names = c("x", "y", "tag"), row.names = c(NA, -24L), class = "data.frame") 

я сгенерировал случайные данные для другого, используя qqxy средство с большим sd.

set.seed(123) 
my_points=data.frame(x=rnorm(n =1000,mean=mean(qq$x),sd=1000), 
y=rnorm(n=1000,mean=mean(qq$y),sd=1000)) 

, если я использую функцию in.out из mgcv пакета, я получаю несколько, что я хочу.

Основные проблемы, связанные с этим подходом, состоят в том, что мой «многоугольник» не закрыт и не будет интерпретироваться как 2 многоугольника. Пакет советует использовать одну строку NA между ними, но я бы предпочел использовать столбец тегов, так как я буду пытаться использовать более двух уровней в моем теге, а именно более двух полигонов). Моя конечная цель - создать таблицу количества точек в каждом из них.

Answer 1

lapply() и sapply() помогут вам использовать уровень функции par.

## a bit edited to make output clear 

library(dplyr); library(mgcv) 

TAG <- unique(qq$tag) 

IN.OUT <- lapply(TAG, function(x) as.matrix(qq[qq$tag==x, 1:2])) %>% # make a matrix par level 
    sapply(function(x) in.out(x, as.matrix(my_points)))  # use in.out() with each matrix 

colnames(IN.OUT) <- TAG 

head(IN.OUT, n = 3) 

#  MPN_right MPN_left 
# [1,]  FALSE FALSE 
# [2,]  FALSE FALSE 
# [3,]  FALSE FALSE 

apply(IN.OUT, 2, table) 

#  MPN_right MPN_left 
# FALSE  983  990 
# TRUE   17  10 

Answer 2

что по этому поводу:

mysppoint <- SpatialPoints(coords = my_points) # create spatial points 
qq$tag <- as.factor(qq$tag) 
polys = list() 

# create one polygon for each factor level 
for (lev in levels(qq$tag)){ 
    first_x <- qq$x[qq$tag == lev][1] 
    first_y <- qq$y[qq$tag == lev][1] 
    qq <- rbind(qq, data.frame(x = first_x, y = first_y, tag = lev)) # "close" the polygon by replicating the first row 
    polys[[lev]] <- Polygons(list(Polygon(matrix(data = cbind(qq$x[qq$tag == lev], # transform to polygon 
                  qq$y[qq$tag == lev]), 
               ncol = 2))), lev) 
} 

mypolys <- SpatialPolygons(polys) # convert to spatial polygons 
inters <- factor(over(mysppoint, mypolys), labels = names(mypolys)) # intersect points with polygons 
table(inters) 

, что дает:

inters 
MPN_left MPN_right 
     10  17 

Преимуществом этого является то, что она дает правильные пространственные объекты для работы. Например:

plotd <- fortify(mypolys) 
p <- ggplot() 
p <- p + geom_point(data = my_points, aes(x = x , y = y), size = 0.2) 
p <- p + geom_polygon(data = plotd, aes(x = long, y = lat, fill = id), alpha = 0.7) 
p 

Answer 3

Я закончил с использованием lapply и сочетание split и более lapply. Итак, вот код, пожалуйста, игнорируйте вспомогательную функцию extract_coords, она в основном дает мне dataframe с x, y и столбцами тегов. Мне также удалось подобрать точки от оригинала your_coords и посчитать их (возвращая их как вектор вместо таблицы).

inside_ROI = function(your_ROI_zip,your_coords){ 

    # Helper function will take list from zip ROIs and merge them into a df 
    qq=extract_coords(your_ROI_zip) 

    # We use tag for splitting by the region 

    lista=split(qq,qq$tag) 

    # We check if they are in or out 
    who_is_in = lapply(lista,function(t) in.out(cbind(t$x,t$y),x=cbind(your_coords$x,your_coords$y))) 

    # We sum to get the by area countings 
    region_sums = unlist(lapply(who_is_in,function(t) sum(as.numeric(t)))) 

    # obtain indices for subset of TRUE values 
    aa=lapply(who_is_in,function(p) which(p==T)) 

    whos_coords=list() 

for (i in aa){ 
    whos_coords = append(whos_coords,values=list(your_coords[i,])) 
    # whos_coords[[i]] = your_coords[i,] 
    } 

    # Change names 
    names(whos_coords) = names(aa) 

    # Put into list for more than one output 
    out=list(region_sums,aa,whos_coords) 

    return(out) 
}