[PROJ] Meridian convergence

[Roger Oberholtzer]

Hi


I was looking back through this thread and realize that I did indeed
write the wrong thing. Sigh. The context for this discussion was
Sweref99  - EPSG:3006 - not EPSG:4326 as I had mysteriously wrote. But
the reason for the question was to ensure if it was possible for more
than EPSG:3006. Read on.

I see that when we do this calculation, our code says the following:

      // Seems the EPSG code does not work in this context
       // tinfo->Meridian = proj_create(0, "EPSG:3006");

       // This is the proj5 string for SWEREF99: +proj=utm +zone=33
+ellps=GRS80 +units=m +no_defs +type=crs
       //      Shouldn't that work too?

       tinfo->Meridian = proj_create(0, "proj=tmerc lon_0=15
k_0=0.9996 x_0=500000 ellps=GRS80");


Then, when we want to know the convergence for a given point, we do:

       PJ_FACTORS factors;
       PJ_COORD vals;

       vals.lp.lam = DegreesToRadians(point->LONGITUDE);
       vals.lp.phi = DegreesToRadians(point->LATITUDE);

       factors = proj_factors(tinfo->Meridian, vals);

       convergence = RadiansToDegrees(factors.meridian_convergence));

We then subtract that convergence from the heading provided by the
GPS/IMU and use that when rotating the point.

So, given that we define a specific projection for calling
proj_factors, and it seems not to be possible to use the one we are
using to project our points (i.e., EPSG:3006), the original question
arose.


On Mon, Mar 25, 2024 at 10:10 AM Kristian Evers <kristianevers at gmail.com> wrote:
>
> Roger,
>
> I believe using proj_factors() is the fastest way, yes. The current implementation has its roots in the very early days of PROJ.
> If we were to implement the same functionality from scratch today they would likely be split into more smaller functions with
> one purpose each. You could always implement your own custom version of it based on the PROJ code. Start here:
> https://github.com/OSGeo/PROJ/blob/master/src/factors.cpp
>
> I doubt that your code is particularly useful, to be honest. If your are indeed using EPSG:4326 I believe you will get a meridian
> convergence of 0.0 no matter what the input. The reason being that EPSG:4326 is not a projected CRS and hence you haven’t
> got a grid north that can deviate from true north. Also, EPSG:4326 is not SWEREF99 - it is the WGS84 datum ensemble given
> in geodetic coordinates. You are probably looking for one of the SWEREF99 TM CRS’s.
>
> The meridian convergence calculation should work for most, if not any, projected CRS.
>
> /Kristian
>
> On 25 Mar 2024, at 08.44, Roger Oberholtzer via PROJ <proj at lists.osgeo.org> wrote:
>
> I have been using proj_factors() to get the meridian convergence for a given point. It has been working great. We need this when projecting LiDAR points that involve multiple rotations. The code is something like this:
>
>        PJ_FACTORS factors;
>        PJ_COORD vals;
>
>        vals.lp.lam = DegreesToRadians(point->LONGITUDE);
>        vals.lp.phi = DegreesToRadians(point->LATITUDE);
>
>        factors = proj_factors(tinfo->Meridian, vals);
>
>        return(RadiansToDegrees(factors.meridian_convergence));
>
> I have been using it with Sweref99 (EPSG:4326).
>
> I have two questions:
>
> Is proj_factors() the fastest way to get this information? A look at the code shows that it is calculating other things as well. Is there is a faster way? When generating point clouds, each contains 10s of millions of points, and we generated 10s of thousands of these over a summer. Needless to say, we are always looking for speedups.
>
> Would this method work for all 'standard' projections? By 'standard' I mean the typical ones used in, say, European countries. Would we need to consider anything when extending this to something other than EPSG:4326?
>
> --
> Roger Oberholtzer
> _______________________________________________
> PROJ mailing list
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj
>
>


-- 
Roger Oberholtzer