[Proj] proj4 string for perfect sphere

Fri Apr 13 19:16:41 PDT 2012

My thoughts:

1. You can compute distance between two lat/lon points on a sphere using the Haversine formula or Vincenty's formula specialized to the sphere.  They are described at:
   http://en.wikipedia.org/wiki/Great-circle_distance
Code for Haversine is below.  Fix it up as needed.

2. If all you need is an R-Tree, just use an R-Tree.  Don't bring in a database just because it has an R-Tree built inside, that will distort your entire programming effort needlessly.  Here is an R-Tree that works in 2 or 3 dimensions.  It worked for me in just a couple of hours of work:
   https://github.com/nushoin/RTree
NOTE: There are multiple versions of this file floating around the Internet.  Use what works for you.

3. Your real problem is that you want an R-Tree on the surface of a sphere, not on a plane.  Projecting your points to a plane and then using an R-Tree will only work if the points are within a localized area (and you have a good bound on how much your distances might get distored).  The reason it's not a global solution is that any projection has certain areas of the sphere where nearby points get projected to far-away points on the plane.  Because a sphere is periodic and a plane is not.

There are two ways to solve this:
 a) Develop an R-Tree algorithm on the sphere --- i.e. for a periodic type of geometry.  This would involve some work, and there's no obvious evidence that anyone has done it.  But I believe it would work.

 b) Use a 3-D R-Tree, based on the Cartesian coordinates of your points in 3-D.  That approach is described here, and apparently works:
   http://blog.cleverelephant.ca/2009/11/postgis-gets-spherical-directors-cut.html

I quote: "How does this magic index work? The "trick" is to be as lazy as possible. I didn't want to write a whole new indexing scheme, and I already had a serviceable R-Tree index. What I needed to do was convert the lat/lon coordinates to a domain where the problems of the singularities at the poles and dateline would go away. The answer is to convert from spherical coordinates (lat/lon) relative to Greenwich into geocentric coordinates (x/y/z) relative to the center of earth. It's easy then to build a 3D R-Tree on the geocentric bounds of the features. Calculating the bounds in 3D is tricky, because the curvature of the features over the spherical surface must be respected, but once that is done, the index performs admirably."

Unfortunately, the author (Paul Ramsey) does not explain how he calculates his bounds in 3D.  I would pester him and find out.  Would be interested in hearing back if you do.

Either way, projections aren't needed or involved.  Hence, you shouldn't be using proj.4.  Using proj.4 on a "no-op" spherical projection doesn't make mush sense.  If all you want is a distance formula, that's really easy without proj.4.  Nor does a database really buy you much either, as far as I can tell.

Or, you could use the GEOGRAPHY data type that Paul Ramsey wrote about implementing in PostgreSQL:
    http://blog.opengeo.org/2009/11/04/postgis-gets-spherical/

If you end up using PostgreSQL in this way, I would not look at is a black box.  I would instead look at it as a convenient package of algorithms, already coded up, that you basically understand.

-- Bob

============================================
D2R = M_PI / 180
inline double sqr(double x) { return x*x; }

/** See: http://www.cs.nyu.edu/visual/home/proj/tiger/gisfaq.html
@return Distance (in degrees) */
extern double haversine_distance(
double lon1_deg, double lat1_deg,
double lon2_deg, double lat2_deg)
{
	// Convert inputs to degrees
	double lon1 = lon1_deg * D2R;
	double lat1 = lat1_deg * D2R;
	double lon2 = lon2_deg * D2R;
	double lat2 = lat2_deg * D2R;

	// Apply the Haversine Formula
	double dlon = lon2 - lon1;
	double dlat = lat2 - lat1;
	double a = sqr(sin(dlat/2)) + cos(lat1) * cos(lat2) * sqr(sin(dlon/2));
	double c = 2 * atan2( sqrt(a), sqrt(1-a) );
	return c * R2D;		// Convert to degrees
}

________________________________________
From: proj-bounces at lists.maptools.org [proj-bounces at lists.maptools.org] On Behalf Of Demitri Muna [thatsanicehatyouhave at me.com]
Sent: Friday, April 13, 2012 9:24 PM
To: PROJ.4 and general Projections Discussions
Subject: Re: [Proj] proj4 string for perfect sphere

Hi,

Quick update. Once I recognized 6371000 as the radius of the Earth in meters, I created my own SRID with these parameters (r = 1 radian as degrees):

"+proj=longlat +a=57.2957795 +b=57.2957795 +no_defs"

This worked perfectly:

> SELECT GreatCircleLength(GeomFromText('LINESTRING(2.47254166667 25.9237777778, 2.77729166667 -12.1073136111)',50000));
38.0322481765926

The correct answer! I would still appreciate your thoughts/advice on what I'm trying to do more generally.

Cheers,
Demitri
_______________________________________________
Proj mailing list
Proj at lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj