Appropriate distance metric for spatial clustering of geographic coordinates

Question

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?

Answer 1

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).

library(sp)
library(rgdal)
library(geosphere)

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

## 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.