How do I find the minimum distance between zip codes in R?

Question

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.

AK<-subset(db,STAABBRV.x=="AK")
AKPO<-subset(AK,ZipCodeType!="STANDARD",select=c("ZIP_CODE","ZipCodeType","Long","Lat"))
AKPO<-within(AKPO,{IS_PO=ifelse(ZipCodeType!="STANDARD",1,0)})
AKSTANDARD<-subset(AK,ZipCodeType=="STANDARD",select=c("ZIP_CODE","ZipCodeType","Long","Lat"))
AKSTANDARD<-within(AKSTANDARD,{IS_PO=ifelse(ZipCodeType!="STANDARD",1,0)})
table<-rbind(AKSTANDARD,AKPO)
table$ZipCodeType<-NULL
rm(AK,AKPO,AKSTANDARD)

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),
                                    x[5:6],
                                    longlat=TRUE)) 



mins<-apply(Lat,1,function(x)return(array(which.min(x))))
mins<-data.frame(row=names(mins),col=mins) 
Lat$mins<-apply(mins,1,FUN=function(x)return(paste(x["row"],colnames(Lat[as.numeric(x["col"])]),Lat[x["row"],as.numeric(x["col"])],sep="/")))

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

greatCircleDistance(coord_data$latitude[1],
coord_data$longitude[1],
coord_data$latitude[2],
coord_data$longitude[2])

 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)
      }
      nearestNeighbors
    }

    coord_data$nearest_neighbor = bruteForceNearestNeighbor(coord_data)