patchGNU Octave - Patches: patch #10022, [octave forge] (mapping) Auxillary...

 
 

patch #10022: [octave forge] (mapping) Auxillary Latitudes

Submitter:  None
Submitted:  Wed 17 Feb 2021 06:56:44 AM UTC
   
 
Category:  Forge : new function Priority:  5 - Normal
Status:  Done Privacy:  Public
Assigned to:  philipnienhuis Originator Email:  -email is unavailable-
Open/Closed:  Closed
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Tue 01 Feb 2022 09:40:47 PM UTC, comment #20: 

I've pushed the adapted geo2iso.m here:
http://hg.code.sf.net/p/octave/mapping/rev/86ccecdf46c9

and iso2geo.m here:
http://hg.code.sf.net/p/octave/mapping/rev/ca467157372a

I've adapted Texinfo, style, simplified input checks (esp. ellipsoid, as we now have one function for that that does it all), adapted tests and examples (values) and some more things.

I am very pleased with your contributions, Anonymous. There was some work to do for me there but I enjoyed it (hobby after all) and I enjoyed working with you.

Thanks again!

Closing report.

Philip Nienhuis <philipnienhuis>
Group Member
Tue 01 Feb 2022 09:05:42 PM UTC, comment #19: 

Thanks very much.
Now I get:

>> iso2geo (geo2iso ([-90 : 10 : 90]))
ans =
 Columns 1 through 17:
  -90  -80  -70  -60  -50  -40  -30  -20  -10    0   10   20   30   40   50   60   70
 Columns 18 and 19:
   80   90


which looks good!
But:
I'm astounded because of this in geo2iso:

:
  elseif (strncmpi (angleUnit, "degrees", min (length (angleUnit), 7)))
:   phi = deg2rad (phi);

:
  chi = geo2con (phi);
:


where in geo2iso.m:

  • I'd expect the call 2 geo2con to also contain input args E and "radians", as geo2con also needs the ellipsoid and higher up, phi was morphed into radians;
  • geo2con.m expects "degrees" by default.


If I change the call 2 geocon() into

  chi = geo2con (phi, E, "radians");


I get:

>> geo2iso ([-90 : 10 : 90], "wgs84", "degrees")
ans =
 Columns 1 through 6:
         -Inf  -1.3921e+02  -9.9071e+01  -7.5123e+01  -5.7614e+01  -4.3465e+01
 Columns 7 through 12:
  -3.1281e+01  -2.0288e+01  -9.9846e+00  -6.3611e-15   9.9846e+00   2.0288e+01
 Columns 13 through 18:
   3.1281e+01   4.3465e+01   5.7614e+01   7.5123e+01   9.9071e+01   1.3921e+02
 Column 19:
          Inf


which is the exact same output as Matlab.

So the new iso2geo.m still yields erroneous values. But that is easily solved as in iso2geo you likewise need to feed args 'E, "radians"' to the input args of con2geo() ...

OK I'll polish style and Texinfo and push geo2iso.m & iso2geo.m later on. With that done it's time for a new mapping pkg release.

Philip Nienhuis <philipnienhuis>
Group Member
Mon 31 Jan 2022 10:33:27 AM UTC, comment #18: 

Cleaner version of geo2iso.
New version of iso2geo.

(file #52781, file #52782)

Anonymous
Wed 05 Jan 2022 08:31:13 PM UTC, comment #17: 

Thanks.
I know little about these latitude conversions; I just checked if **2geo & geo2** were each other's inverse.
Also, I tend to think of Matlab as maybe authorative in principle  but no guarantees it is absolutely correct :-)

Apart from having a fixed iso2geo, I want to replace the various reference ellipsoid checks in the code by sph_chk. In the mean time I'm busy with fixes for shaperead and shapedraw.
Once all that has been done it's definitely time for a new mapping package release.

Philip Nienhuis <philipnienhuis>
Group Member
Wed 05 Jan 2022 02:01:34 PM UTC, comment #16: 

Started looking.

Because the source for iso2geo algorithm was known.
Stared with geo2iso

L82 came from
https://en.wikipedia.org/wiki/Latitude#Isometric_latitude

Following Synder pg [15]

psi = log(tan(pi/4+phi/2)*((1-ecc.*sin(phi))/(1+e*sin(phi)))^(ecc/2));


