proj_factors() only seems to work with a proj string

[ojwb]

Example of problem

Patching the PROJ test suite as below demonstrates the issue:

diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp
index 5fc8f0f1..2ef339cc 100644
--- a/test/unit/gie_self_tests.cpp
+++ b/test/unit/gie_self_tests.cpp
@@ -454,7 +454,7 @@ TEST(gie, info_functions) {
                 1e-7);
 
     /* test proj_derivatives_retrieve() and proj_factors_retrieve() */
-    P = proj_create(PJ_DEFAULT_CTX, "+proj=merc +ellps=WGS84");
+    P = proj_create(PJ_DEFAULT_CTX, "EPSG:3395");
     a = proj_coord(0, 0, 0, 0);
     a.lp.lam = proj_torad(12);
     a.lp.phi = proj_torad(55);

Problem description

With the patch above applied the PROJ test suite fails like so:

$ make check

[....]

[ RUN      ] gie.info_functions
Unrecognized vertical grid format for filename 'proj.db'
Unrecognized horizontal grid format for filename 'proj.db'
Invalid latitude or longitude
gie_self_tests.cpp:463: Failure
Value of: proj_errno(P)
  Actual: true
Expected: false
[  FAILED  ] gie.info_functions (2 ms)

More generally I can't seem to get proj_factors() to work in any case unless I've called proj_create() or similar with a PROJ string.

This is a problem because current PROJ docs say:

The use of proj-string to describe a CRS is discouraged. It is a legacy means of conveying CRS descriptions: use of object codes (EPSG:XXXX typically) or WKT description is recommended for better expressivity.

It appears I'm currently prevented from following this advice if I want to use proj_factors().

Expected Output

I'd expect your test suite to still pass with the above patch applied (assuming that EPSG:3395 is actually equivalent to +proj=merc +ellps=WGS84 - it appears to be to me at least).

Environment Information

• PROJ version: git master as of 47fb85b
• Operation System Information: Debian unstable amd64

Installation method

I'm testing as above in a clone of your git repo.

(Issue originally encountered with debian proj package version 7.2.1-1 - I checked then with git master to avoid reporting if it was already fixed.)

[kbevers]

That's definitely not intended. The reason you are seeing this behaviour is that proj_factors() was a fairly early addition to the proj.h API which basically just bridged the gap between the old projects.h and the new proj.h. The function should either be removed or adapted to fit the more advanced functionality we have today.

[ojwb]

I'd certainly prefer it not to be removed.

I'm using it to get the grid convergence (and nothing else) to be able to convert magnetic bearings to bearings relative to grid north (in combination with the IGRF magnetic model).

I could essentially implement what pj_factors() and pj_deriv() do internally to calculate this in my own code, but that means others will also have to implement this themselves. I know at least grass currently uses PROJ to get meridian_convergence, and also seems to only uses pj_factors()/proj_factors() for that.

Calculating this numerically myself also means missing out when there's an analytical solution (the docs for proj_factors() say "Depending on the underlying projection values will be calculated either numerically (default) or analytically" - it looks to me like currently it's actually always done numerically, but presumably it could be done analytically for some projections in a future PROJ version).

Incidentally, if this API does get reworked, it seems odd that it still takes radians when in general PROJ now seems to work in degrees for latitude and longitude. If it does stay using radians, it would be helpful to note this point in the API docs.

[busstoptaktik]

