proj_factors(): accept P to be a projected CRS (fixes #2854)

[rouault]

Updated doc:

Starting with PROJ 8.2, the P object can be a projected CRS, for example
instantiated from a EPSG CRS code. The factors computed will be those of the
map projection implied by the transformation from the base geographic CRS of
the projected CRS to the projected CRS.

The input geodetic coordinate lp should be such that lp.lam is the longitude
in radian, and lp.phi the latitude in radian (thus independently of the
definition of the base CRS, if P is a projected CRS).

CC @busstoptaktik

[busstoptaktik]

I have no objections to this - as I have stated in the discussion in #2854 I believe the code "as is" is not broken; I do, however recognize the use cases presented by the participants in the discussion, and laid out especially clearly by @snowman2, so I think this patch holds merit and am happy to endorse it, and hope it fulfills the needs of the intended audience, although it is not entirely clear to me whether this will be the case: The fixed lp.lam/lp.phi interpretation of the input may be a problem - I cannot comment on that: To me at least, this looks good!

[rouault]

The fixed lp.lam/lp.phi interpretation of the input may be a problem - I cannot comment on that: To me at least, this looks good!

yeah I was unsure about the best approach for that. I'm afraid there's no clear winner. I also considered the option to use the axis unit and order of the base geographic CRS of the projected CRS, as expected with proj_trans().
That said we do have somehow a problem with the PJ_COORD structure. It was thought prior to all the PROJ 6 work, where the axis order of the authority wasn't taken into account, and if you for example do proj_trans() from EPSG:4326 to something else, you need to provide input as latitude, longitude degree. But in which member of the PJ_COORD structure to do that ? A "naive" user might be tempted to use the lp / lpz / lpzt, but he would get things wrong as the implied semantics of lp.lam and lp.phi is infact completely lost when you just get a PJ_COORD. So the most reasonable is to use the xyzt or v[4] members.

[busstoptaktik]

That said we do have somehow a problem with the PJ_COORD structure. It was thought prior to all the PROJ 6 work, where the axis order of the authority wasn't taken into account

Yes - over in https://github.com/busstoptaktik/geodesy I try a different approach, by simply calling the coordinates "first, second, third and fourth", i.e. the v[4] solution, and proclaiming that internally all operators work in "gis" order with radians for angular units, then having "adaptors" handling the ormat of input and output systems in both ends of pipelines. In Rust Geodesy, a coordinate-tuple is "just a 4-tuple", its interpretation is up to the operator.

I think it's fine to merge this, but perhaps @snowman2 or @ojwb would have an opinion on this?

[ojwb]

Thanks for working on this.

How the coordinates are specified here probably makes it simpler for my existing code, but I can certainly see the point that it's not entirely consistent with the modern PROJ approach. It's clearly and explicitly documented here at least.

I'll aim to find time to try it out over the weekend and see if I have any feedback in the light of actually trying to use it.

[ojwb]

I've now tested and this branch does indeed solve my problem with getting convergence values from PROJ, so thanks for that.

I don't think I really have any further useful insights about the coordinate order and units. They seems OK from my code, though that was written already.