Which gives


>> geo2iso(0:10:90)
ans =

     -2.72      7.33     17.70     28.75     40.99     55.19     72.73     96.71    136.86   2136.24

>> iso2geo (geo2iso (0:10:90))
ans =

   -2.74    7.36   17.54   27.80   38.14   48.54   59.01   69.51   80.05   90.61


In Synder pg [15] it also states

chi = geo2con(phi);
psi = log ( tan(pi/4 + chi/2) );


could be used but also states

although there is no particular advantage over using the previous equation


The results are

>> geo2iso(0:10:90)
ans =

    -0.00     9.98    20.28    31.24    43.36    57.39    74.66    98.08   136.58   300.72

>> iso2geo (geo2iso (0:10:90))
ans =

       0   10.00   20.00   30.00   40.00   50.00   60.00   70.00   80.00   90.00



Next will look at iso2geo


Anonymous
Sun 02 Jan 2022 09:10:19 PM UTC, comment #15: 

Sorry for long delay.

In the meantime I checked with Matlab and AFAICS con2geo.m/geo2con.m, auth2geo.m/geo2auth.m and rect2geo.m/geo2rect.m seem to work fine. They've already been pushed, as was geo2iso.m
However, iso2geo seems to yield incorrect answers.

According to Matlab docs and tried here in r2020b:


>> phi = [-90 : 10 : 90];
conv = map.geodesy.IsometricLatitudeConverter(wgs84Ellipsoid);
psi = forward(conv,phi)   % geodetic to isometric
psi =
  Columns 1 through 9
      -Inf   -2.4296   -1.7291   -1.3112   -1.0055   -0.7586   -0.5460   -0.3541   -0.1743
  Columns 10 through 18
         0    0.1743    0.3541    0.5460    0.7586    1.0055    1.3112    1.7291    2.4296
  Column 19
       Inf
>> psi = rad2deg (forward(conv,phi))
psi =
  Columns 1 through 9
      -Inf -139.2081  -99.0708  -75.1234  -57.6137  -43.4647  -31.2810  -20.2878   -9.9846
  Columns 10 through 18
         0    9.9846   20.2878   31.2810   43.4647   57.6137   75.1234   99.0708  139.2081
  Column 19
       Inf


and conversely (using inverse rather than forward).

Your geo2iso does it fairly well, apart from the Inf extremes:

>> format bank     % to avoid scientific notation
>> pssi = geo2iso ([-90 : 10 : 90])
pssi =
 Columns 1 through 8:
  -2178.29   -139.21    -99.07    -75.12    -57.61    -43.46    -31.28    -20.29
 Columns 9 through 16:
     -9.98         0      9.98     20.29     31.28     43.46     57.61     75.12
 Columns 17 through 19:
     99.07    139.21   2178.29


Before invoking iso2geo the latitude interval checks have to be dropped, and then:

>> iso2geo (geo2iso (0:10:90))
ans =
       0   10.00   20.01   30.03   40.08   50.15   60.24   70.35   80.47   90.61


... so iso2geo and geo2iso don't seem to be exactly each other's inverse. AFAICS it's esp. iso2geo.m which doesn't seem to be sufficiently accurate.
It was already months ago that I first noted this but I simply hadn't time to delve deeper into it, sorry.

Philip Nienhuis <philipnienhuis>
Group Member
Fri 03 Dec 2021 09:18:53 PM UTC, comment #14: 

Not specifically.
I merely absorbed the patches you kindly supplied into the mapping package. The last one (iso2geo) is pending. So I suppose all coordinate / geodesy stuff is about complete, right?

The wiki might need an update. Several functions you supplied are still marked as "unimplemented". If you want you can do that.

A next thing I'd like to have in mapping is projections and projection transformations; I have some basic code invoking the proj package as backend but presently lack time.

Philip Nienhuis <philipnienhuis>
Group Member
Fri 03 Dec 2021 02:35:48 PM UTC, comment #13: 

In the wiki it states this:

The current roadmap is to further integrate Felipe Nievinsky's geodesy toolbox, updated by M. Hirsch.

-----
However M. Hirsch's functions have been included.
Are there certain functions from Felipe Nievinsky's work you are looking for?

