[GIS] How to compute distance between 4500 points fast

arcgis-10.0arcgis-desktopgoogle earthkmzspherical-geometry

I am working on a location problem with 150 customers (points) and 30 facilities (points). I have determined the location of any customer and facility by google earth and then saved them in a KMZ file. I need to compute the distance between any pair of customers and facilities. So I should compute 150*30=4500 values by ruler tool in Google Earth. Is there any other way to obtain these values faster and easier? Even by converting the KMZ file to a layer in ArcGIS 10?

Best Answer

One of the fastest and easiest possible solutions uses a short program written with the free open source program R (the R project for statistical computing). The following code computes the distance matrix (using spherical distances) between two arrays of (lon, lat) coordinates named customers and facilities and stores it in an array distances (with rows for customers and columns for facilities).

dist <- function(x, y=c(0,0), radius=6366710) {
  # `x` and `y` are both (lon, lat) in radians.
  # `x` is an `n` by 2 matrix.
  # Returns distances on a sphere of radius `radius`; the default value
  # is a spherical model in meters.
  return (2 * radius * asin( sqrt(sin((x[,2]-y[2])/2)^2 + 
                          cos(x[,2])*cos(y[2])*(sin(x[,1]-y[1])/2)^2)))
}
distances <- apply(facilities * pi/180, 1, function(y) dist(customers * pi/180, y)

It runs reasonably quickly: your problem with 150 customers and 30 facilities generates a matrix of all 150*30 = 4500 distances in 2.5 milliseconds.

To extract the coordinates from a Google Earth file, save the locations directly in KML format or unzip the KMZ file (which will produce a KML version). This is an ASCII file containing (lon, lat) coordinates.. The maptools package reads such files directly. Here is a working example that reads KML files of customers and facilities and computes all distances between them.

read.kml <- function(s) {
  require(maptools)
  x <- getKMLcoordinates(s, ignoreAltitude=TRUE)
  return (matrix(unlist(x), ncol=2))
}
customers <- read.kml("F:/temp/customers.kml")
facilities <- read.kml("F:/temp/facilities.kml")
distances <- apply(facilities * pi/180, 1, function(y) dist(customers * pi/180, y))

At this point you can perform just about any kind of calculation you might like with the distances and--of course--you can write them to a file for post-processing on another platform.

Related Question