# How to calculate great circle distance between unprojected SpatialPolygons?

I'm working (in R) with a SpatialPolygonsDataFrame describing US historical counties, and I would like to calculate great circle distances between centroids of counties (polygons).

I've looked into functions like `spDists` from the `sp` package, `gDistance` from `rgeos`, and `distGeo` from `geosphere`, but I'm confused about what exactly they're calculating. For instance, the documentation for `spDists` states that it calculates Euclidean distance if the spatial object is not projected, and great circle distance if the spatial object is projected. Similarly, `gDistance` rejects my data because it "expects planar coordinates". This is my first source of confusion because I would think it should be the opposite, eg isn't it more natural to calculate Euclidean distance from projected data because it is embedded in two-dimensional space? Similarly, shouldn't it be easier to calculate great circle distance from longitudes and latitudes (my data) than from coordinates for some projection that may have distorted distance?

I don't mind projecting my data to use these functions, but I'm not sure exactly what comprises a "projection." The output of `proj4string` for my data is

``````> proj4string(US)
 "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
``````

I've seen in online tutorials that there's usually a "+init=..." component. Is that why it's not projected? Does it matter what "init" value I use to project it?

Here's what ?spDists says:

`````` longlat: logical; if FALSE, Euclidean distance, if TRUE Great Circle
(WGS84 ellipsoid) distance; if ‘x’ is a Spatial object,
longlat should not be specified but will be derived from
is.projected‘(x)’
``````

So if the data is not in latitude-longitude then it uses Euclidean distance in 2-dimensions. If the data is in lat-long then it uses great circle distances.

This is not unreasonable. If I am working in the UK, I will probably be using our local grid system because its a smallish area, and the approximation of the sphere to the plane is good enough. I can compute distances using Pythagoras and not get too much error with distance.

If I'm using lat-long coordinates I might be working on the whole globe, and so distances computed by Pythagoras' theorem on lat-long coordinates would be grossly misleading. I need to use great circle distances.

(The `+init` parameter is mostly used as a shortcut for specifying projections and coordinate systems. eg. "+init=epsg:4326" is equivalent to "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0". The projection string may or may not include the `+init` along with the other parameters - this is usually a bit redundant but tolerated...)

• Thanks, that makes sense. I had assumed that `longlat` was taking the literal value of `is.projected(x)`, but looking at the code, I can see that it's using the negation. So is it correct to say that my projection string (or the equivalent "+init=epsg:4326") specifies a coordinate system but not a projection? If I wanted to project it, what would I need to add to my projection string? – user119779 May 1 '18 at 21:29
• Yes, that's an unprojected geographic coordinate system definition. You can't just add to a projection string, you have to transform to a completely new projection string, using `spTransform`, a process which changes the coordinates stored in the object. Which projection to use depends on where your data is and what distortion you can cope with. But your question did say you wanted to do great circle distances, so why transform? – Spacedman May 1 '18 at 21:48
• Yeah I don't actually want to transform, just wanted to understand more. Thanks, that clarifies things – user119779 May 1 '18 at 23:56