创建密度栅格并按多边形特征提取总和
create density raster and extract sum by polygon feature
我有一个多边形 (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()
我尝试使用 {spatstat
} 转换为 ppp
,然后使用 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)
我不确定这是否回答了您的所有问题,但应该是一个好的开始。如果您需要不同类型的输出,请在评论或问题中澄清。
它会删除不在 'zone' 多边形之一内的所有点,按区域对它们进行计数,并绘制按落入其中的点数着色的区域。
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 网格的栅格,raster::focal 具有 9x9 的 window 使用均值对数据进行插值。
当您直接使用地理坐标(经度、纬度)时,很难找到单位。如果可能,您应该转换为平面坐标(这是 spatstat
的要求)并从那里继续。平面坐标通常以米为单位,但我想这取决于特定的投影和底层椭球等。你可以看到 for how to project to planar coordinates with sf
and export to spatstat
format using maptools
. Note: You have to manually choose a sensible projection (you can use http://epsg.io 找到一个)并且你必须同时投影多边形和点。
一旦所有内容都采用 spatstat
格式,您就可以使用 density.ppp
进行内核平滑。生成的网格值(class im
的对象)是点的强度,即每平方单位(例如平方米)的点数。如果您想在某个区域进行聚合,您可以使用 integral.im(..., domain = ...)
来获取具有给定强度的点过程模型在该区域中的预期点数。
我有一个多边形 (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()
我尝试使用 {spatstat
} 转换为 ppp
,然后使用 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)
我不确定这是否回答了您的所有问题,但应该是一个好的开始。如果您需要不同类型的输出,请在评论或问题中澄清。
它会删除不在 'zone' 多边形之一内的所有点,按区域对它们进行计数,并绘制按落入其中的点数着色的区域。
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 网格的栅格,raster::focal 具有 9x9 的 window 使用均值对数据进行插值。
当您直接使用地理坐标(经度、纬度)时,很难找到单位。如果可能,您应该转换为平面坐标(这是 spatstat
的要求)并从那里继续。平面坐标通常以米为单位,但我想这取决于特定的投影和底层椭球等。你可以看到 sf
and export to spatstat
format using maptools
. Note: You have to manually choose a sensible projection (you can use http://epsg.io 找到一个)并且你必须同时投影多边形和点。
一旦所有内容都采用 spatstat
格式,您就可以使用 density.ppp
进行内核平滑。生成的网格值(class im
的对象)是点的强度,即每平方单位(例如平方米)的点数。如果您想在某个区域进行聚合,您可以使用 integral.im(..., domain = ...)
来获取具有给定强度的点过程模型在该区域中的预期点数。