создать растр плотности и извлечь сумму с помощью полигона

Question

У меня есть многоугольник (zones) и набор координат (points). Я хотел бы создать пространственный растр плотности ядра для всего многоугольника и извлечь сумму плотности по зонам. Точки за пределами многоугольника следует отбрасывать.

library(raster)
library(tidyverse)
library(sf)
library(spatstat)
library(maptools)

load(url("https://www.dropbox.com/s/iv1s5butsx2v01r/example.RData?dl=1"))
# alternatively, links to gists for each object
# https://gist.github.com/ericpgreen/d80665d22dfa1c05607e75b8d2163b84
# https://gist.github.com/ericpgreen/7f4d3cee3eb5efed5486f7f713306e96

ggplot() +
  geom_sf(data = zones) +
  geom_sf(data = points) +
  theme_minimal()

Я попытался преобразовать в ppp с помощью {spatstat} и затем использовать density(), но меня смущают единицы в результате. Я считаю, что проблема связана с единицами карты, но я не знаю, как действовать дальше.

Обновление

Вот код для воспроизведения карты плотности, которую я создал:

zones_owin <- as.owin(as_Spatial(zones))
pts <- st_coordinates(points)
p <- ppp(pts[,1], pts[,2], window=zones_owin, unitname=c("metre","metres"))
ds <- density(p) 
r <- raster(ds)
plot(r)

Answer 1

Единицы сложны, когда вы работаете напрямую с географическими c координатами (долгота, широта). Если возможно, вы должны преобразовать в плоские координаты (что является обязательным требованием для spatstat) и перейти оттуда. Плоские координаты обычно бывают в единицах метров, но я думаю, что это зависит от заданной проекции c и лежащего в основе эллипсоида et c. Вы можете увидеть этот ответ о том, как проецировать на плоские координаты с помощью sf и экспортировать в формат spatstat с использованием maptools. Примечание: вам нужно вручную выбрать разумную проекцию (вы можете использовать http://epsg.io, чтобы найти ее), и вы должны спроецировать как многоугольник, так и точки.

Когда все в формате spatstat, вы можете использовать density.ppp для сглаживания ядра. Результирующие значения сетки (объект класса im) представляют собой интенсивность точек, то есть количество точек на квадратную единицу (например, квадратный метр). Если вы хотите агрегировать по некоторому региону, вы можете использовать integral.im(..., domain = ...), чтобы получить ожидаемое количество точек в этом регионе для модели точечного процесса с заданной интенсивностью.

Answer 2

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

Удаляет все точки, не входящие в один из полигонов 'зоны', считает их по зонам и наносит на график зоны, окрашенные количеством точек, попадающих в зону.

library(raster)
library(tidyverse)
library(sf)
#> Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3
library(spatstat)

library(maptools)
#> Checking rgeos availability: TRUE

load(url("https://www.dropbox.com/s/iv1s5butsx2v01r/example.RData?dl=1"))
# alternatively, links to gists for each object
# https://gist.github.com/ericpgreen/d80665d22dfa1c05607e75b8d2163b84
# https://gist.github.com/ericpgreen/7f4d3cee3eb5efed5486f7f713306e96

p1 <- ggplot() +
  geom_sf(data = zones) +
  geom_sf(data = points) +
  theme_minimal()

#Remove points outside of zones
points_inside <- st_intersection(points, zones)
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
nrow(points)
#> [1] 308
nrow(points_inside)
#> [1] 201

p2 <- ggplot() + 
  geom_sf(data = zones) + 
  geom_sf(data = points_inside)

points_per_zone <- st_join(zones, points_inside) %>%
  count(LocationID.x)
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar

p3 <- ggplot() + 
  geom_sf(data = points_per_zone, 
          aes(fill = n)) + 
  scale_fill_viridis_c(option = 'C')

points_per_zone
#> Simple feature collection with 4 features and 2 fields
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 34.0401 ymin: -1.076718 xmax: 34.17818 ymax: -0.9755066
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +ellps=WGS84 +no_defs
#> # A tibble: 4 x 3
#>   LocationID.x     n                                                    geometry
#> *        <dbl> <int>                                               <POLYGON [°]>
#> 1           10   129 ((34.08018 -0.9755066, 34.0803 -0.9757393, 34.08046 -0.975…
#> 2           20    19 ((34.05622 -0.9959458, 34.05642 -0.9960835, 34.05665 -0.99…
#> 3           30    29 ((34.12994 -1.026372, 34.12994 -1.026512, 34.12988 -1.0266…
#> 4           40    24 ((34.11962 -1.001829, 34.11956 -1.002018, 34.11966 -1.0020…

cowplot::plot_grid(p1, p2, p3, nrow = 2, ncol = 2)

Кажется, я недооценил сложность вашей проблемы. Является ли что-то похожее на график ниже (и лежащие в его основе данные), что вы ищете?

Используется растр с сеткой ~ 50x50, растр :: фокус с окном 9x9, используя среднее значение для интерполяции данных.