Блог

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

Для начала у меня есть SpatialPointsDataFrame с несколькими переменными, измеренными для каждой точки. Я хотел бы получить среднее значение всех переменных для точек в пределах указанного расстояния. Например, получение средних значений кадмия из meuse данных для точек в пределах 100 м друг от друга:

 library(sf)
library(sp)
data(meuse)
pts <- st_as_sf(meuse, coords = c("x", "y"),remove=FALSE)
pts100 <- st_is_within_distance(pts, dist = 100)
# can use sapply to get mean of a variable. E.g., cadmium
sapply(pts100, function(x){ mean(pts$cadmium[x]) })
  

Итак, я выяснил, как использовать sapply для выполнения этого от переменной к переменной. Таким образом, я мог бы, если бы захотел, вычислить среднее значение для каждой переменной, сгенерировать центроид для каждой точки, а затем SpatialPointsDataFrame уникальных значений. Например, для первых нескольких переменных:

 res <- data.frame(id=1:length(pts100),
                  x=NA, y=NA,
                  cadmium=NA, copper=NA, lead=NA)
res$x <- sapply(pts100, function(p){ mean(pts$x[p]) })
res$y <- sapply(pts100, function(p){ mean(pts$y[p]) })
res$cadmium <- sapply(pts100, function(p){ mean(pts$cadmium[p]) })
res$copper <- sapply(pts100, function(p){ mean(pts$copper[p]) })
res$lead <- sapply(pts100, function(p){ mean(pts$lead[p]) })
res2 <- res[duplicated(res$cadmium),]
coordinates(res2) <- c("x","y")
bubble(res2,"cadmium")
  

Это работает, но кажется громоздким, и, похоже, должен быть более эффективный способ.

Похоже, что для sf пакета есть агрегатная функция, которая имеет аргумент join, в котором вы можете указать тип соединения.

 ibrary(sf)
library(sp)
data(meuse)
pts <- st_as_sf(meuse, coords = c("x", "y"),remove=FALSE)

# This will give lots of warnings since there are non-numeric columns
pts_agg <- aggregate(pts,
                     pts,
                     FUN = mean, 
                     join = function(x, y) st_is_within_distance(x, y, dist = 100))

head(pts_agg)

Simple feature collection with 6 features and 14 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 181025 ymin: 333260 xmax: 181390 ymax: 333611
CRS:            NA
         x        y cadmium copper lead   zinc  elev        dist   om ffreq soil lime landuse dist.m
1 181048.5 333584.5   10.15     83  288 1081.5 7.446 0.006791165 13.8    NA   NA   NA      NA     40
2 181048.5 333584.5   10.15     83  288 1081.5 7.446 0.006791165 13.8    NA   NA   NA      NA     40
3 181165.0 333537.0    6.50     68  199  640.0 7.800 0.103029000 13.0    NA   NA   NA      NA    150
4 181298.0 333484.0    2.60     81  116  257.0 7.655 0.190094000  8.0    NA   NA   NA      NA    270
5 181307.0 333330.0    2.80     48  117  269.0 7.480 0.277090000  8.7    NA   NA   NA      NA    380
6 181390.0 333260.0    3.00     61  137  281.0 7.791 0.364067000  7.8    NA   NA   NA      NA    470
               geometry
1 POINT (181072 333611)
2 POINT (181025 333558)
3 POINT (181165 333537)
4 POINT (181298 333484)
5 POINT (181307 333330)
6 POINT (181390 333260)
  

Выборочная проверка 9-й строки pts, так как в pts100 было несколько совпадений:

 > pts[pts100[[9]], 'cadmium'] %>% st_drop_geometry %>% summarise(mean = mean(cadmium))
  mean
1 2.25


> pts_agg[9,'cadmium']
Simple feature collection with 1 feature and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 181060 ymin: 333231 xmax: 181060 ymax: 333231
CRS:            NA
  cadmium              geometry
9    2.25 POINT (181060 333231)