Anonymous
Thu 21 Oct 2021 06:18:42 PM UTC, comment #12: 

Thank you for adding these:

From comment #9 found the issue.
Attached the hg diff for con2geo.

Also iso2geo is attached.



(file #52141, file #52142)

Anonymous
Sun 17 Oct 2021 07:17:57 PM UTC, comment #11: 

As to comment #7, hold on, I'll retract that comment, I chose the wrong converter.

I've pushed con2geo.m & auth2geo.m here:
http://hg.code.sf.net/p/octave/mapping//rev/bae1ee23f10c

with similar edits as for rect2geo.

Will you also supply iso2geo?

Philip Nienhuis <philipnienhuis>
Group Member
Sun 17 Oct 2021 06:52:45 PM UTC, comment #10: 

I pushed rect2geo here:
http://hg.code.sf.net/p/octave/mapping//rev/5d625af9eea8

I edited style (Texinfo & code) a bit.
Note that your rect2geo.m only processed scalars or row vectors properly, other shapes resulted in scalar (and incorrect) output. Where you had this:

    A = [mu;
        ((3 / 2) * n - (27 / 32) * n3) * sin(2 * mu);
        ((21 / 16) * n2 - (55 / 32) * n4) * sin(4 * mu);
        ((151 / 48) * n3) * sin(6 * mu);
        ((1097 / 512) * n4) * sin(8 * mu)];
    phi = sum (A);


... I changed it into:

  phi = mu .+ ...
        (((3 / 2) * n - (27 / 32) * n3) * sin (2 * mu)) .+ ...
        (((21 / 16) * n2 - (55 / 32) * n4) * sin (4 * mu)) .+ ...
        (((151 / 48) * n3) * sin (6 * mu)) .+ ...
        (((1097 / 512) * n4) * sin (8 * mu));


... that way rect2geo will accept about any shape of input.

Thanks so far.

Philip Nienhuis <philipnienhuis>
Group Member
Sun 17 Oct 2021 06:44:26 PM UTC, comment #9: 

I compared con2geo with Matlab r2021b:

>> spheroid = referenceEllipsoid ('wgs84');
>> converter = map.geodesy.ConformalLatitudeConverter (spheroid)
converter =
  ConformalLatitudeConverter with properties:

    Eccentricity: 0.0818
>> chi = [0:15:90]
chi =
     0    15    30    45    60    75    90
>> inverse (converter, chi)
ans =
         0   15.0965   30.1669   45.1923   60.1662   75.0958   90.0000
>> forward (converter, chi)
ans =
         0   14.9041   29.8337   44.8077   59.8332   74.9036   90.0000

In Octave, with con2geo:

>> con2geo (chi)
ans =
         0   15.1117   30.1821   45.1923   60.1510   75.0806   90.0000
>> geo2con (chi)
ans =
         0   14.9041   29.8337   44.8077   59.8332   74.9036   90.0000


so it seems geo2con.m works properly but con2geo.m's output doesn't match Matlab.


Philip Nienhuis <philipnienhuis>
Group Member
Sat 07 Aug 2021 07:11:07 PM UTC, comment #8: 

Pushed geo2auth.m, geo2con.m, geo2iso.m and geo2rect.m here:
http://hg.code.sf.net/p/octave/mapping/rev/61a916b69206

Some of the functions had slight edits.
I've simply entered a for loop in meridianarc.m to allow vector/array input values and so could "vectorize" geo2rect.m.
Not beautiful but It Just Works.

On to the next contributions from comment #6 ... (may take a little time.)

It should be possible to "classdef-ify" the lot for max. Matlab compatibility but let's ignore that for now.

Thanks very much so far :-)

Philip Nienhuis <philipnienhuis>
Group Member
Sat 07 Aug 2021 12:56:36 PM UTC, comment #7: 

Thanks, I've fixed sph_chk and added some Texinfo.

I'm busy with the actual fuctions.

Philip Nienhuis <philipnienhuis>
Group Member
Sat 31 Jul 2021 05:13:12 PM UTC, comment #6: 

Added auth2geo, con2geo rect2geo

Also found
sph_check is missing 2 closing
L22 needs a }
L29 needs a )

Still working on iso2geo

