I have a set of locations in geographic coordinates, and I would like to group the points using hierarchical clustering followed by tree-cutting at various "heights" in order to calculate group-wise means of variables at recorded at each location.

Hierarchical clustering of the distance matrix of geographic coordinates, I presume, may be a misleading way to form groups because latitude and longitude are not equally spaced.

I can then imagine two ways forward:

  1. Using the great circle distance for the distance metric.
  2. Converting the geographic coordinates to an equally-scaled projection and then finding the Euclidean distance.

Apart from option two being more complicated to perform, are these approaches equivalent? And what exactly is the meaning of the tree cutting height in these cases?

  • 2
    Re 2: There's no such thing as an "equally-scaled projection": the distance matrix obtained with any projection will differ from that using geodesic distances. However, if all points are nearby--say, within a few thousand kilometers of each other--there will be many projections whose distance matrices are reasonably close approximations to the correct one. For working across the entire earth, a third option is available: consider using the 3D distance matrix. It is easy to compute and should, in general, give results similar to the geodesic distance. – whuber Apr 9 '12 at 0:36

Thanks to @whuber for setting me on the right track here. Looks as if there will be no additional answers forthcoming, so will settle this question by posting my own observations that may be useful for others learning about distances, clustering, and projections.

The following R code, using the geosphere, rgdal, and sp packages demonstrates that careful selection of the right projection can give an accurate distance matrix (where accurate is defined as geodesic distance) when points are up to 2000 km apart (axes are in metres).


## Produce 200 randomly positioned geographic coordinates
## in central Canada
xyLatLon <- data.frame(lon=(runif(200)*-30)-85,

## Convert to a Lambert Conformal Conic projection that should
## reasonably approximate the true distance
newProj <- "+proj=lcc +lat_1=49 +lat_2=77 +lat_0=63.390675
            +lon_0=-91.86666666666666 +x_0=6200000 +y_0=3000000
            +ellps=GRS80 +units=m +no_defs" 
xyLcc <- spTransform(SpatialPoints(xyLatLon, proj4string=CRS("+proj=longlat")), CRS(newProj))

## Find the geodesic distance matrix from geographic coordinates
## assuming the WGS84 ellipsoid
xyDist1 <- distm(xyLatLon, fun=distMeeus)

## Find the Euclidean distance matrix from the projection
xyDist2 <- as.matrix(dist(coordinates(xyLcc)))

## Find the Euclidean distance matrix of the geographic coordinates
xyDist3 <- as.matrix(dist(xyLatLon))

Plots of the elements of these three distance matrices are shown below. The plot on the left indicates that the projection selected is highly correlated with the geodesic distance across the range of distances used here. While the right plot demonstrates the considerable error that would be expected if unprojected geographic coordinates were to be used.

enter image description here

  • Hi @digitalmaps I'm trying to cluster temperature data in a regularly spaced lat-lon grid. I would like to take geographic distances into account, not only temp values. Do I need to calculate a distance matrix as you explain? Thanks – pacomet Nov 11 '16 at 12:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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