如何有效地将纬度/经度对映射到邮政编码（无需 API）答案

【问题标题】：How to efficiently map lat/long pairs to zipcodes (without an API)如何有效地将纬度/经度对映射到邮政编码（无需 API）
【发布时间】：2015-08-11 17:12:47
【问题描述】：

我正在将大量纬度/经度坐标对映射到相关的邮政编码。由于调用限制，记录数量太大而无法使用 Google Maps 或 geonames 等 API。

我有一个查找表，其中包含邮政编码和每个邮政编码的纬度/经度质心。您可以在此处获取查找表：

# zipcode data with lat/lon coordinates
url <- "http://www.boutell.com/zipcodes/zipcode.zip"
fil <- "ziplatlong.zip"

# download an unzip
if (!file.exists(fil)) { download.file(url, fil) }
unzip(fil, exdir="zips")

library(readr)
ziplkp<-read_csv("zips/zipcode.csv")

对于我数据中的每个纬度/经度对，我想通过查找该纬度/经度对与查找表中每个质心之间的最小绝对差来将其与最近的邮政编码质心匹配。

将这种“查找”函数逐行应用于大量记录的最有效方法是什么？

示例数据：经纬度坐标列表：

latlongdata <- 
  structure(list(dropoff_longitude = c(-73.981705, -73.993553, 
-73.973305, -73.988823, -73.938484, -74.015503, -73.95472, -73.9571, 
-73.971298, -73.99794), dropoff_latitude = c(40.760559, 40.756348, 
40.762646, 40.777504, 40.684692, 40.709881, 40.783371, 40.776657, 
40.752148, 40.720535)), row.names = c(8807218L, 9760613L, 3175671L, 
10878727L, 2025038L, 5345659L, 14474481L, 1650223L, 684883L, 
9129975L), class = "data.frame", .Names = c("dropoff_longitude", 
"dropoff_latitude"))

    print(latlongdata)
         dropoff_longitude dropoff_latitude
8807218          -73.98171         40.76056
9760613          -73.99355         40.75635
3175671          -73.97330         40.76265
10878727         -73.98882         40.77750
2025038          -73.93848         40.68469
5345659          -74.01550         40.70988
14474481         -73.95472         40.78337
1650223          -73.95710         40.77666
684883           -73.97130         40.75215
9129975          -73.99794         40.72053

**ZipLooker 函数：查找从输入坐标对到最近的邮政编码质心的最小绝对距离并返回该邮政编码

library(dplyr)
ZipLooker<-function(dropoff_longitude,dropoff_latitude){
  if(is.na(dropoff_longitude)|is.na(dropoff_latitude)){
    z<-NA_character_
  } else {

    tryCatch({
      x<-ziplkp1 
      x$latdiff=abs(dropoff_latitude-x$Latitude)
      x$londiff=abs(dropoff_longitude-x$Longitude)
      x$totdiff=x$latdiff+x$londiff
      z<-head(top_n(x,1,-totdiff),n=1)$Postal
      return(z)
    }, error=function(e) NA)
  }
}

使用 dplyr 的 rowwsie() 函数应用 Ziplooker 函数

  latlongdata %>% 
  rowwise() %>% 
  mutate(zipcode=ZipLooker(dropoff_longitude,dropoff_latitude)
         )

【问题讨论】：

如果你是学生，我推荐 Smarty Streets。他们为学生/学者提供免费帐户，没有通话限制。
谢谢，迈克尔。不是学生，现在每月 1,000 美元稍微超出预算！ ;) 此外，Smarty Streets 看起来对于在给定邮政编码的情况下查找纬度/经度坐标很有用，但我正在尝试另一种方式
就ZipLooker 和mutate 遇到的问题而言：在这种情况下，如果if 没有else，如果您在你的dropoff 变量之一，所以用else 包裹整个TryCatch。此外，由于您最终会在此处返回一个字符变量，因此在您的 if 语句中使用 "NA" 或 NA_character_ 会有所帮助。使用ZipLooker 时不要忘记定义ziplkp 参数，否则您将获得所有NA 值。
感谢 aosmith，所有伟大的建议和更正。 ZipLooker 现在似乎可以工作了，所以我将问题改成了更具体的效率问题。