(the docs for proj_factors() say "Depending on the underlying projection values will be calculated either numerically (default) or analytically" - it looks to me like currently it's actually always done numerically, but presumably it could be done analytically for some projections in a future PROJ version)

The docs are wrong. Originally, there were analytical solutions for, if I remember correctly, two projections, and no snowball's chance in hell any other were to join the exclusive team of two. For these two, an entire sub-component of the source tree was maintained (or rather: Was left unmaintained for 25 years), for the dubious benefit of having an alternative way of computing virtually bit-identical values.

So while cleaning up PJ_FACTORS (some time around PR #683), I removed all traces of analytical solutions: They complicated the code for no real gain, and had the unfortunate side effect of overshadowing the numerical solutions on the API-side of things, making it hard to test for consistency.

Testing for consistency is important, because the numerical version is guaranteed to return the characteristics of the actually implemented code, while the analytical stuff was an entirely different branch of the code, which could return anything.

So all-in-all: Much complication for no real gain, testing headaches, and no sign of interest from anyone for maintaining and extending the analytical derivatives code - that made it an easy decision to axe that part of the code.

But a pull-request correcting the documentation would be most welcome :-)

[busstoptaktik]

@kbevers opined:

That's definitely not intended. The reason you are seeing this behaviour is that proj_factors() was a fairly early addition to the proj.h API which basically just bridged the gap between the old projects.h and the new proj.h. The function should either be removed or adapted to fit the more advanced functionality we have today.

This is actually an interesting case, neatly illustrating some serious shortcommings of the representation of geodetic coordinate systems in the ISO/OGC geoinformatics model, but also some tripwires in general nomenclature.

Essentially, the code works as intended: it should neither be amended nor removed, and Olly's expected output should not be expected. But it will take quite a few words to explain why...

@ojwb said about the expected output:

I'd expect your test suite to still pass with the above patch applied (assuming that EPSG:3395 is actually equivalent to +proj=merc +ellps=WGS84 - it appears to be to me at least).

For a start, let's try to state the case in a nutshell: No, the coordinate reference system EPSG:3395 is not equivalent to the coordinate operation proj=merc ellps=WGS84. They are two fundamentally different things, and non-equivalent in the same sense as apples and gravity are non-equivalent.

Confusing the two is understandable, however, especially since the plain text name of EPSG:3395 is "WGS84 / World Mercator", making it obvious to believe that the WGS84 of EPSG:3395 is the same as the WGS84 of proj=merc ellps=WGS84. But the former is the datum ensemble WGS84, while the latter is the identically named reference ellipsoid. Again, two different things.

The confusion of these matters is encouraged by the mistaken foundation of the ISO/OGC geospatial standards series, which somehow asserts that a CRS is definable in absolute terms. Being able to define a CRS in absolute terms would be nice, since once you have an absolute definition of two CRS', you would be able to determine an infallible transformation between the two.

That's possible in mathematics, where coordinate systems are Platonic ideals. In geodesy, coordinate systems are much more messy: You can only define a reference system, by writing a book describing how to realize that system on the physical Earth.

The reference system is the book. The realization is the corresponding reference frame. And the reference frame (i.e. a collection of physical points with associated coordinate and velocity information) is what you can measure point coordinates with respect to.

So the definition (i.e. the book) may guide us toward constructing a transformation involving a given CRS/reference frame. But we cannot determine any coordinates of physical features with respect to the system - only with respect to the frame.

The projected CRS EPSG:3395 is related to the geographical base CRS EPSG:4326 by the coordinate operation described by proj=merc ellps=WGS84 (or actually its inverse form). So the closest you can get to a "definition" of EPSG:3395 is to say that EPSG:3395 is the CRS for which coordinates gets related to EPSG:4326 by subjecting them to the coordinate operation given by inv proj=merc ellps=WGS84. Or in other words: The definition of a CRS is the coordinate operation which brings us to its base CRS. Once we arrive at the base CRS, that's the end of the definition in absolute terms: You have arrived from your trip from the platonic gardens of mathematics to the messy moors of geodesy.

To continue the journey from there and onto another base CRS, you will have to rely on empirically determined transformations - you are in the waste lands of approximations, where a meter is not a meter, a radian not a radian, and the distance between two points is not the same as the difference between their coordinates. Welcome to geodesy!

In one sense, however, things are much simpler in geodesy: A CRS is really just a label, without any internal state. While that sounds strange comming from the ISO/OGC world, it really simplifies a lot of things, since that label is the key to looking up the transformation to any other CRS in the little black book all geodesists are secretly given upon graduation (or, having lost the book: Looking it up at the EPSG website).

Now back to proj_factors(): What you should remember here is that proj_factors() provides information about the geometric properties of a projection, i.e. a purely mathematical concept. The projection is given with respect to purely mathematical geographical coordinates which in turn are given with respect to a purely mathematical ellipsoid. The messy details of geodesy is left out here.

Obviously, by feeding proj_factors() a diet of sufficiently perverse PROJ pipelines, you will be able to force the geodetic mess onto it, but even in that case, proj_factors() tells you about the geometric properties of a coordinate operation: It has nothing to say about the properties of a CRS (although, from a geodetic point of view, a CRS and the coordinate operation relating coordinates to it, are closely related).

So all-in-all: proj_factors() works as expected, and I will close this issue in a few days unless @ojwb or @kbevers have additional input on the subject.

[rouault]

I guess it could be reasonable to accept a PJ* object it is a Projected CRS and to use the conversion from its base CRS as the PROJ pipeline. The question is: what about axis order and units of the input ? Should we use long / lat in radians as with a raw "+proj==merc +ellps=WGS84", or lat / long in degrees as implied by the base CRS of EPSG:3395 being EPSG:4326

[busstoptaktik]

The question is: what about axis order and units of the input ?

Yes, that's a huge question, which I think needs input from actual use cases, so perhaps @ojwb can advise about what would be expected?

[snowman2]

pyproj exports the CRS to a PROJ string ref

And it converts from radians to degrees for the user ref

Since you are exporting to a PROJ CRS string, the axis order will always be x,y correct?

[rouault]

Since you are exporting to a PROJ CRS string, the axis order will always be x,y correct?

if using the conversion attached to the projected CRS, then the "natural" pipeline would be the one returned by projinfo -s EPSG:4326 -t EPSG:3395 which is

+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84

that it it expects input to be in lat, long order and degree.

But yes, if using instead the PROJ.4 string of the CRS, we have "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs", so if we remove the +type=crs, then it would be a valid pipeline with expected input in long, lat order and radians

[snowman2]

That would be nice! The proposed workflow would be:

>>> from pyproj import CRS, Transformer
>>> crs = CRS("EPSG:3395")
>>> transformer = Transformer.from_crs(crs.geodetic_crs, crs)
>>> transformer.definition
'proj=pipeline step proj=axisswap order=2,1 step proj=unitconvert xy_in=deg xy_out=rad step proj=merc lon_0=0 k=1 x_0=0 y_0=0 ellps=WGS84'

This is always guaranteed to only have one transformation with a definition, correct?

[ojwb]

Obviously, by feeding proj_factors() a diet of sufficiently perverse PROJ pipelines, you will be able to force the geodetic mess onto it, but even in that case, proj_factors() tells you about the geometric properties of a coordinate operation: It has nothing to say about the properties of a CRS (although, from a geodetic point of view, a CRS and the coordinate operation relating coordinates to it, are closely related).

If I follow what you're saying then this should work shouldn't it?

diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp
index 2ef339cc..080a0c6e 100644
--- a/test/unit/gie_self_tests.cpp
+++ b/test/unit/gie_self_tests.cpp
@@ -454,7 +454,7 @@ TEST(gie, info_functions) {
                 1e-7);
 
     /* test proj_derivatives_retrieve() and proj_factors_retrieve() */
-    P = proj_create(PJ_DEFAULT_CTX, "EPSG:3395");
+    P = proj_create_crs_to_crs(PJ_DEFAULT_CTX, "EPSG:4326", "EPSG:3395", nullptr);
     a = proj_coord(0, 0, 0, 0);
     a.lp.lam = proj_torad(12);
     a.lp.phi = proj_torad(55);

But that also fails.

So all-in-all: proj_factors() works as expected, and I will close this issue in a few days unless @ojwb or @kbevers have additional input on the subject.

Are you saying proj_factors() is really only expected to work if a proj string is used?

If not, then this issue seems valid even if my initial reduced test case isn't.

Yes, that's a huge question, which I think needs input from actual use cases, so perhaps @ojwb can advise about what would be expected?

My use case is getting the grid convergence so I can convert magnetic bearings to be relative to grid north (in combination with the IGRF model which gets me from magnetic to true north).

I know the coordinates of the point in question in some source CRS and also have them already converted to WGS84 (as that's what the IGRF model takes as input). And there's a target CRS (which can be the same or different to the source CRS) I want the grid convergence for.

At the moment I feed proj_factors() a PJ from proj_create() passing the proj string for the target CRS and WGS84 lat/long (in radians).

I don't really care what the units or coordinate order are so long as they are clearly documented, though it seems a bit odd that radians are required here still. Working from the source CRS coordinates would be fine too (though I'll still need to convert them to WGS84 for IGRF anyway).

I've ported the rest of the codebase away from the old PROJ API, and this is the remaning blocker.

[busstoptaktik]

Are you saying proj_factors() is really only expected to work if a proj string is used?

More or less: Originally, the code was designed to work with "classic PROJ", i.e. say prior to 4.9, which supported nothing but basic projections. Once you mix EPSG-codes into the cocktail, entirely different parts of the code base gets invoked under the hood (I'll leave commenting on that part to @rouault who designed it). So feed it a classic PJ object (most probably including a handcrafted pipeline, although I have never had the need for that), and it should work. Otherwise, probably not, although @rouault describes a possible road to implementation above.

But I repeat: proj_factors() is not designed to tell you anything about a CRS, but about the projection implicitly defining the CRS in terms of its base-CRS. To proj_factors() the world is smooth, and does not include any of the geometrical bumps on the road a practically realized CRS is (and especially was, for older systems) born with.

though it seems a bit odd that radians are required here still

The PROJ string parameterizes the PROJ-internal operators, which work in radians. The question of units is a matter of definition for the CRS at hand (although they will typically be degrees or gons/grads), not for the internal operations, which work in radians for practical reasons (since typically projection code involves a large number of trigonometric functions, and the C standard library expects radians as input for those).

[busstoptaktik]

If I follow what you're saying then this should work shouldn't it?

Probably not, since what you get out from the proj_create_crs_to_crs() takes you from EPSG:4326, which is latitude/longitude in degrees, while proj_factors() iirc will expect input in "GIS order" of longitude/latitude in radians.

[snowman2]

If I follow what you're saying then this should work shouldn't it?

Not yet. Changes are required for this to work.