[PROJ] Projection type required for proj_lp_dist

[Kristian Evers]

Sorry, I thought your problem was getting set up before the call to proj_lp_dist. 
I’ve tried to reproduce your problem and I believe you’ve come across a bug.
Below is a bit of test code I wrote to demonstrate the behaviour. As far as I can tell
you get similar results.


#include <stdio.h>
#include "proj.h"

int main() {

    PJ_CONTEXT *C = proj_context_create();
    PJ *P;
    //P = proj_create(C, "EPSG:4258");                                  // returns inf
    //P = proj_get_ellipsoid(C, proj_create(C, "EPSG:4258"));           // returns inf
    //P = proj_create(C, "+proj=latlong +ellps=GRS80");                 // returns 12836.740931
    //P = proj_create(C, "+proj=latlong +ellps=GRS80 +type=crs");       // returns inf
    //P = proj_create(C, "+proj=utm +zone=32 +ellps=GRS80");            // returns 12836.740931
    //P = proj_create_crs_to_crs (C,"EPSG:4326", "EPSG:4258", NULL);    // returns 12836.740931
    P = proj_create_crs_to_crs (C,"EPSG:4326", "EPSG:32632", NULL);     // returns 12836.740931

    PJ_COORD c1 = proj_coord(proj_torad(12.0), proj_torad(55.0), 0, 0);
    PJ_COORD c2 = proj_coord(proj_torad(12.1), proj_torad(55.1), 0, 0);

    double d = proj_lp_dist(P, c1, c2);

    printf("Distance: %f", d);

    proj_destroy(P);

    return 0;

}


I would expect all of the above (perhaps with the exception of the proj_get_ellipsoid
line) to return the same distance.

This is most likely because internally P->geod is not initialized when creating a CRS object.
I haven’t been able to pinpoint exactly where but someone more familiar with that part of the
code might be able to do so quickly.



> On 10 Mar 2021, at 15:27, Jason Dick <JasonD at TheTranstecGroup.com> wrote:
> 
> On 3/9/2021 5:06 PM, Kristian Evers wrote:
>> Try taking a look at this function from the testsuite:
>> https://github.com/OSGeo/PROJ/blob/f278b5bee641dd94245ed54f4eb75e29c3d4d993/test/unit/test_c_api.cpp#L166-L203
> Thanks for the pointer, but it's not getting through to me. From what I gather, Infinity is the correct answer when using EPSG:4326, but it doesn't explain why. Why does distance require an operation PJ and not a CRS? I've looked at the source code for the distance function before, but it wasn't helpful either. Keep in mind, I have no geo background. I just use libraries for transforming coordinates. I've learned a lot, but PROJ is much more "technical" than what I've used in the past.
>