标签： r dplyr

【解决方案1】：

这里有一个完整的解决方案：

library(sp)
library(maptools)
library(zipcode)

# grab the zip code boundaries
url <- "http://www2.census.gov/geo/tiger/GENZ2014/shp/cb_2014_us_zcta510_500k.zip"
fil <- "ztca.zip"

# don't waste bandwidth
if (!file.exists(fil)) { download.file(url, fil) }
unzip(fil, exdir="ztca")

# read them in (this takes a bit)
ztca <- readShapePoly("ztca/cb_2014_us_zcta510_500k.shp", verbose=TRUE)

# extract NY
ny <- ztca[as.character(ztca$ZCTA5CE10) %in% as.character(zipcode[zipcode$state=="NY",]$zip),]

# your points
latlongdata <- 
  structure(list(dropoff_longitude = c(-73.981705, -73.993553, 
-73.973305, -73.988823, -73.938484, -74.015503, -73.95472, -73.9571, 
-73.971298, -73.99794), dropoff_latitude = c(40.760559, 40.756348, 
40.762646, 40.777504, 40.684692, 40.709881, 40.783371, 40.776657, 
40.752148, 40.720535)), row.names = c(8807218L, 9760613L, 3175671L, 
10878727L, 2025038L, 5345659L, 14474481L, 1650223L, 684883L, 
9129975L), class = "data.frame", .Names = c("dropoff_longitude", 
"dropoff_latitude"))

# make them all super spatial-like (must be in lon,lat format)
pts <- SpatialPoints(as.matrix(latlongdata[,1:2]))

# figure out where they are (this can take a bit)
dat <- pts %over% ny

# merge your data back in (there are many ways to do this)
dat$lon <- latlongdata$dropoff_longitude
dat$lat <- latlongdata$dropoff_latitude
rownames(dat) <- rownames(latlongdata)

# boom
dat
##          ZCTA5CE10     AFFGEOID10 GEOID10 ALAND10 AWATER10       lon      lat
## 8807218      10019 8600000US10019   10019 1780742        0 -73.98171 40.76056
## 9760613      10018 8600000US10018   10018  836253        0 -73.99355 40.75635
## 3175671      10022 8600000US10022   10022 1107169        0 -73.97330 40.76265
## 10878727     10069 8600000US10069   10069  249044        0 -73.98882 40.77750
## 2025038      11221 8600000US11221   11221 3582803        0 -73.93848 40.68469
## 5345659      10280 8600000US10280   10280  300652    38759 -74.01550 40.70988
## 14474481     10128 8600000US10128   10128 1206195        0 -73.95472 40.78337
## 1650223      10028 8600000US10028   10028  811363        0 -73.95710 40.77666
## 684883       10017 8600000US10017   10017  820953        0 -73.97130 40.75215
## 9129975      10013 8600000US10013   10013 1425085        0 -73.99794 40.72053

【讨论】：

这是完美的，谢谢！我知道纬度/经度对是特定于纽约地区的，有什么方法可以减少 ztca 的大小以加快速度？

【解决方案2】：

我使用这些方法将 lon-lat 转换为包含多边形：

library(maptools)
points.file<-readShapePoints("path.to.pts.shp")
poly.file<-read.ShapePoy("path.to.poly.shp")
points.file %over% poly.file

【讨论】：

谢谢。关于将纬度/经度对转换为形状文件的任何建议？
我确信在 R 中有一种方法可以做到这一点，但我不知道。我所做的是使用QGIS（开源 GIS 程序），它能够读取具有 lat-lon 规范的 .csv 文件，然后您可以将其保存为 shapefile。
另外，不确定readShapePoints 是否可以处理这种格式，但this 问题表明R 可以很容易地将您的经纬度转换为rgdal 文件，不管是什么...

