I have a dataset that lists all zip codes in the U.S., their types (standard, po box, university, etc). I want to replace po box and university zip codes with the next closest standard zip code. I broke down the dataset by state so that R wouldn't have to make as many calculations. In theory, I would like to have standard zip codes in the first column and zip codes that need replacement in the first row, and have the distance between the two be the intersection value.

For example,

                REP 1   REP 2   REP 3   REP 4
STD 1           0.215   0.152   0.025   0.124   
STD 2           0.365   0.410   0.074   0.234
STD 3           0.234   0.201   1.322   0.683
STD 4           0.543   0.282   0.483   0.094 
MINS            STD 1   STD 1   STD 2   STD 4

where STD 1 is a standard zip code with its own latitude and longitude, and REP 1 is a zip code that needs to be replaced (is a university/po box zip) with its own latitude and longitude. I only have about 5 weeks of experience in R, so please bear with me if something doesn't quite make sense to me immediately. I have tried to do this in excel and having a sheet with close to 10,000 columns by 40,000 rows crashes every time that I try to calculate all of the distances because there are just too many calculations.

I have a feeling that either the apply() or mapply() functions are needed here. I want to calculate the distances using a formula that considers the curvature of the earth, (euclidean, etc) like dist() or the geosphere package to maintain accuracy and be reproducible.

If there is anything else that would be helpful to add on here, let me know and I'll upload it asap. Here is my R code for Alaska, the first state in alphabetical order.


This sets up a table that has column names "ZIP_CODE", "Long", "Lat", and "IS_PO". "IS_PO" is a numerical indicator for whether or not the zip code is standard or po/university. 1 indicates that the zip code is a po/univ zip and 0 indicates a standard zip. I did this because some functions required that the data in the dataset be the same type (numerical).

Here are some of my failed attempts at writing code to calculate the minimum distances.

 lapply(bit::chunk(1, nrow(zipcode), 1e2), function(ridx) {
  merge(zipcode, zipcode[ridx[1]:ridx[2]], by = "dum", allow.cartesian = T)[
    , dist := distGeo(matrix(c(longitude.x, latitude.x), ncol = 2), 
                      matrix(c(longitude.y, latitude.y), ncol = 2))/1609.34 # meters to miles
    ][dist <= 5 # necessary distance treshold
      ][, dum := NULL]
}) %>% rbindlist -> zip_nearby_dt

DOESITWORK<-apply(db, 1, function(x) spDistsN1(matrix(x[3:4], nrow=1),

  • $\begingroup$ Zip codes are polygons/contiguous objects, so they don't have a natural distance from each other. If you want to find distances between zip codes you'll need to first identify a "center". This is obviously subjective to the use case. If you want to use curvature of the earth you'll need to project the data. $\endgroup$
    – Jon
    Jun 24, 2016 at 15:48
  • $\begingroup$ I was given latitudes and longitudes in my dataset that I will treat as the centers of the zip codes for calculating distances. $\endgroup$ Jun 24, 2016 at 18:56

1 Answer 1


I think I have read your question correctly it looks like you need a nearest neighbor implementation. If you are unfamiliar with the concept you can find the wiki article here https://en.wikipedia.org/wiki/Nearest_neighbor_search.

I went ahead a wrote an example implementation you can use as a guide. Please note that this is a brute force method and not useful for big data sets. Once you have a grasp of the material I suggest checking out some libraries like RANN that have "real" implementations.

Read in some random test data and clean For this test let us assume that we want to find the closest AMERICAN city for each location

coord_data = read.csv("~/Downloads/SalesJan2009.csv", stringsAsFactors = F)
coord_data$id = c(1:nrow(coord_data))
coord_data$is_usa = ifelse(coord_data$Country == "United States", 1, 0)
coord_data = coord_data[ , c("id", "Latitude", "Longitude", "is_usa")]
names(coord_data) = tolower(names(coord_data))

Define your distance function. Here we have geo-coordinates over long distance so Euclidean will not do. I am using the Law of Cosines to calculate great circle distance but Haversine and Vincenty should be considered given your needs. To learn more start here: https://en.wikipedia.org/wiki/Great-circle_distance.

greatCircleDistance = function(latAlpha, longAlpha, latBeta, longBeta, radius = 6371) {
  ## Function taken directly from Wikipedia
  ## Earth radius in km is default (6371)
  ## Long/Lats are in degrees so need helper function to convert to radians
  degreeToRadian = function(degree) (degree * pi / 180)
  deltaLong = degreeToRadian(longBeta) - degreeToRadian(longAlpha)
  sinLat = sin(degreeToRadian(latAlpha)) * sin(degreeToRadian(latBeta))
  cosLat = cos(degreeToRadian(latAlpha)) * cos(degreeToRadian(latBeta))
  ## acos is finicky with precision so we will assume if NA is thrown
  ## the argument was very close to 1 and therefore will return 0
  ## acos(1) == 0
  acosRaw = suppressWarnings(acos(sinLat + cosLat * cos(deltaLong)))
  acosSafe = ifelse(is.na(acosRaw), 0, acosRaw)
  acosSafe * radius

Distance between Basildon, UK and Parkville, US


 Returns [1] [1] 6929.351 km.

It matches Google's calc so we are good to go!

Brute Force Example: As you noticed in your Excel sheet this will blow up quickly as data set gets larger. There are much more efficient ways of implementing the search. One idea is to start with the the geo-data structure itself and write an R-Tree, but I'll leave that for you.

 bruteForceNearestNeighbor = function(geoData) {
      makeCoordinate = function(idx) {
        c("id" = idx, "latitude" = geoData$latitude[idx], "longitude" = geoData$longitude[idx])
      singleCoordMinDistance = function(coordinate, locations) {
        locationsUS = locations[locations$is_us == 1 & locations$id != coordinate["id"], ]
        distances = mapply(greatCircleDistance,
              latAlpha = coordinate["latitude"],
              longAlpha = coordinate["longitude"],
              latBeta = locationsUS$latitude,
              longBeta = locationsUS$longitude)
        closestIndex = which(distances == min(distances))
        locations[closestIndex, "id"]
      nearestNeighbors = vector("numeric", nrow(geoData))
      for ( i in 1:nrow(geoData) ) {
        coord = makeCoordinate(i)
        nearestNeighbors[i] = singleCoordMinDistance(coord, geoData)

    coord_data$nearest_neighbor = bruteForceNearestNeighbor(coord_data)
  • $\begingroup$ I believe that if I run this (or a similar function) on my dataset broken down by state that it would be a sane computational load. This is a great example and I can clearly see the support behind the idea. I'm still new to R, so I will have to learn how to write a user defined function, but this is a great start! $\endgroup$ Jun 27, 2016 at 14:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.