映射R中每个位置的半径内有多少点答案

【问题标题】：Mapping how many points are within a radius of every location in R映射R中每个位置的半径内有多少点
【发布时间】：2021-03-18 09:50:34
【问题描述】：

在 R 中，我正在尝试创建一个等值线图。我已经建立了一个企业数据库，有些是连锁店的一部分（例如麦当劳），有些是独立的。我想计算地图上每个点 30 公里范围内有多少商家，但将连锁店的不同位置视为一个商家。

例如，如果你的一个点是：

距离麦当劳 5 公里，
距塔可钟 10 公里
距离 Chick-Fil-A 15 公里
离肯德基20公里
距离麦当劳 25 公里
距离五个人 35 公里

颜色表示30公里内有4家快餐店。

我很乐意使用任何 R 包，但我最熟悉 tmaps 和 ggplot2 地图。

在这个阶段，我能想到的最好方法是为每个链创建多边形并将它们堆叠为相同颜色的透明层。我认为这不会非常有效，也不会创建一个非常漂亮的等值线。

我能找到的其他答案是计数点（例如 https://gis.stackexchange.com/questions/229066/counting-how-many-times-a-point-is-inside-a-set-of-intersecting-polygons-in-r）或 GIS 软件。

编辑： 我设法从每条链的每个位置创建了一个 30 公里的半径（使用 rgeos gIntersection）。我现在有一系列多边形。

为了解决我的问题，我需要做的额外事情是为 where 创建多边形：

只有一个多边形覆盖该区域，
两个多边形覆盖该区域，
等

为了尝试视觉化，我使用了https://gis.stackexchange.com/questions/229066/counting-how-many-times-a-point-is-inside-a-set-of-intersecting-polygons-in-r的答案

在链接的问题中，他们试图计算有多少多边形覆盖了编号点（右图）。我要做的是在左侧创建图像，其中有不重叠的多边形（1），两个重叠的多边形（2）等等。

【问题讨论】：

只是一个技术问题 - 您要在周围绘制多边形的点是否与商家的位置不同？
嗨贾斯汀，我基本上想要一张新西兰地图，根据每个点 30 公里内有多少快餐店进行着色。我认为实际上这意味着在每个链中围绕每个业务建立半径（我可以做到）。然后每次有不同的链重叠时改变颜色 - 这个我还没想通。我认为可能有更简洁的方法。
从评论中听起来您想创建一个网格并将每个单元格中的远食餐厅制成表格。光栅方法在这里会更合适。
嗨，Brian，你能告诉我在哪里可以找到有关如何使用 Raster 方法的教程或示例吗？
Robin Lovelace 在 R 中的地理计算是一个很好的参考 geocompr.robinlovelace.net/spatial-class.html#raster-data

标签： r ggplot2 choropleth tmap

【解决方案1】：

我认为最好使用光栅方法而不是叶绿素来完成您要完成的工作。要制作 chorlopleth，您需要定义一组（通常是不规则的）多边形，总结每个多边形内的某些内容，然后根据属性为多边形着色。如果您想知道每个州或县有多少家快餐店，或者每个州人均快餐店的数量，这将是一个好方法。
但是，根据您的描述，您正在寻找所有点的设定半径内有多少快餐店。这更像是一个栅格问题，因为您可以在常规网格上表示您的数据。
raster 包是处理栅格数据的良好开端，并且与sf 包配合得很好。

您需要确定实现目标所需的密度，然后使用它来确定栅格的分辨率。一旦你得到了，你可以使用 raster::rasterize() 来总结你的（我假设）点数据。

我假设你有一个包含每家餐厅位置的对象，我将把这个对象称为“点”。

library(raster)
library(sf)

# create raster template with 30km resolution (I'm assuming your projection is in meters)
raster_template = raster((extent(points), 
                         resolution = 30000,
                         crs = st_crs(points)
                         )
# rasterize your point data
r = rasterize(points, raster_template, fun = "count")

这应该创建一个网格，其中每个单元格都有每个 30 公里单元格内的点数。然后您应该能够绘制栅格，但可能需要对其进行剪辑或遮罩以仅显示新西兰境内的部分

【讨论】：