【解决方案3】：

比较 dplyr rowwise 选项与 hrbrmstr 的 maptools 解决方案的速度，看起来 dplyr 胜出（至少在较小的数据集上）

# Test lat/long data
latlongdata <- 
structure(list(dropoff_longitude = c(-73.981705, -73.993553, 
-73.973305, -73.988823, -73.938484, -74.015503, -73.95472, -73.9571, 
-73.971298, -73.99794), dropoff_latitude = c(40.760559, 40.756348, 
40.762646, 40.777504, 40.684692, 40.709881, 40.783371, 40.776657, 
40.752148, 40.720535)), row.names = c(8807218L, 9760613L, 3175671L, 
10878727L, 2025038L, 5345659L, 14474481L, 1650223L, 684883L, 
9129975L), class = "data.frame", .Names = c("dropoff_longitude", 
"dropoff_latitude"))


# zipcode data with lat/lon coordinates
url <- "http://www.boutell.com/zipcodes/zipcode.zip"
fil <- "ziplatlong.zip"

# download an unzip
if (!file.exists(fil)) { download.file(url, fil) }
unzip(fil, exdir="zips")

library(readr)
ziplkp<-read_csv("zips/zipcode.csv")


# Method 1: dplyr + ZipLooker function

ZipLooker<-function(dropoff_longitude,dropoff_latitude){
  if(is.na(dropoff_longitude)|is.na(dropoff_latitude)){
    z<-NA_character_
  } else {

    tryCatch({
      x<-ziplkp
      x$latdiff=abs(dropoff_latitude-x$latitude)
      x$londiff=abs(dropoff_longitude-x$longitude)
      x$totdiff=x$latdiff+x$londiff
      z<-head(top_n(x,1,-totdiff),n=1)$zip
      return(z)
    }, error=function(e) NA_character_)
  }
}


latlongdata %>% 
  rowwise() %>% 
  mutate(zipcode=ZipLooker(dropoff_longitude,dropoff_latitude)
  )




# Method 2: maptools + sp
library(sp)
library(maptools)

# grab the zip code boundaries
url <- "http://www2.census.gov/geo/tiger/GENZ2014/shp/cb_2014_us_zcta510_500k.zip"
fil <- "ztca.zip"

# don't waste bandwidth
if (!file.exists(fil)) { download.file(url, fil) }
unzip(fil, exdir="ztca")

# read them in (this takes a bit)
ztca <- readShapePoly("ztca/cb_2014_us_zcta510_500k.shp", verbose=TRUE)


# extract NY
ny <- ztca[as.character(ztca$ZCTA5CE10) %in% as.character(ziplkp[ziplkp$state=="NY",]$zip),]

# make them all super spatial-like (must be in lon,lat format)
pts <- SpatialPoints(as.matrix(latlongdata[,1:2]))

# figure out where they are (this can take a bit)
dat <- pts %over% ny

# merge your data back in (there are many ways to do this)
dat$lon <- latlongdata$dropoff_longitude
dat$lat <- latlongdata$dropoff_latitude
rownames(dat) <- rownames(latlongdata)



# comparing the two (only the bulkiest parts)

library(microbenchmark)
microbenchmark(
dat <- pts %over% ny
,
latlongdata %>% 
  rowwise() %>% 
  mutate(zipcode=ZipLooker(dropoff_longitude,dropoff_latitude)
  )
,times = 10)

输出：

Unit: milliseconds

expr
dat <- pts %over% ny
latlongdata %>% rowwise() %>% mutate(zipcode = ZipLooker(dropoff_longitude,      dropoff_latitude))

       min        lq   median       uq      max neval
 275.89494 286.38187 297.9254 421.8727 445.7165    10
  95.18166  97.09873 101.8102 108.8677 122.0515    10

【讨论】：

我不应该说，经过几次抽查，maptools方法更准确。使用 dplyr 方法选择的一些邮政编码有一点偏差；有时它会选择一个相邻的拉链。