Investigating PostGIS functions to see what they are doing?

Question

Let's say you have a function like ST_Distance and I want to know how it works. The docs don't say much at all, how do I go about debugging this or figuring out how it ticks.

Answer 1

Not sure if this is helpful, let's go down the rabbit hole. First we need to know what's happening when you call the function in SQL. To do this we reference the output of \dS+. \dS+ shows the override table with one entry for every function and prototype, as well as the function that it is dispatching to.

1. You're calling ST_Distance($1::geog,$2::geog).
2. ST_Distance($1::geog,$2::geog) is calling _st_distance($1::geog, $2::geog, 0.0, true)
3. _st_distance($1::geog, $2::geog, 0.0, true) is calling a C function called geography_distance

From there we need to start looking into the C functions,

4. geography_distance is a self-named function on :201 in postgis/postgis/geography_measurement.c
5. geography_distance then calls lwgeom_distance_spheroid(lwgeom1, lwgeom2, &s, tolerance);
6. lwgeom_distance_spheroid is defined in the liblwgeom used by PostGIS on :2091 of lwgeodetic.c.
7. Because you're giving it two points, that calls ptarray_distance_spheroid which is defined here
8. pt_array_distance_spheroid internally calls liblwgeom's spheroid_distance and finally..

we have code

/**
* Computes the shortest distance along the surface of the spheroid
* between two points, using the inverse geodesic problem from
* GeographicLib (Karney 2013).
*
* @param a - location of first point
* @param b - location of second point
* @param s - spheroid to calculate on
* @return spheroidal distance between a and b in spheroid units
*/
double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
{
    struct geod_geodesic gd;
    geod_init(&gd, spheroid->a, spheroid->f);
    double lat1 = a->lat * 180.0 / M_PI;
    double lon1 = a->lon * 180.0 / M_PI;
    double lat2 = b->lat * 180.0 / M_PI;
    double lon2 = b->lon * 180.0 / M_PI;
    double s12; /* return distance */
    geod_inverse(&gd, lat1, lon1, lat2, lon2, &s12, 0, 0);
    return s12;
}

So that's what eventually get's called and does the lifting when you call ST_Distance. The beef of that is geod_inverse. Not seeing that in lwspheroid.c, I googled and it seems to come from GeographicLib. Looking back at lwspheroid.c you can see that being included here,

/* GeographicLib */
#if PROJ_GEODESIC
#include <geodesic.h>
#endif

And GeographicLib provides an implementation in C of the geodesic algorithms described in C. F. F. Karney, Algorithms for geodesics, J. Geodesy 87, 43–55 (2013); DOI: 10.1007/s00190-012-0578-z you can see geod_inverse documented here. You can find the implementation here on line :1054

This also matches a patch accepted only 17 months ago, here to move to use to GeographicLib.

So that's what it looks like to debug PostGIS and Postgres from the top to the bottom.

Answer 2

You can also download and debug the source code from Git repository. I personally looked up the function definitions by debugging the source code in Visual Studio and viewing the collaboration diagrams available at http://postgis.net/docs/doxygen/2.1/annotated.html