postgis/postgis
postgis/postgis
2
6
Fork 4
Code Issues Pull Requests Releases Activity
postgis/postgis/lwgeom_spheroid.c
Paul Ramsey a80b805b09 Move to snprintf where easily done
2022-07-25 14:58:45 -07:00

570 lines
15 KiB
C
Raw Permalink Blame History

/**********************************************************************
*
* PostGIS - Spatial Types for PostgreSQL
* http://postgis.net
*
* PostGIS is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* PostGIS is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
*
**********************************************************************
*
* Copyright (C) 2001-2003 Refractions Research Inc.
*
**********************************************************************/
#include "postgres.h"
#include <math.h>
#include <float.h>
#include <string.h>
#include <stdio.h>
#include "access/gist.h"
#include "access/itup.h"
#include "fmgr.h"
#include "utils/elog.h"
#include "../postgis_config.h"
#include "liblwgeom.h"
#include "lwgeom_pg.h"
#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);
Datum ellipsoid_out(PG_FUNCTION_ARGS);
Datum LWGEOM_length2d_ellipsoid(PG_FUNCTION_ARGS);
Datum LWGEOM_length_ellipsoid_linestring(PG_FUNCTION_ARGS);
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
* ie. SPHEROID["name",A,rf] or SPHEROID("name",A,rf)
* SPHEROID["GRS_1980",6378137,298.257222101]
* wkt says you can use "(" or "["
*/
PG_FUNCTION_INFO_V1(ellipsoid_in);
Datum ellipsoid_in(PG_FUNCTION_ARGS)
{
char *str = PG_GETARG_CSTRING(0);
SPHEROID *sphere = (SPHEROID *) palloc(sizeof(SPHEROID));
int nitems;
double rf;
memset(sphere,0, sizeof(SPHEROID));
if (strstr(str,"SPHEROID") != str )
{
elog(ERROR,"SPHEROID parser - doesn't start with SPHEROID");
pfree(sphere);
PG_RETURN_NULL();
}
nitems = sscanf(str,"SPHEROID[\"%19[^\"]\",%lf,%lf]",
sphere->name, &sphere->a, &rf);
if ( nitems==0)
nitems = sscanf(str,"SPHEROID(\"%19[^\"]\",%lf,%lf)",
sphere->name, &sphere->a, &rf);
if (nitems != 3)
{
elog(ERROR,"SPHEROID parser - couldnt parse the spheroid");
pfree(sphere);
PG_RETURN_NULL();
}
sphere->f = 1.0/rf;
sphere->b = sphere->a - (1.0/rf)*sphere->a;
sphere->e_sq = ((sphere->a*sphere->a) - (sphere->b*sphere->b)) /
(sphere->a*sphere->a);
sphere->e = sqrt(sphere->e_sq);
PG_RETURN_POINTER(sphere);
}
PG_FUNCTION_INFO_V1(ellipsoid_out);
Datum ellipsoid_out(PG_FUNCTION_ARGS)
{
SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(0);
char *result;
size_t sz = MAX_DIGS_DOUBLE + MAX_DIGS_DOUBLE + 20 + 9 + 2;
result = palloc(sz);
snprintf(result, sz, "SPHEROID(\"%s\",%.15g,%.15g)",
sphere->name, sphere->a, 1.0/sphere->f);
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));
}
/*
* Find the "length of a geometry"
* length2d_spheroid(point, sphere) = 0
* length2d_spheroid(line, sphere) = length of line
* length2d_spheroid(polygon, sphere) = 0
* -- could make sense to return sum(ring perimeter)
* uses ellipsoidal math to find the distance
* x's are longitude, and y's are latitude - both in decimal degrees
*/
PG_FUNCTION_INFO_V1(LWGEOM_length2d_ellipsoid);
Datum LWGEOM_length2d_ellipsoid(PG_FUNCTION_ARGS)
{
GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0);
SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(1);
LWGEOM *lwgeom = lwgeom_from_gserialized(geom);
double dist = lwgeom_length_spheroid(lwgeom, sphere);
lwgeom_free(lwgeom);
PG_FREE_IF_COPY(geom, 0);
PG_RETURN_FLOAT8(dist);
}
/*
* Find the "length of a geometry"
*
* length2d_spheroid(point, sphere) = 0
* length2d_spheroid(line, sphere) = length of line
* length2d_spheroid(polygon, sphere) = 0
* -- could make sense to return sum(ring perimeter)
* uses ellipsoidal math to find the distance
* x's are longitude, and y's are latitude - both in decimal degrees
*/
PG_FUNCTION_INFO_V1(LWGEOM_length_ellipsoid_linestring);
Datum LWGEOM_length_ellipsoid_linestring(PG_FUNCTION_ARGS)
{
GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0);
LWGEOM *lwgeom = lwgeom_from_gserialized(geom);
SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(1);
double length = 0.0;
/* EMPTY things have no length */
if ( lwgeom_is_empty(lwgeom) )
{
lwgeom_free(lwgeom);
PG_RETURN_FLOAT8(0.0);
}
length = lwgeom_length_spheroid(lwgeom, sphere);
lwgeom_free(lwgeom);
PG_FREE_IF_COPY(geom, 0);
/* Something went wrong... */
if ( length < 0.0 )
{
elog(ERROR, "lwgeom_length_spheroid returned length < 0.0");
PG_RETURN_NULL();
}
/* Clean up */
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 comming 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)
*/
PG_FUNCTION_INFO_V1(geometry_distance_spheroid);
Datum geometry_distance_spheroid(PG_FUNCTION_ARGS)
{
GSERIALIZED *geom1 = PG_GETARG_GSERIALIZED_P(0);
GSERIALIZED *geom2 = PG_GETARG_GSERIALIZED_P(1);
SPHEROID *sphere = (SPHEROID *)PG_GETARG_POINTER(2);
int type1 = gserialized_get_type(geom1);
int type2 = gserialized_get_type(geom2);
bool use_spheroid = PG_GETARG_BOOL(3);
LWGEOM *lwgeom1, *lwgeom2;
double distance;
gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
/* Calculate some other parameters on the spheroid */
spheroid_init(sphere, sphere->a, sphere->b);
/* Catch sphere special case and re-jig spheroid appropriately */
if ( ! use_spheroid )
{
sphere->a = sphere->b = sphere->radius;
}
if ( ! ( type1 == POINTTYPE || type1 == LINETYPE || type1 == POLYGONTYPE ||
type1 == MULTIPOINTTYPE || type1 == MULTILINETYPE || type1 == MULTIPOLYGONTYPE ))
{
elog(ERROR, "geometry_distance_spheroid: Only point/line/polygon supported.\n");
PG_RETURN_NULL();
}
if ( ! ( type2 == POINTTYPE || type2 == LINETYPE || type2 == POLYGONTYPE ||
type2 == MULTIPOINTTYPE || type2 == MULTILINETYPE || type2 == MULTIPOLYGONTYPE ))
{
elog(ERROR, "geometry_distance_spheroid: Only point/line/polygon supported.\n");
PG_RETURN_NULL();
}
/* Get #LWGEOM structures */
lwgeom1 = lwgeom_from_gserialized(geom1);
lwgeom2 = lwgeom_from_gserialized(geom2);
/* We are going to be calculating geodetic distances */
lwgeom_set_geodetic(lwgeom1, LW_TRUE);
lwgeom_set_geodetic(lwgeom2, LW_TRUE);
distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, sphere, 0.0);
PG_RETURN_FLOAT8(distance);
}
PG_FUNCTION_INFO_V1(LWGEOM_distance_ellipsoid);
Datum LWGEOM_distance_ellipsoid(PG_FUNCTION_ARGS)
{
SPHEROID s;
/* No spheroid provided */
if (PG_NARGS() == 2) {
/* Init to WGS84 */
spheroid_init(&s, 6378137.0, 6356752.314245179498);
PG_RETURN_DATUM(DirectFunctionCall4(geometry_distance_spheroid,
PG_GETARG_DATUM(0),
PG_GETARG_DATUM(1),
PointerGetDatum(&s),
BoolGetDatum(true)));
}
PG_RETURN_DATUM(DirectFunctionCall4(geometry_distance_spheroid,
PG_GETARG_DATUM(0),
PG_GETARG_DATUM(1),
PG_GETARG_DATUM(2),
BoolGetDatum(true)));
}
PG_FUNCTION_INFO_V1(LWGEOM_distance_sphere);
Datum LWGEOM_distance_sphere(PG_FUNCTION_ARGS)
{
SPHEROID s;
/* Init to WGS84 */
spheroid_init(&s, 6378137.0, 6356752.314245179498);
if (PG_NARGS() == 3) {
s.radius = PG_GETARG_FLOAT8(2);
}
s.a = s.b = s.radius;
PG_RETURN_DATUM(DirectFunctionCall4(geometry_distance_spheroid,
PG_GETARG_DATUM(0), PG_GETARG_DATUM(1), PointerGetDatum(&s), BoolGetDatum(false)));
}
View Git Blame Copy Permalink