[PROJ] Horner polynomial evaluation

Tue Nov 21 12:24:26 PST 2023

Gabor,

too bad you didn't raise this issue about working with geographic 
coordinates before implementing the workaround. The fix to allow them is 
really trivial: https://github.com/OSGeo/PROJ/pull/3960

Even

Le 21/11/2023 à 20:53, molnar via PROJ a écrit :
> Dear list members, Dear Jochem,
>
> thanks for the detailed description, it would have been very useful 
> for me two years ago, when I prepared a paper, where I used Horner 
> Polynomials for transforming ellipsoidal coordinates to spherical 
> coordinates.
> The oblique Mercator projection used in Hungary (EOV, epsg:23700), is 
> defined in GIS environment as Hotine (or Swiss) Oblique Mercator.  As 
> it is a two-step projection, using the pipeline operator, an 
> alternative method can be defined,that coincides the original 
> definition of the projection: first ellipsoidal coordinates are 
> transformed to spherical coordinates on the Gauss sphere (aposphere) 
> and in the second step the spherical form of the oblique Mercator 
> projection equations are applied on these spherical coordinates to get 
> the projected coordinates.
> The pipeline operator enables us to apply this two step method, and we 
> wanted to use the Horner method, for the first step, to transform the 
> ellipsoidal coordinates to spherical coordinates. (The coefficients 
> were published in the standardization document, published in 1972, in 
> the calculator era).
> The manual of Horner Polynomial Evaluation on proj.org states, that 
> both input and output coordinates can be either geodetic, either 
> projected, but we could not apply the method on geodetic coordinates 
> in a cct pipeline.
> As a roundabout process, in the pipeline,
> 1:
> the input geodetic coordinates were converted to fake 'projected' 
> coordinates, using the
> +step +proj=eqc +a=57.295779513082323
> that literally does not change(!) the figures, but proj regards the 
> output as projected coordinates. The equidistant cylindrical 
> projection with +a=57.29... Earth radius ('a' is actually 1 radian in 
> degrees) converts 47 degree North latitude and 19 degree East 
> longitude to 47 meter Northings, and 19 meter Eastings...
> 2:
> We used Horner Polynomial Evolution to modify the coordinates. Proj 
> believes that it transforms projected coordinates to projected 
> coordinates.
> 3:
> We used the inverse projection, to get ellipsoidal (actually 
> spherical) coordinates:
> +step +proj=eqc +a=57.295779513082323 +inv
> and we finally used the spherical form of the oblique Mercator 
> projection.
>
> The details are in the paper:
> https://www.fig.net/resources/proceedings/fig_proceedings/fig2022/papers/ts01c/TS01C_molnar_toth_11649.pdf 
>
>
> Regards, Gabor
>
>
>
> On 2023-11-21 16:23, Lesparre, Jochem via PROJ wrote:
>> Dear list members,
>>
>> The implementation of +proj=horner is nice. I used it to reproduce a
>> 99 year old polynomial transformation of an even older Dutch CRS
>> published in 1827.  However, I encountered some problems with the
>> documentation on
>> https://proj.org/en/9.3/operations/transformations/horner.html. The
>> order of the coefficients in the formulas and in the description of
>> the PROJ parameters (+fwd_u and +fwd_u) seem mutually inconsistent.
>>
>> PROBLEM 1: The matrix formula (after formula 3) seems inconsistent
>> with formula (1) and the actual implementation.
>>
>> \begin{bmatrix}
>>
>>      u_{0,1} + u_{0,2} U + ... & u_{1,0} + u_{1,1} U + u_{2,0} V + ...
>> \\
>>
>>      v_{1,0} + v_{1,1} V + v_{2,0} U + ... & v_{0,1} + v_{0,2} V \\
>>
>> \end{bmatrix}
>>
>> should be changed to:
>>
>> \begin{bmatrix}
>>
>>      u_{1,0} + u_{2,0} U + ... & u_{0,1} + u_{1,1} U + u_{0,2} V + ...
>> \\
>>
>>     v_{0,1} + v_{1,1} V + v_{0,2} U + ... & v_{1,0} + v_{2,0} V \\
>>
>> \end{bmatrix}
>>
>> PROBLEM 2: The description of the parameters seems inconsistent with
>> formula (1) and the actual implementation.
>>
>> +fwd_u=<u_11,u_12,...,u_ij,..,u_mn>
>>
>> and
>>
>> +fwd_v=<v_11,v_12,...,v_ij,..,v_mn>
>>
>> should be changed to:
>>
>> +fwd_u=<u_00,u_10,...,u_ij,..,u_nn>
>>
>> and
>>
>> +fwd_v=<v_00,v_01,...,v_ij,..,v_nn>
>>
>> As a result it was difficult to determine the correct order of the
>> coefficients for +fwd_u and +fwd_v parameters. I think it would be
>> clarifying to add a table for the order of the coefficients in the
>> documentation (example for +deg=3):
>>
>> Order of coefficients of +fwd_u:
>>
>>     | x^0 | x^1 | x^2 | x^3
>>
>> y^0 | 1   | 2   | 3   | 4
>>
>> y^1 | 5   | 6   | 7   | -
>>
>> y^2 | 8   | 9   | -   | -
>>
>> y^3 | 10  | -   | -   | -
>>
>> Order of coefficients of +fwd_v:
>>
>>     | x^0 | x^1 | x^2 | x^3
>>
>> y^0 | 1   | 5   | 8   | 10
>>
>> y^1 | 2   | 6   | 9   | -
>>
>> y^2 | 3   | 7   | -   | -
>>
>> y^3 | 4   | -   | -   | -
>>
>> SUGGESTION FOR A DIFFERENT IMPLEMENTATION: Actually, I would have
>> preferred a different order of the coefficients in the implementation,
>> in a way that the order would be the same for +fwd_u and +fwd_v, and
>> it would keep the first coefficients the same when increasing the
>> degree (+deg) of the polynomial:
>>
>> Order of coefficients of +fwd_u and +fwd_v:
>>
>>     | x^0 | x^1 | x^2 | x^3
>>
>> y^0 | 1   | 2   | 4   | 7
>>
>> y^1 | 3   | 5   | 8   | -
>>
>> y^2 | 6   | 9   | -   | -
>>
>> y^3 | 10  | -   | -   | -
>>
>> Regards, Jochem
>>
>> Disclaimer:
>> De inhoud van deze e-mail is vertrouwelijk en uitsluitend bestemd voor
>> de geadresseerde(n).
>> Gebruik, openbaarmaking, vermenigvuldiging, verspreiding en/of
>> verstrekking van deze informatie aan derden is niet toegestaan.
>> Op al onze producten en diensten zijn onze algemene
>> leveringsvoorwaarden van toepassing
>> [https://www.kadaster.nl/algemene-leveringsvoorwaarden].
>>
>> Disclaimer:
>> This email and any files transmitted with it are confidential and
>> intended solely for the use of the individual or entity to whom they
>> are addressed.
>> If you are not the intended recipient, you are notified that
>> disclosing, copying, distributing or taking any action in reliance on
>> the contents of this information is strictly prohibited.
>> Our general terms and conditions of delivery apply to all our products
>> and services
>> [https://www.kadaster.com/general-terms-and-conditions].
>> _______________________________________________
>> PROJ mailing list
>> PROJ at lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/proj
> _______________________________________________
> PROJ mailing list
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj

-- 
http://www.spatialys.com
My software is free, but my time generally not.