[PROJ] Projection type required for proj_lp_dist

[Kristian Evers]

Thanks, Even! I’ll take care of the docs later tonight or tomorrow.

/Kristian

> On 10 Mar 2021, at 19:07, Even Rouault <even.rouault at spatialys.com> wrote:
> 
> Fix pending in https://github.com/OSGeo/PROJ/pull/2570
> 
> Documentation of those functions would deserve enhancements though. See my comment in the PR
> 
> Le 10/03/2021 à 16:39, Kristian Evers a écrit :
>> 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.
>>> 
>> _______________________________________________
>> PROJ mailing list
>> PROJ at lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/proj
> 
> -- 
> http://www.spatialys.com
> My software is free, but my time generally not.
>