如何使用SF在R中进行空间分析

您在哪里投票？ 您是谁的立法者？ 你的邮政编码是多少？ 这些问题在地理空间上有共同点：答案涉及确定点落入哪个多边形。

通常使用专门的GIS软件进行此类计算。 但是在R中也很容易做到。您需要三件事：

1. 对地址进行地址解析以找到经度和纬度的一种方法；
2. 概述邮政编码多边形边界的Shapefile； 和
3. sf软件包。

对于地理编码，我通常使用geocod.io API 。 它每天免费进行2500次查询，并且具有不错的R包 ，但是您需要一个（免费）API密钥才能使用它。 为了解决本文的复杂性，我将使用免费的开源Open Street Map Nominatim API 。 它不需要密钥。 tmaptools程序包具有使用该API的函数geocode_OSM() 。

[也在InfoWorld上： 在我们的“用R做更多”的视频系列中获取Sharon Machlis的R技巧 | 按任务，主题或程序包搜索“使用R进行更多操作”操作视频 ]

导入和准备地理空间数据

我将使用sf，tmaptools，tmap和dplyr软件包。 如果要继续，请使用pacman::p_load()加载每个文件，或使用install.packages()安装尚未安装在系统上的文件，然后使用library()加载每个文件。

在此示例中，我将创建一个带有两个地址的矢量，这是我们位于马萨诸塞州弗雷明汉的IDG办公室和位于波士顿的RStudio办公室。

addresses <- c("492 Old Connecticut Path, Framingham, MA", 
               "250 Northern Ave., Boston, MA")

使用geocode_OSM可以轻松进行地址解析。 您可以通过打印出前三列（包括纬度和经度）来查看结果：

geocoded_addresses <- geocode_OSM(addresses)
print(geocoded_addresses[,1:3])
                                     query      lat       lon
# 1 492 Old Connecticut Path, Framingham, MA 42.31348 -71.39105
# 2            250 Northern Ave., Boston, MA 42.34806 -71.03673

有几种获取邮政编码shapefile的方法。 最简单的可能是美国人口普查局的邮政编码列表区域，该区域与美国邮政服务的边界相似，甚至不完全相同。

您可以直接从美国人口普查局下载ZCTA文件 ，但这是整个国家/地区的文件。 仅当您不介意大数据文件时才这样做。

Census Reporter是下载单个状态的ZCTA文件的地方。 按州搜索任何数据（例如人口），然后将邮政编码添加到地理区域，然后选择下载数据作为shapefile 。

我可以手动解压缩下载的文件，但在R中更容易。在这里，我对下载的文件使用基本R的unzip()函数，并将其解压缩到名为ma_zip_shapefile的项目子目录中。 该junkpaths = TRUE参数表示，我不想基于zip文件的名称对另一个子目录进行解压缩。

unzip("data/acs2017_5yr_B01003_86000US02648.zip", 
      exdir = "ma_zip_shapefile", junkpaths = TRUE, 
      overwrite = TRUE)

使用SF进行地理空间导入和分析

现在终于完成了一些地理空间工作。 我将使用sf的st_read()函数将shapefile导入R。

zipcode_geo <- st_read("ma_zip_shapefile/acs2017_5yr_B01003_86000US02648.shp")
# Reading layer `acs2017_5yr_B01003_86000US02648' from data source `/Users/smachlis/Documents/MoreWithR/ma_zip_shapefile/acs2017_5yr_B01003_86000US02648.shp' using driver `ESRI Shapefile'
# Simple feature collection with 548 features and 4 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -73.50821 ymin: 41.18705 xmax: -69.85886 ymax: 42.95774
# epsg (SRID):    4326
# proj4string:    +proj=longlat +datum=WGS84 +no_defs

运行st_read()时，我包括了控制台响应，因为那里显示了一些信息：epsg。 那就是说使用了什么坐标参考系统来创建文件 。 在这里是4326。在不深入杂草的情况下，epsg基本上可以指示使用什么系统将三维地球（地球）上的区域转换为二维坐标（纬度和经度） 。 这很重要，因为有很多不同的坐标参考系。 我希望我的邮政编码多边形和地址点使用相同的多边形，以便它们正确对齐。

注意：此文件恰好包含了马萨诸塞州整个州的多边形，我不需要。 所以我过滤掉马萨诸塞州的那一行

zipcode_geo <- dplyr::filter(zipcode_geo, 
                             name != "Massachusetts")

使用tmap映射shapefile

映射多边形数据不是必需的，但这是对我的shapefile的很好检查，以查看几何形状是否符合我的期望。 您可以使用tmap的qtm() （快速主题地图的缩写）功能对sf对象进行快速绘图。

qtm(zipcode_geo) +
  tm_legend(show = FALSE)

屏幕由IDG的Sharon Machlis拍摄

带有tmap R软件包的马萨诸塞州邮政编码列表区域的地图。

看起来我确实确实拥有马萨诸塞州的几何图形，而多边形可能是邮政编码。

接下来，我要使用地理编码的地址数据。 当前这是一个普通的数据框，但需要将其转换为具有正确坐标系的sf地理空间对象。

我们可以使用sf的st_as_sf()函数来做到这一点。 （注意：对空间数据进行操作的sf包函数以st_开头， st_代表“ spatial”和“ temporal”。）

[ 通过InfoWorld的机器学习和分析报告时事通讯来了解机器学习，人工智能和大数据分析的最新进展 ]

st_as_sf()接受几个参数。 在下面的代码中，第一个参数是要转换的对象-我的地理编码地址。 第二个自变量向量告诉函数哪些列具有x（经度）和y（纬度）值。 第三个将坐标参考系统设置为4326，因此它与我的邮政编码多边形相同。

point_geo <- st_as_sf(geocoded_addresses, 
                      coords = c(x = "lon", y = "lat"), 
                      crs = 4326)

SF地理空间连接

现在，我已经设置了两个数据集，使用sf的st_join()函数可以轻松计算每个地址的邮政编码。 语法：

st_join(point_sf_object, polygon_sf_object, join = join_type)

在此示例中，我想首先在经过地理编码的点上运行st_join() ，然后在邮政编码区域上运行st_join() 。 这是一种所谓的左联接格式：包括第一个数据（经地理编码的地址）中的所有点，但仅包含第二个（ZIP代码）数据中匹配的点。 最后，我的st_within类型是st_within ，因为我希望匹配点在其中。

my_results <- st_join(point_geo, zipcode_geo, 
                      join = st_within)

而已！ 现在，如果我通过打印几个最重要的列来查看结果，您将看到每个地址都有一个邮政编码（在“名称”列中）。

print(my_results[,c("query", "name", "geometry")])
# Simple feature collection with 2 features and 2 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: -71.39105 ymin: 42.31348 xmax: -71.03673 ymax: 42.34806
# epsg (SRID):    4326
# proj4string:    +proj=longlat +datum=WGS84 +no_defs
#                                      query  name                   geometry
# 1 492 Old Connecticut Path, Framingham, MA 01701 POINT (-71.39105 42.31348)
# 2            250 Northern Ave., Boston, MA 02210 POINT (-71.03673 42.34806)

使用tmap映射点和多边形

如果要映射点和多边形，这是使用tmap的一种方法：

tm_shape(zipcode_geo) +
tm_fill() +
tm_shape(my_results) +
tm_bubbles(col = "red", size = 0.25)

屏幕快照由IDG的Sharon Machlis摄

用tmap映射的马萨诸塞州点和多边形。

是否需要更多R技巧？ 前往“用R做更多”页面！

From: https://www.infoworld.com/article/3505897/how-to-do-spatial-analysis-in-r-with-sf.html