[SCM] PostGIS branch master updated. 3.6.0rc2-328-g760c7eb39
git at osgeo.org
git at osgeo.org
Thu Feb 5 16:03:59 PST 2026
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "PostGIS".
The branch, master has been updated
via 760c7eb3951d48547c947f00cac81bc1e2ad397a (commit)
from bf36aef0ca85b254b8576f0936d377612a4668c0 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 760c7eb3951d48547c947f00cac81bc1e2ad397a
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Thu Feb 5 16:03:49 2026 -0800
Remove dead spherical code
diff --git a/postgis/lwgeom_spheroid.c b/postgis/lwgeom_spheroid.c
index 8c11edff5..572470fe2 100644
--- a/postgis/lwgeom_spheroid.c
+++ b/postgis/lwgeom_spheroid.c
@@ -44,12 +44,6 @@
#define SHOW_DIGS_DOUBLE 15
#define MAX_DIGS_DOUBLE (SHOW_DIGS_DOUBLE + 6 + 1 + 3 +1)
-/*
- * distance from -126 49 to -126 49.011096139863
- * in 'SPHEROID["GRS_1980",6378137,298.257222101]'
- * is 1234.000
- */
-
/* PG-exposed */
Datum ellipsoid_in(PG_FUNCTION_ARGS);
@@ -60,15 +54,6 @@ Datum LWGEOM_distance_ellipsoid(PG_FUNCTION_ARGS);
Datum LWGEOM_distance_sphere(PG_FUNCTION_ARGS);
Datum geometry_distance_spheroid(PG_FUNCTION_ARGS);
-/* internal */
-double distance_sphere_method(double lat1, double long1,double lat2,double long2, SPHEROID *sphere);
-double distance_ellipse_calculation(double lat1, double long1, double lat2, double long2, SPHEROID *sphere);
-double distance_ellipse(double lat1, double long1, double lat2, double long2, SPHEROID *sphere);
-double deltaLongitude(double azimuth, double sigma, double tsm,SPHEROID *sphere);
-double mu2(double azimuth,SPHEROID *sphere);
-double bigA(double u2);
-double bigB(double u2);
-
/*
* Use the WKT definition of an ellipsoid
@@ -131,196 +116,6 @@ Datum ellipsoid_out(PG_FUNCTION_ARGS)
PG_RETURN_CSTRING(result);
}
-/*
- * support function for distance calc
- * code is taken from David Skea
- * Geographic Data BC, Province of British Columbia, Canada.
- * Thanks to GDBC and David Skea for allowing this to be
- * put in PostGIS.
- */
-double
-deltaLongitude(double azimuth, double sigma, double tsm,SPHEROID *sphere)
-{
- /* compute the expansion C */
- double das,C;
- double ctsm,DL;
-
- das = cos(azimuth)*cos(azimuth);
- C = sphere->f/16.0 * das * (4.0 + sphere->f * (4.0 - 3.0 * das));
-
- /* compute the difference in longitude */
- ctsm = cos(tsm);
- DL = ctsm + C * cos(sigma) * (-1.0 + 2.0 * ctsm*ctsm);
- DL = sigma + C * sin(sigma) * DL;
- return (1.0 - C) * sphere->f * sin(azimuth) * DL;
-}
-
-
-/*
- * support function for distance calc
- * code is taken from David Skea
- * Geographic Data BC, Province of British Columbia, Canada.
- * Thanks to GDBC and David Skea for allowing this to be
- * put in PostGIS.
- */
-double
-mu2(double azimuth,SPHEROID *sphere)
-{
- double e2;
-
- e2 = sqrt(sphere->a*sphere->a-sphere->b*sphere->b)/sphere->b;
- return cos(azimuth)*cos(azimuth) * e2*e2;
-}
-
-
-/*
- * Support function for distance calc
- * code is taken from David Skea
- * Geographic Data BC, Province of British Columbia, Canada.
- * Thanks to GDBC and David Skea for allowing this to be
- * put in PostGIS.
- */
-double
-bigA(double u2)
-{
- return 1.0 + u2/256.0 * (64.0 + u2 * (-12.0 + 5.0 * u2));
-}
-
-
-/*
- * Support function for distance calc
- * code is taken from David Skea
- * Geographic Data BC, Province of British Columbia, Canada.
- * Thanks to GDBC and David Skea for allowing this to be
- * put in PostGIS.
- */
-double
-bigB(double u2)
-{
- return u2/512.0 * (128.0 + u2 * (-64.0 + 37.0 * u2));
-}
-
-
-
-double
-distance_ellipse(double lat1, double long1,
- double lat2, double long2, SPHEROID *sphere)
-{
- double result = 0;
-#if POSTGIS_DEBUG_LEVEL >= 4
- double result2 = 0;
-#endif
-
- if ( (lat1==lat2) && (long1 == long2) )
- {
- return 0.0; /* same point, therefore zero distance */
- }
-
- result = distance_ellipse_calculation(lat1,long1,lat2,long2,sphere);
-
-#if POSTGIS_DEBUG_LEVEL >= 4
- result2 = distance_sphere_method(lat1, long1,lat2,long2, sphere);
-
- POSTGIS_DEBUGF(4, "delta = %lf, skae says: %.15lf,2 circle says: %.15lf",
- (result2-result),result,result2);
- POSTGIS_DEBUGF(4, "2 circle says: %.15lf",result2);
-#endif
-
- if (result != result) /* NaN check
- * (x==x for all x except NaN by IEEE definition)
- */
- {
- result = distance_sphere_method(lat1, long1,
- lat2,long2, sphere);
- }
-
- return result;
-}
-
-/*
- * Given 2 lat/longs and ellipse, find the distance
- * note original r = 1st, s=2nd location
- */
-double
-distance_ellipse_calculation(double lat1, double long1,
- double lat2, double long2, SPHEROID *sphere)
-{
- /*
- * Code is taken from David Skea
- * Geographic Data BC, Province of British Columbia, Canada.
- * Thanks to GDBC and David Skea for allowing this to be
- * put in PostGIS.
- */
-
- double L1,L2,sinU1,sinU2,cosU1,cosU2;
- double dl,dl1,dl2,dl3,cosdl1,sindl1;
- double cosSigma,sigma,azimuthEQ,tsm;
- double u2,A,B;
- double dsigma;
-
- double TEMP;
-
- int iterations;
-
-
- L1 = atan((1.0 - sphere->f ) * tan( lat1) );
- L2 = atan((1.0 - sphere->f ) * tan( lat2) );
- sinU1 = sin(L1);
- sinU2 = sin(L2);
- cosU1 = cos(L1);
- cosU2 = cos(L2);
-
- dl = long2- long1;
- dl1 = dl;
- cosdl1 = cos(dl);
- sindl1 = sin(dl);
- iterations = 0;
- do
- {
- cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosdl1;
- sigma = acos(cosSigma);
- azimuthEQ = asin((cosU1 * cosU2 * sindl1)/sin(sigma));
-
- /*
- * Patch from Patrica Tozer to handle minor
- * mathematical stability problem
- */
- TEMP = cosSigma - (2.0 * sinU1 * sinU2)/
- (cos(azimuthEQ)*cos(azimuthEQ));
- if (TEMP > 1)
- {
- TEMP = 1;
- }
- else if (TEMP < -1)
- {
- TEMP = -1;
- }
- tsm = acos(TEMP);
-
-
- /* (old code?)
- tsm = acos(cosSigma - (2.0 * sinU1 * sinU2)/(cos(azimuthEQ)*cos(azimuthEQ)));
- */
-
- dl2 = deltaLongitude(azimuthEQ, sigma, tsm,sphere);
- dl3 = dl1 - (dl + dl2);
- dl1 = dl + dl2;
- cosdl1 = cos(dl1);
- sindl1 = sin(dl1);
- iterations++;
- }
- while ( (iterations<999) && (fabs(dl3) > 1.0e-032));
-
- /* compute expansions A and B */
- u2 = mu2(azimuthEQ,sphere);
- A = bigA(u2);
- B = bigB(u2);
-
- /* compute length of geodesic */
- dsigma = B * sin(sigma) * (cos(tsm) +
- (B*cosSigma*(-1.0 + 2.0 * (cos(tsm)*cos(tsm))))/4.0);
- return sphere->b * (A * (sigma - dsigma));
-}
/*
@@ -385,95 +180,6 @@ Datum LWGEOM_length_ellipsoid_linestring(PG_FUNCTION_ARGS)
PG_RETURN_FLOAT8(length);
}
-/*
- * For some lat/long points, the above method doesn't calculate the distance very well.
- * Typically this is for two lat/long points that are very very close together (<10cm).
- * This gets worse closer to the equator.
- *
- * This method works very well for very close together points, not so well if they're
- * far away (>1km).
- *
- * METHOD:
- * We create two circles (with Radius R and Radius S) and use these to calculate the distance.
- *
- * The first (R) is basically a (north-south) line of longitude.
- * Its radius is approximated by looking at the ellipse. Near the equator R = 'a' (earth's major axis)
- * near the pole R = 'b' (earth's minor axis).
- *
- * The second (S) is basically a (east-west) line of latitude.
- * Its radius runs from 'a' (major axis) at the equator, and near 0 at the poles.
- *
- *
- * North pole
- * *
- * *
- * *\--S--
- * * R +
- * * \ +
- * * A\ +
- * * ------ \ Equator/centre of earth
- * *
- * *
- * *
- * *
- * *
- * *
- * South pole
- * (side view of earth)
- *
- * Angle A is lat1
- * R is the distance from the centre of the earth to the lat1/long1 point on the surface
- * of the Earth.
- * S is the circle-of-latitude. Its calculated from the right triangle defined by
- * the angle (90-A), and the hypotenuse R.
- *
- *
- *
- * Once R and S have been calculated, the actual distance between the two points can be
- * calculated.
- *
- * We dissolve the vector from lat1,long1 to lat2,long2 into its X and Y components (called DeltaX,DeltaY).
- * The actual distance that these angle-based measurements represent is taken from the two
- * circles we just calculated; R (for deltaY) and S (for deltaX).
- *
- * (if deltaX is 1 degrees, then that distance represents 1/360 of a circle of radius S.)
- *
- *
- * Parts taken from PROJ - geodetic_to_geocentric() (for calculating Rn)
- *
- * remember that lat1/long1/lat2/long2 are coming in a *RADIANS* not degrees.
- *
- * By Patricia Tozer and Dave Blasby
- *
- * This is also called the "curvature method".
- */
-
-double distance_sphere_method(double lat1, double long1,double lat2,double long2, SPHEROID *sphere)
-{
- double R,S,X,Y,deltaX,deltaY;
-
- double distance = 0.0;
- double sin_lat = sin(lat1);
- double sin2_lat = sin_lat * sin_lat;
-
- double Geocent_a = sphere->a;
- double Geocent_e2 = sphere->e_sq;
-
- R = Geocent_a / (sqrt(1.0e0 - Geocent_e2 * sin2_lat));
- /* 90 - lat1, but in radians */
- S = R * sin(M_PI_2 - lat1) ;
-
- deltaX = long2 - long1; /* in rads */
- deltaY = lat2 - lat1; /* in rads */
-
- /* think: a % of 2*pi*S */
- X = deltaX/(2.0*M_PI) * 2 * M_PI * S;
- Y = deltaY/(2.0*M_PI) * 2 * M_PI * R;
-
- distance = sqrt((X * X + Y * Y));
-
- return distance;
-}
/*
* distance (geometry,geometry, sphere)
-----------------------------------------------------------------------
Summary of changes:
postgis/lwgeom_spheroid.c | 294 ----------------------------------------------
1 file changed, 294 deletions(-)
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list