3

I`m trying to implement postgis St_Project analog on C++ (calculation method is taken from here - http://www.movable-type.co.uk/scripts/latlong.html):

#include <cmath>

struct GeoPoint
{
    double longitude;
    double latitude;
};

static const double EARTH_RADIUS = 6378137.0;
static double degToRad(double deg) { return (deg * M_PI) / 180.0; }
static double radToDeg(double rad) { return (rad * 180.0) / M_PI; }
static GeoPoint projectGeoPoint(const GeoPoint& gpt, double distance, double bearing)
    {
        double angular = distance / EARTH_RADIUS;
        double bearingRad = degToRad(bearing);
        double latRad = degToRad(gpt.latitude);
        double lonRad = degToRad(gpt.longitude);

        double lat = std::asin(std::sin(latRad) * std::cos(angular) +
                                       std::cos(latRad) * std::sin(angular) * std::cos(bearingRad));

        double lon = lonRad + std::atan2(sin(bearingRad) * std::sin(angular) * std::cos(latRad),
                                                          std::cos(angular) - std::sin(latRad) * std::sin(lat));

        return  { radToDeg(lon), radToDeg(lat) };
    }

int main()
{
    GeoPoint gp = {0, 0};
    GeoPoint gp_new = MapsMath::projectGeoPoint(gp, 100000, 45);
    std::cout.precision(14);
    std::cout << "longitude_new - " << gp_new.longitude << " latitude_new - " << gp_new.latitude << std::endl;
} 

Result is: longitude_new - 0.63523085373227 latitude_new - 0.6351918164606

I expect that postgis St_Project function will give exatcly same result:

select st_astext(st_project(st_makepoint(0,0,4326), 100000, radians(45)))

However, this query issues different result: longitude_new - 0.635231029125537 latitude_new - 0.639472334729198

What I`m doing wrong and how achieve exactly same result akin postgis St_Project?

  • ST_Project seems to take geography as input postgis.net/docs/ST_Project.html so I believe that it knows that Earth is somewhat flat while your formula probably does not because it has only one constant postgis.net/docs/ST_Project.html. The difference is in sixth decimal place and thus about ten centimeters which is not huge for 100 kilometer distance if you have some practical implementation in your mind. – user30184 Sep 1 '16 at 13:48
  • It could also be the difference between a great circle and a geodesic. – mkennedy Sep 1 '16 at 19:06
  • The second paragraph on the linked page basically says it all. The calculations do not take into account that the world is not a perfect sphere, but a wee bit flattened (an oblate ellipsoid). The other page from the same site may provide a better method – MickyT Sep 4 '16 at 21:51
  • 2
    Digging a little deeper it appears that PostGIS uses this library – MickyT Sep 4 '16 at 22:03
1

PostGIS uses the C implementation of GeographicLib, so I suggest using the native C++ implementation.

ST_Project uses a direct geodesic calculation, e.g. see these docs.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.