0

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)
[1] "+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?

Improve this question

1 Answer 1

Reset to default
1

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

Improve this answer
3
  • 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, 2018 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, 2018 at 21:48
  • Yeah I don't actually want to transform, just wanted to understand more. Thanks, that clarifies things
    – user119779
    May 1, 2018 at 23:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.