R中的批量地理编码街道交叉口答案

【问题标题】：Batch geocoding street intersections in RR中的批量地理编码街道交叉口
【发布时间】：2022-10-24 03:19:40
【问题描述】：

我目前正在处理格式如下的数据：

    tribble(
      ~street1, ~street2, ~county, ~state
      N BENTON WY, W TEMPLE ST, LOS ANGELES, CA,
      11TH PL, BLAINE ST, LOS ANGELES, CA,
      W 6TH ST, HOPE ST, LOS ANGELES, CA,
      S GRAND AV, W 18TH ST, LOS ANGELES, CA,
      BROADWAY, 5TH ST, LOS ANGELES, CA,
    )

这对应于包含大约 825,000 个缺少坐标的观测值的数据集。这些数据只有最近的十字路口的名称、县和州信息（注意它们不包括街道编号）。我需要对这些观察结果进行地理编码并恢复坐标，以便我的最终数据看起来像这样：

   tribble(
     ~street1, ~street2, ~county, ~state, ~latitude, ~longitude
     N BENTON WY, W TEMPLE ST, LOS ANGELES, CA, XX.XXXX, -YY.YYYY,
     11TH PL, BLAINE ST, LOS ANGELES, CA, XX.XXXX, -YY.YYYY,
     W 6TH ST, HOPE ST, LOS ANGELES, CA, XX.XXXX, -YY.YYYY,
     S GRAND AV, W 18TH ST, LOS ANGELES, CA, XX.XXXX, -YY.YYYY,
     BROADWAY, 5TH ST, LOS ANGELES, CA, XX.XXXX, -YY.YYYY,
   )

我已经研究了一些可能的解决方案，但还没有找到可行的方法。

虽然 Google Maps API（ggmap 包）非常擅长识别来自十字路口的坐标作为输入，但对这么多观察进行地理编码的成本（根据他们的website，每 1000 次查询 4.00 美元）使得该选项不可行。

我查看了其他软件包的文档，例如 RDSTK 和 tidygeocoder，但它们似乎不支持使用两个街道名称作为输入的 API 查询。人口普查地理编码器同样没有该选项，只允许单个地址输入。

在阅读了this 非常详细的 StackOverflow 答案之后，通过osmdata 包使用 OpenStreetMap API 似乎是一个有前途的选择，但是尝试使用更大的边界框复制此代码每次都会产生运行时错误。

例如，请参阅以下代码，使用洛杉矶县，遵循上述帖子中用户 Hugh-allan 的格式：

library(sf)
library(tidyverse)
library(osmdata)

tribble(
      ~point, ~lat, ~lon, 
      1, 32.75004, -118.951721, 
      2, 34.823302, -118.951721, 
      3, 34.823302, -117.646374, 
      4, 32.75004, -117.646374,
    ) %>% 
      st_as_sf(
        coords = c('lon', 'lat'), 
        crs = 4326
      ) %>% 
      {. ->> LA_bounds}
    
    st_bbox(LA_bounds) %>% 
      opq %>% 
      add_osm_feature(key = 'highway') %>% 
      osmdata_sf %>% 
      `[[`('osm_lines') %>% 
      {. ->> LA_streets}

如果有人知道如何使用 OpenStreetMaps 解决此错误，或者以其他方式调整另一个包的语法以适应交叉街道和县作为输入，我将不胜感激。

【问题讨论】：

标签： r google-maps geospatial openstreetmap geocode

【解决方案1】：

我没有 osmdata 的解决方案。但是，我确实在tidygeocoder 上尝试过。如果你正在寻找批在不需要 API 密钥的情况下进行编码，唯一免费的方法是美国人口普查局tidygeocoder，但计算成本很高。为此，我将street1 和street2 与与符号& 结合在一起。然后将其与 county 和 state 组合成一个名为 line_address 的单列，而不是多列：

examples_address <- tibble(line_address= c("N BENTON WY & W TEMPLE ST, LOS ANGELES, CA", "11TH PL & BLAINE ST, LOS ANGELES, CA", "W 6TH ST & HOPE ST, LOS ANGELES, CA", "S GRAND AV & W 18TH ST, LOS ANGELES, CA", "BROADWAY & 5TH ST, LOS ANGELES, CA"))
examples_address1 <- examples_address %>% 
    tidygeocoder::geocode(address = line_address, method = "census", verbose = TRUE)
examples_address1

我得到的输出：

line_address lat long
N BENTON WY & W TEMPLE ST, LOS ANGELES, CA  34.07289 -118.2757
11TH PL & BLAINE ST, LOS ANGELES, CA    NA  NA
W 6TH ST & HOPE ST, LOS ANGELES, CA  34.04944 -118.2563
S GRAND AV & W 18TH ST, LOS ANGELES, CA  34.03420 -118.2673
BROADWAY & 5TH ST, LOS ANGELES, CA  34.04808 -118.2507

不幸的是，正如你在上面看到的，并不是所有的行都给了我们一个纬度和长从批量查询返回。

我们可以在函数内部使用method = "argis" 来为我们提供所有结果，但由于某些原因，返回的纬度和长可能不同。见最后一条：

line_address lat long
N BENTON WY & W TEMPLE ST, LOS ANGELES, CA   34.07290 -118.2757
11TH PL & BLAINE ST, LOS ANGELES, CA     34.04517 -118.2716
W 6TH ST & HOPE ST, LOS ANGELES, CA  34.04946 -118.2564
S GRAND AV & W 18TH ST, LOS ANGELES, CA  34.03417 -118.2673
BROADWAY & 5TH ST, LOS ANGELES, CA   34.01587 -118.4927

arcgis 不支持批查询tidygeocoder。

【讨论】：