[PROJ] Calculating meridian convergence

[Roger Oberholtzer]

On Wed, Jun 1, 2022 at 1:27 PM Even Rouault <even.rouault at spatialys.com> wrote:
>
> Roger,
>
> I believe the issue here is that you use a transformation pipeline that
> expects input coordinates to be in degrees.
>
> You should just use your code with
>
> PJ* proj = proj_create(0, "EPSG:3006");

I tried that, and it did not help:

proj_create(0, "EPSG:3006");

LAT 55.600783
LONG 13.100962
angular_distortion 3.141593
meridian_parallel_angle -89.999999
meridian_convergence -1.543445

LAT 58.618118
LONG 15.468616
angular_distortion 3.141593
meridian_parallel_angle -89.999999
meridian_convergence -1.577779

I also tried a suggestion in a different answer here, and it did not help:

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

LAT 55.600783
LONG 13.100962
angular_distortion 0.000000
meridian_parallel_angle 89.999999
meridian_convergence -0.027351

LAT 58.618118
LONG 15.468616
angular_distortion 0.000000
meridian_parallel_angle 90.000000
meridian_convergence 0.006982

Both methods gave different results.

I think the calculation should be

γ = arctan [tan (λ - λ0) × sin φ]

where

λ is the longitude
λ0 is the central meridian of the projection
φ is the latitude

So, an alternative question is: is there a way to find the central
meridian in a PJ CRS? Or wherever I would look for it in proj?

>
> Even
>
> Le 01/06/2022 à 11:49, Roger Oberholtzer a écrit :
> > I am trying to determine the meridian convergence for a point. The
> > projection is defined as:
> >
> > PJ *proj = proj_normalize_for_visualization(0,
> >             proj_create_crs_to_crs(0, "EPSG:4326", "EPSG:3006", 0));
> >
> > The source lat/longs are wgs84 from a GPS receiver. The projection is
> > Sweref99. That is working just fine. In my case, the LONG/LAT are in
> > degrees.
> >
> >         PJ_COORD ivals, vals;
> >         vals.xyz.x = LONG;
> >         vals.xyz.y = LAT;
> >         vals.xyz.z = HEIGHT;
> >         ivals = proj_trans(tinfo->Proj, PJ_FWD, vals);
> >
> > So far, so good.
> >
> > I need to do some rotations of the projected points. Towards that end,
> > I would like to determine the meridian convergence. I am expecting to
> > use this as a correction to a yaw value that my GPS receiver also
> > provides via an integrated IMU.
> >
> >         PJ_FACTORS factors;
> >
> >         vals.lp.lam = proj_torad(LONG);
> >         vals.lp.phi = proj_torad(LAT);
> >
> >         factors = proj_factors(tinfo->Proj, vals);
> >
> > I must be doing something incorrectly, as I get similar results for
> > data from different locations. Or I am missing a step. Or I am just
> > doing something totally wrong.
> >
> > Data from west of the Sweref99 central meridian (convergence should be
> > -1.56713° according to the Swedish authority):
> >
> > LAT 55.600783
> > LONG 13.100962
> > ALT 13.463000
> > E 380351.081178
> > N 6163285.609844
> > angular_distortion 0.567652
> > meridian_parallel_angle 1.570796
> > meridian_convergence -0.004468
> >
> > Data from east of the Sweref99 central meridian (convergence should be
> > 0.40007° according to the Swedish authority):
> >
> > LAT 58.618118
> > LONG 15.468616
> > ALT 79.686000
> > E 527220.410693
> > N 6497625.569591
> > angular_distortion 0.645715
> > meridian_parallel_angle 1.570796
> > meridian_convergence -0.004696
> >
> > angular_distortion, meridian_parallel_angle, and meridian_convergence
> > are the values returned by proj_factors().
> >
> > This is with proj 8.2.1 on Linux. It is generally working fine. It's just this.
> >
> > Any pointers would be greatly appreciated.
> >
> --
> http://www.spatialys.com
> My software is free, but my time generally not.
>


-- 
Roger Oberholtzer