proj_factors produces wrong results with EPSG:2222 (ft)

[jjimenezshaw]

Looking for the behaviour of proj_factors with with a CRS in feet, I found that it returns a huge value

Example of problem

I guess there is no option in proj CLI to pass a CRS identifier. So I modified the file gie_self_tests.cpp. Sorry for the hack.

diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp
index 434822b0..9d7ef11e 100644
--- a/test/unit/gie_self_tests.cpp
+++ b/test/unit/gie_self_tests.cpp
@@ -509,9 +509,11 @@ TEST(gie, info_functions) {
 
     // Test with a projected CRS
     {
-        P = proj_create(PJ_DEFAULT_CTX, "EPSG:3395");
-
+        P = proj_create(PJ_DEFAULT_CTX, "EPSG:2222");
+        a.lp.lam = proj_torad(-110.1666);
+        a.lp.phi = proj_torad(31);
         const auto factors2 = proj_factors(P, a);
+        std::cout << "******" << factors2.meridional_scale << std::endl;
 
         EXPECT_EQ(factors.angular_distortion, factors2.angular_distortion);
         EXPECT_EQ(factors.meridian_parallel_angle,

With this input (using a system in ft) the meridional_scale is 3.51837e+07. However, with EPSG:26948 (the equivalent CRS in m), the escale is 0.9999 (the expected value)

Problem description

The huge meridional_scale value makes no sense. The other factors are nonsense as well.
Debuging the code, I saw that with EPSG:2222 (and NOT with EPSG:26948) it gets into the function fwd_prepare in fwd.cpp. In the line 108

coo.lp.lam = (coo.lp.lam - P->from_greenwich) - P->lam0;

it changes coo.lp.lam from a small value (near zero) to something like 1.9 radians.

I guess that the problem is that it was previously "adjusted" in factors.cpp, line 65

lp.lam -= P->lam0;

Reached this point I do not know how to fix it without breaking anything else.

Expected Output

proj_factors produce meaningful values for any CRS, including CRS in feet.

Environment Information

• PROJ version (proj): master
• Operation System Information: Ubuntu 22.04

Installation method

• compiled

[jjimenezshaw]

Using PR #3825 to use a CRS in proj CLI, I get these results:

$ echo -110.16666666666 31 | ./bin/proj -C EPSG:26948 -V
#Transverse Mercator
#	Cyl, Sph&Ell
#	approx
# +proj=tmerc +lat_0=31 +lon_0=-110.166666666667 +k=0.9999 +x_0=213360 +y_0=0
# +ellps=GRS80
#Final Earth figure: ellipsoid
#  Major axis (a): 6378137.000
#  1/flattening: 298.257222
#  squared eccentricity: 0.006694380023
Longitude: 110d10'W [ -110.16666667 ]
Latitude:  31dN [ 31 ]
Easting (x):   213360.00
Northing (y):  0.00
Meridian scale (h) : 0.99990000  ( -0.01 % error )
Parallel scale (k) : 0.99990000  ( -0.01 % error )
Areal scale (s):     0.99980001  ( -0.02 % error )
Angular distortion (w): 0.000
Meridian/Parallel angle: 90.00000
Convergence : 0d [ 0.00000000 ]
Max-min (Tissot axis a-b) scale error: 0.99990 0.99990

$ echo -110.16666666666 31 | ./bin/proj -C EPSG:2222 -V
#Transformation pipeline manager
# +proj=pipeline +step +proj=tmerc +lat_0=31 +lon_0=-110.166666666667
# +k=0.9999 +x_0=213360 +y_0=0 +ellps=GRS80 +step +proj=unitconvert +xy_in=m
# +xy_out=ft
#Final Earth figure: ellipsoid
#  Major axis (a): 6378137.000
#  1/flattening: 298.257222
#  squared eccentricity: 0.006694380023
Longitude: 110d10'W [ -110.16666667 ]
Latitude:  31dN [ 31 ]
Easting (x):   700000.00
Northing (y):  0.00
Meridian scale (h) : 35183775.10362367  ( 3.518e+09 % error )
Parallel scale (k) : 35183775.10965775  ( 3.518e+09 % error )
Areal scale (s):     1237898030754670.25000000  ( 1.238e+17 % error )
Angular distortion (w): 0.000
Meridian/Parallel angle: 90.00000
Convergence : 125d14'54.806" [ 125.24855720 ]
Max-min (Tissot axis a-b) scale error: 35183775.46019 35183774.75309

[jjimenezshaw]

Also using PR #3825 to use a CRS in proj CLI, I get these results with a CRS that is northing-easting, compared with easting-northing:

echo 9 0 | ./bin/proj -C EPSG:25832 -V
#Universal Transverse Mercator (UTM)
#	Cyl, Ell
#	zone= south approx
# +proj=utm +zone=32 +ellps=GRS80
#Final Earth figure: ellipsoid
#  Major axis (a): 6378137.000
#  1/flattening: 298.257222
#  squared eccentricity: 0.006694380023
Longitude: 9dE [ 9 ]
Latitude:  0dN [ 0 ]
Easting (x):   500000.00
Northing (y):  0.00
Meridian scale (h) : 0.99960000  ( -0.04 % error )
Parallel scale (k) : 0.99960000  ( -0.04 % error )
Areal scale (s):     0.99920016  ( -0.07998 % error )
Angular distortion (w): 0.000
Meridian/Parallel angle: 90.00000
Convergence : 0d [ -0.00000000 ]
Max-min (Tissot axis a-b) scale error: 0.99960 0.99960

echo 9 0 | ./bin/proj -C EPSG:3044 -V
#Transformation pipeline manager
# +proj=pipeline +step +proj=utm +zone=32 +ellps=GRS80 +step +proj=axisswap
# +order=2,1
#Final Earth figure: ellipsoid
#  Major axis (a): 6378137.000
#  1/flattening: 298.257222
#  squared eccentricity: 0.006694380023
Longitude: 9dE [ 9 ]
Latitude:  0dN [ 0 ]
Easting (x):   0.00
Northing (y):  500000.00
Meridian scale (h) : 6375585.74552306  ( 6.376e+08 % error )
Parallel scale (k) : 6375585.74498689  ( 6.376e+08 % error )
Areal scale (s):     -40648093595098.46093750  ( -4.065e+15 % error )
Angular distortion (w): -nan
Meridian/Parallel angle: -90.00000
Convergence : -90d [ -90.00000000 ]
Max-min (Tissot axis a-b) scale error: -nan -nan

[busstoptaktik]

I do not think you should expect any meaningful result in this case: the geometry factor functionality relates to a projection, not to a CRS, so specifying a CRS and asking for factors should be mutually exclusive

[rouault]

I'm on it

[jjimenezshaw]

so specifying a CRS and asking for factors should be mutually exclusive

I do not think so. In fact, I am very interested on the distortion produced in a certain CRS.
This interesting presentation from the NOAA talks about the distortions of the state planes in US: https://geodesy.noaa.gov/web/science_edu/webinar_series/changes-afoot-after-2022.shtml (it includes the horizontal distortion due to the elevation, something not included in proj... yet ;)

[jjimenezshaw]

@rouault Do you know why was it failing? I tried removing one of the adjustments of the longitude (so it is done only once) but still the distortion factors were huge. Computing the derivatives numerically should work.

[rouault]

Do you know why was it failing?

pj_deriv() calls directly the fwd function pointer of a PJ* object. When it is just a single step pipeline of a map projection, this goes directly to the forward map projection method, which doesn't any pre-processing like done in the fwd_prepare() method, whih is called on the different steps by pipeline_forward()
fwd_prepare() also multiplies the x,y values by P->a, which is the reason for the really huge values. And pj_deriv() isn't also ready for the axis swap or unitconvert operations.

[jjimenezshaw]

I see. P->a was the "culprit" (in addition to the adjustment of the longitude twice). The derivative "should" be independent of the units, and those transformations are conformal, so the axis swap is not noticeable.

That means that proj_factors (and subsequently proj -S) do work only with pipelines with a single step. I was thinking on supporting DERIVEPROJECTED, but I guess it is going to be complicated (or impossible).

[busstoptaktik]

I am very interested on the distortion produced in a certain CRS.

The distortion tells something about the relation between two CRS's. Typically between a projected CRS and a geographic base CRS. The scale factor, e.g. is the ratio between (an infinitesimal) distance measured in the first, respectively the second CRS. Talking about distortions is meaningless if not specifying in relation to what.

[jjimenezshaw]

Typically between a projected CRS and a geographic base CRS

That is exactly what #3825 does

[busstoptaktik]

That is exactly what #3825 does

But not exactly what you argued for above