(file #51717, file #51718, file #51719)

Anonymous
Tue 27 Jul 2021 08:50:35 PM UTC, comment #5: 

Aren't there ready-baked formulas Out There? Someone / several people must have thought about it, not just The Mathworks.

I suspect that the Boost library and/or the Proj library have these functions. But then we'd have to create an .oct file. Not the hardest thing to do, I've worked at those for the rasterread functions and for cell2sv/csv2cell in the io package. But that needs a bit more more work.

Philip Nienhuis <philipnienhuis>
Group Member
Tue 27 Jul 2021 06:46:36 PM UTC, comment #4: 

There is an inverse function.
https://www.mathworks.com/help/map/ref/map.geodesy.authaliclatitudeconverter.inverse.html

But this was not on the Octave wiki.

Using Newton-Raphson or something similar.
https://en.wikipedia.org/wiki/Latitude#Inverse_formulae_and_series

Will look into this.

Anonymous
Sun 25 Jul 2021 08:44:52 PM UTC, comment #3: 

Do we also need the inverse operations? (con2geo, iso2geo, rect2geo, auth2geo)

Philip Nienhuis <philipnienhuis>
Group Member
Fri 16 Jul 2021 09:30:29 PM UTC, comment #2: 

Well, vectorizing meridianarc is hard as quad() (a built-in Fortran-based function) isn't reentrant.
I tried with quadv (in Octave's scripting language) but that also doesn't accept vectors as it stands; that is, intensive surgery is required there. For a core function that seems a bridge to far for me.

It seems we cannot avoid a for loop in meridianarc.

Philip Nienhuis <philipnienhuis>
Group Member
Wed 17 Feb 2021 08:52:47 PM UTC, comment #1: 

I'll look into it.
Vectorizing should be easy but can be a bit of an art. Sometimes one needs a little luck to see a way out.

S/th for next weekend, if he weather is lousy enough :-)

Philip Nienhuis <philipnienhuis>
Group Member
Wed 17 Feb 2021 06:56:44 AM UTC, original submission:  

These should the functions behind, the classdef functions:

map.geodesy.AuthalicLatitudeConverter etc.

geo2rect cannot accept vectors.
meridianarc needs to be updated to accept vectors first.
How is this done with special functions?

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #52781:  geo2iso.diff added by None (2KiB - application/octet-stream - version 2.)
file #52782:  iso2geo.m added by None (3KiB - application/octet-stream - version 2.)
file #52141:  con2geo.m added by None (944B - application/octet-stream)
file #52142:  iso2geo.m added by None (3KiB - application/octet-stream)
file #51717:  auth2geo.m added by None (3KiB - application/octet-stream)
file #51718:  con2geo.m added by None (3KiB - application/octet-stream)
file #51719:  rect2geo.m added by None (3KiB - application/octet-stream)
file #50867:  geo2auth.m added by None (3KiB - application/octet-stream)
file #50868:  geo2con.m added by None (3KiB - application/octet-stream)
file #50869:  geo2iso.m added by None (3KiB - application/octet-stream)
file #50870:  geo2rect.m added by None (3KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by philipnienhuis
  • -email is unavailable- added by siko1056 (Updated the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only logged-in users can vote.

     

    Follow 17 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-02-01 philipnienhuis StatusIn Progress Done
        Open/ClosedOpen Closed
    2022-01-31 None Attached File- Added geo2iso.diff, #52781
        Attached File- Added iso2geo.m, #52782
    2021-10-21 None Attached File- Added con2geo.m, #52141
        Attached File- Added iso2geo.m, #52142
    2021-08-07 philipnienhuis StatusNone In Progress
    2021-07-31 None Attached File- Added auth2geo.m, #51717
        Attached File- Added con2geo.m, #51718
        Attached File- Added rect2geo.m, #51719
    2021-07-20 philipnienhuis Carbon-Copy- Added philipnienhuis
    2021-07-16 philipnienhuis Assigned toNone philipnienhuis
    2021-04-28 siko1056 Summary[octave forge](mapping) Auxillary Latitudes [octave forge] (mapping) Auxillary Latitudes
    2021-02-17 None Attached File- Added geo2auth.m, #50867
        Attached File- Added geo2con.m, #50868
        Attached File- Added geo2iso.m, #50869
        Attached File- Added geo2rect.m, #50870

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code