Projecting latitude and longitude such that the Euclidean distance is a good approximation of geographic distance

Question

Goal:

I want to run various algorithms in Shapely on some geographic data that I have. It is acceptable to have some measure of error, as long as it is "reasonable". e.g. find closest pair.

Method:

I would like to project the points from latitude-longitude to cartesian space such that

euclidan_distance(P(p0), P(p1)) ~= geodesic_distance(p0, p1)

I was thinking of using PyProj but it may be overkill, and it doesn't seem the easiest thing to do.

Naive approach:

y = latitude * 110574
x = longitude * 11320 * cos(radians(latitude))

(numbers from https://en.wikipedia.org/wiki/Latitude#Length_of_a_degree_of_latitude )

Results of naive approach:

I put in a set of coordinates nearby, and got relative errors of 0.5%-15.3%. 15% seems excessive.

Answer 1

You could use one of the "AUTO" projections described by Appendix E of the WMS specification

1. AUTO:42001 - Universal Transverse Mercator
2. AUTO:42002 - Transverse Mercator
3. AUTO:42003 - Orthographic
4. AUTO:42004 - Equirectangular

The specification gives WKT projection strings for them all (all you have to do is add the coordinates of your centre point). They are supported in MapServer so I would expect them to be available in all Proj4+ based systems.