patchGNU Octave - Patches: patch #9958, [octave forge](mapping) gcxgc

 
 

patch #9958: [octave forge](mapping) gcxgc

Submitted by:  None
Submitted on:  Mon 20 Jul 2020 01:55:03 PM UTC  
 
Category:  Forge : new function Priority:  5 - Normal
Status:  Done Privacy:  Public
Assigned to:  Philip Nienhuis <philipnienhuis> Originator Email:  -email is unavailable-
Open/Closed:  Closed

Add a New Comment (Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

( Jump to the original submission)

Fri 21 Aug 2020 09:43:28 PM UTC, comment #30: 

I just realized the great circle pair coincidence tests could be even much simpler by just comparing longitudes of intersections with the equator, after first deselecting & if needed processing those which already coincide with the equator.
Modification + extra test pushed here;
http://hg.code.sf.net/p/octave/mapping/rev/1c9031e0f161

and adapted gcxgc.m attached.

(file #49694)

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Mon 17 Aug 2020 05:54:29 PM UTC, comment #29: 

Thanks for the heads up, fixed in repo.
I decreased the tolerance as this gcxgc.m seems fairly accurate.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Mon 17 Aug 2020 01:26:41 PM UTC, comment #28: 

The asserts need to be fixed

+%!test

+%! [lat3, lon3] = gcxgc ( 51.8853, 0.2545, 108.55, 49.0034, 2.5735, 32.44);

+%! assert (degrees2dms (lat3(1)), [50 54 27], 10-3)

+%! assert (degrees2dms (lon3(1)), [04 30 31], 10-3)

So it needs to be 10e-3
See comment #24

Anonymous
Sun 16 Aug 2020 08:24:23 PM UTC, comment #27: 

Just to be clear, your latest version of comment #15 was also more accurate then Matlab. Especially when computing intersection points at the equator.
I'm still intrigued why it yields false positives for coinciding great circles. Maybe some unfortunate sin/cos or +- swap somewhere in the formulas?

Anyway after some more adaptations I've pushed gcxgc3.m fro comment #22, renamed to gcxgc.m, here:
http://hg.code.sf.net/p/octave/mapping/rev/4bad0a437b90
I still perceive it as a first go, hopefully the version from comment #1 can be fixed one day; but for now I won't spend more time on it.
Things I adapted:

  • Option for a single output (cf. Matlab)
  • Dropped the requirement of Longitudes and Azimuths to be in the range [-180,180] 0r [pi,pi], as Matlab happily accepts those and they are implicitly wrapped into that interval in our version.

Matlab still accepts azimuths or 90 degrees and yields finite answers; a bug if you ask me. But maybe that has been fixed in later releases; I have only r2014a and didn't have time to install the prerelease 2020b yet.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Fri 14 Aug 2020 02:26:21 PM UTC, comment #26: 

What's your opinion?

This is no problem for me, accurate is important.

You also have the ability to compare against matlab.

Anonymous
Fri 14 Aug 2020 12:59:41 PM UTC, comment #25: 

I'm reviewing & adapting your latest contribution in more detail.
Some code needs to be properly "vectorized" (actually: adapted to allow vector inputs), a.o. L.182 & 185 (".^" rather than just "^") and the cross product in L.134 that needs a 3rd DIM argument (in this case: 2).
While similarly adapting the if block in L.154-166 it turned out that the invoked computation method is still unreliable:

## Correct:
>> [a, b] = gcxgc (45, 45, 90, -40, -135, 90)
warning: gcxgc: Non-unique intersection.
a =
   NaN   NaN
b =
   NaN   NaN

>> [a, b] = gcxgc (45, 45, 90, -40, -135, 90)
warning: gcxgc: Non-unique intersection.
a =
   NaN   NaN
b =
   NaN   NaN
## Incorrect. (45, 45) and (-40, -135) are no antipodes so
## even though the great circles have horizontal W tangents in
## those coordinates they do not coincide. Same of course for:
>> [a, b] = gcxgc (40, 45, 90, -40, -135, 90)
warning: gcxgc: Non-unique intersection.
a =
   NaN   NaN

b =
   NaN   NaN

The gcxgc3.m I uploaded yesterday works fine and gives more accurate results than Matlab. I tend to -for the time being- use that one until/unless you want to improve the code from the comment #15 version.

What's your opinion?

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Fri 14 Aug 2020 12:54:04 PM UTC, comment #24: 

Yes the error should be 10e-3 for everything

The first test changed to

%!test
%! [lat3, lon3] = gcxgc ( 51.8853, 0.2545, 108.55, 49.0034, 2.5735, 32.44);
%! assert (degrees2dms (lat3(1)), [50 54 27.387], 10e-3)
%! assert (degrees2dms (lon3(1)), [04 30 30.8687], 10e-3)

All others should pass with the 10e-3 change.

I will look at your other approach, to see where I can improve.

Anonymous
Wed 12 Aug 2020 09:06:21 PM UTC, comment #23: 

gcxgc3.m attached.

(file #49662)

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Wed 12 Aug 2020 09:05:18 PM UTC, comment #22: 

In the mean time I played around about detecting coinciding great circles with my idea from comment #18.
I got it working, see attached gcxgc3.m (based on your gcxgc.m version of comment #6). Just interpret it as a proof-of-concept for now.

Turned out it isn't even necessary to compute & then compare polar axes, it suffices to find & compare the (Lat,Lon) coordinates of the location where the great circle tangents run exactly EW (azimuth 90 degrees).
I wish I had realized this earlier on. Very easy to grasp.

In that gcxgc3.m I've also made it return the indices of coinciding great circle pairs, if any. A sort of Matlab extension so to speak.

I'll proceed with the version of comment #15 (your last version).

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Wed 12 Aug 2020 08:52:42 PM UTC, comment #21: 

One of those corrected test results is

>> test gcxgc
***** test
[lat3, lon3] = gcxgc ( 51.8853, 0.2545, 108.55, 49.0034, 2.5735, 32.44);
assert (degrees2dms (lat3(1)), [50 54 27], 10e-3);
assert (degrees2dms (lon3(1)), [04 30 31], 10e-3);
!!!!! test failed
ASSERT errors for:  assert (degrees2dms (lat3 (1)),[50, 54, 27],10e-3)

  Location  |  Observed  |  Expected  |  Reason
    (3)         27.387         27        Abs err 0.387 exceeds tol 0.01 by 0.4

where your gcxgc.m from comment #15 actually gives accurate results.

I'll fix those.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Wed 12 Aug 2020 08:34:45 PM UTC, comment #20: 

Not all tests pass well.
You've added a few tests like this:

%!test
%! [lat3, lon3] = gcxgc_1 (51.8853, 0.2545, 108.55, 49.0034, 2.5735, 32.44);
%! assert (degrees2dms (lat3(1)), [50 54 27], 10-3)
%! assert (degrees2dms (lon3(1)), [04 30 31], 10-3)

where I suppose you mean the tolerance to be 10e-3 (0.0001) but actually the tolerance is 10-3 = 7.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Mon 10 Aug 2020 12:00:49 PM UTC, comment #19: 

Step 5 obviously refers to identifying coinciding great circles.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Mon 10 Aug 2020 12:00:05 PM UTC, comment #18: 

It occurred to me that finding the polar axis for any great circle defined by a (Lat, Lon, Azimuth) trio is actually very easy, but involved.

  1. Find intersection with the equator (Lat 0, Lon 0, azi = 90).
  2. Get Lon values halfway between the intersections, and define great circle through those points and the N- an S-poles.
  3. Find intersection point between original great circle and the last great circle found in step 2; that's where the azimuth of the orig. great circle = exactly 90 degress, or E-W. Add 90 degrees to found Lat (wrap around 180 degrees or pi if needed & also wrap Lon around 360 or 2.pi), then the resulting (Lat, Lon) = coordinates of intersection of polar axis with sphere.
  4. Convert that to ECEF if needed.
  5. But, computing inner products may not even be needed, a simple comparison of (Lat, Lon) values + antipodes would do.

There are probably several optimizations and shortcuts possible.

If your original gcxgc is more accurate than the one based on converting to ECEF (XYZ), see comment #8, it might be worth pursuing this way; at least to know how much more accurate the answers would be; obviously there's a trade-off as regards overall performance. If you want/agree I'll give it a try when I find a spell.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Mon 10 Aug 2020 09:38:06 AM UTC, comment #17: 

"Copyright (C) 2020 The Octave Project Developers"
Would that suit you?

Perfect will add next time.

Anonymous
Fri 07 Aug 2020 08:54:59 PM UTC, comment #16: 

Thanks! I tried a few cases and the results seem fine now.

Later this weekend/ early next week I'll need to check what happens with input vectors of e.g., mixed coinciding and non-coinciding great circle pairs (L.154 and up). AFAICS that piece applies to scalar values only and won't work (properly) with vectors.
But if you want to beat me, be my guest :-)
Do we indicate index numbers of all coinciding pairs (could be a lot) or do we just mention that one or more pairs coincide? I'll check with Matlab.

I also have tests lying around for vector and mixed scalar/vector inputs and error checks for uneven length, I'll add those as well. I've waited with those until you had the algorithm sufficiently polished.

If you still want to remain anonymous, these days we have a new copyright string:
"Copyright (C) 2020 The Octave Project Developers"
Would that suit you? (pushing to the repo is still by me, can't be anonymous)

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Fri 07 Aug 2020 05:40:07 PM UTC, comment #15: 

Thank you for the tests

Added the warning.
Fixed for more NaN

test gcxgc
PASSES 20 out of 20 tests

 

(file #49642)

Anonymous
Fri 07 Aug 2020 11:53:00 AM UTC, comment #14: 

Thanks, we're getting closer I think.

Better morph the test for non-uniqueness from regular test into a warning check:

%!warning <non-unique> [lat3, lon3] = gcxgc (0, 0, 45, 0, 180, -45);

(note semicolons at end, had better be there for error checks as well although errors return no output values).

Another simple non-uniqueness test that still fails:

>> [a, b] = gcxgc (45, 45, 90, -45, -135, 90)
a =
   54.736  -54.736
b =
    45  -135

(note: horizontal EW azimuth at (45, 45) and its antipode (-45, -135)
where Matlab correctly returns NaNs and warns.

Same for about any coordinates & azimuth, their antipodes & (-rem (azimuth, 180))

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Thu 06 Aug 2020 10:23:03 AM UTC, comment #13: 

I added more BIST

Added \n for the warning

(file #49638)

Anonymous
Tue 04 Aug 2020 08:47:27 AM UTC, comment #12: 

Yeah you're absolutely right. What gcxgc and related functions need to accomplish looks easy at first sight but soon turns out to be very complicated, esp. getting sufficiently accurate answers.

As to preliminaries for great circle poles I found this site:
https://www.aa.quae.nl/en/reken/grootcirkel.html
where some apparently useful info is given.
It seems that many great circle problems for decades have been extensively investigated in aviation, but many of that literature is of the "gray" variety and hard to uncover.

I'm of little help here ATM, I simply lack time now although I find it a fascinating subject and very useful for the mapping package.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Tue 04 Aug 2020 07:24:01 AM UTC, comment #11: 

Do you also plan to supply gcxsc, scxsc and rhxrh ?

I need to understand this one first.
Also finding already presented solutions (not matlab's) is difficult.

Anonymous
Thu 30 Jul 2020 08:59:13 PM UTC, comment #10: 

That should read:
(the inner product of the polar axes vectors divided by the square of the radius of the sphere, i.e., the cosine of the angle between them)

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Thu 30 Jul 2020 08:54:16 PM UTC, comment #9: 

This condition at L.129:

:
  elseif (sind (rad2deg (vect(:, 1))) ==  sind (rad2deg (vect(:, 4))) &&
          sind (rad2deg (vect(:, 2))) == -sind (rad2deg (vect(:, 5))))
warning("gcxgc: Non-unique intersection.\n")

is too simple.
I think the only reliable way to detect coinciding great circles is to compare the (Lat, Lon) values of their poles at the surface of the sphere (i.e., the cross product of those polar axes).
Those can easily be deduced by the cross products of the input points and inferred points e.g., a quarter of the circumference away on the great circles. Or maybe directly from the given (Lat, Lon) points + azimuths, but I couldn't find an expression for the latter on the web. (And TBH my own spherical trig textbook is too thick :-) last time I seriously looked in it was 45 years ago at high school when I did a navigation project for geography.)

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Thu 30 Jul 2020 08:36:42 PM UTC, comment #8: 

Hmmm, this doesn't make sense:

>> [a, b] = gcxgc3 (0, 0, 90, 0, 180, 80)
warning: gcxgc: Non-unique intersection.
warning: called from
    gcxgc3 at line 132 column 5

a =
   NaN   NaN

b =
   NaN   NaN

The azimuths are clearly different, aren't they?
Matlab gives:

>> [a, b] = gcxgc3 (0, 0, 90, 0, 180, 80)
a =
   1.0e-6 *
   -0.3632   -0.3632
b =
    0.0000  180.0000
>>

and your previous version yields the same, even much more accurate than Matlab.

But then:

>> [a, b] = gcxgc2 (90, 0, 90, -90.000001, 180, 80) ## previous version
error: atan2: not defined for complex numbers
error: called from
    gcxgc2 at line 126 column 10
## new version
>> [a, b] = gcxgc3 (90, 0, 90, -90.000001, 180, 80)
a =
  -89.99999432871816   89.99999432871817

b =
  -90.00000003544301   89.99999996455699

==> I think there should be a check for latitudes at or beyond the N and S poles. That's a straightforward fix. Probably by chance your new version correctly detects latitudes exactly at the poles through the math, but not so beyond the poles.

(FYI Matlab accepts and gives finite answers for any latitude != 90 or -90 degrees, incl. +90.00001 and -90.00001 - the latter incorrectly IMO. Matlab doesn't document this or any wrapping around the [-90,90] interval)

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Thu 30 Jul 2020 06:29:03 PM UTC, comment #7: 

Could 3D vectors also be derived by e.g., geodetic2ecef using a spherical geoid? (just curious)
Matlab makes it a bit difficult by requiring one set of coordinates + azimuth values per circle; I was thinking maybe great circles input could also be specified using pairs of non-antipodal coordinates, as a sort of Octave-only extension. I have use cases for such occasions. But hmm, maybe later.

Anyway I think you did a fine job, nice!

I'll see if the check for coinciding great circles can be vectorized a little more and I'll ass a few more tests for vector inputs.
IMO the coincidental circles warning might need a "\n" at the end to suppress the backtrace (although Matlab supplies that as well).
I'll also add a stanza in the texinfo about coinciding circles.

Do you also plan to supply gcxsc, scxsc and rhxrh ?

FYI the functionality of polyxpoly is already present in the geometry package, only a wrapper would be required.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Thu 30 Jul 2020 03:15:38 PM UTC, comment #6: 

I like the improvements you did much cleaner

I changed to vectors instead of spherical trigonometry.

It passes the two new tests you provided.

When a test is given and a warning is expected, what should be done?

test gcxgc
warning: gcxgc: Non-unique intersection.
PASSES 17 out of 17 tests

(file #49584)

Anonymous
Tue 28 Jul 2020 04:05:41 PM UTC, comment #5: 

Yeah the answer is right, I don't understand the "ambiguity" the comment refers to.

As to style:

  • consistently parentheses around if conditions (e.g., L.91 and L.119 - in the latter you can just move the rightmost parenthesis behind the "0")
  • consistently spaces around assignment "=" and operators and 1 space between function name and opening parenthesis. So e.g., L.87:

del_12=2*asin(sqrt(sin(del_lat/2)^2+cos(vect(1))*cos(vect(4))*sin(del_lon/2)^2));
==>
del_12 = 2 * asin (sqrt (sin (del_lat / 2) ^ 2 + cos (vect(1)) * cos (vect(4)) * sin (del_lon / 2) ^ 2));
(You can add continuation markers ... and break statements over multiple lines)

  • double quotes " rather than ' around strings (L. 114).

Agreed, all nitpicking but this is the coding style adopted by Octave devs.

In the mean time I've already done this and also adapted the function to accept vector inputs. The intermediate result (gcxgc2.m) is attached.

(file #49560)

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Tue 28 Jul 2020 01:16:27 AM UTC, comment #4: 

Thank you for the response.

When I comment the L119

>> [lat,lon]=gcxgc (20, -5, 45, 30, 5, 15)
lat =

  -28.062   28.062

lon =

  -175.5879     4.4121

Will review code, the tests seen on the documentation gave correct results. Extra examples are good for additional BIST.

Agree with getting full compatibility.
What style issues need to be fixed?

Anonymous
Mon 27 Jul 2020 06:15:50 PM UTC, comment #3: 

Correction:

az = <any number in open interval (0 ... 180)>

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Mon 27 Jul 2020 04:18:09 PM UTC, comment #2: 

For a check on two identical great circles, try:

az = <any number in open interval (-90 ... 90)>
[lat, lon] = gcxgc (0, 0, az, 0, 180, -az)

Matlab r2014a gives:

>> format compact
>> [lat, lon] = gcxgc (0, 0, 45, 0, 180, -45)
warning: Non-unique intersection.
> In scxsc at 80
  In gcxgc at 63
lat =
   NaN   NaN
lon =
   NaN   NaN

The answers from your gcxgc seem correct but not unique. It would be nice if we could get full Matlab compatibility.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Sun 26 Jul 2020 06:38:38 PM UTC, comment #1: 

Thanks for this contribution. Same guy as patch #9492, patch #9634, ... ?

Function looks good, just minor style issues, but many BIST tests, good!

I checked with Matlab, and I want you to review the used formulas. With your function:

>> [a, b] = gcxgc (20, -5, 45, 30, 5, 15)
a =
  -28.062   28.062
b =
     0   180

while Matlab gives:

a =
   28.0620  -28.0620
b =
    4.4121 -175.5879

AFAICS the culprit is something in L.119 and up where a comment suggests "Causes ambiguity" and an if statement screws up the otherwise correct answer.

Also, I see no checks for identical great circle inputs (should return NaNs, cf. the on-line Matlab help text). As far as my own spherical trig proficiency goes (warning: rusty because long time ago!) this could be a simple comparison of the poles of the great circles. Perhaps / probably the answer is already in one of the intermediate results in your function.

Thanks

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Mon 20 Jul 2020 01:55:03 PM UTC, original submission:  

Finds intersection points of two great circles

Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #49694:  gcxgc.m added by philipnienhuis (9KiB - application/octet-stream)
file #49662:  gcxgc3.m added by philipnienhuis (7KiB - application/octet-stream)
file #49642:  gcxgc.m added by None (8KiB - text/x-objcsrc - version 4)
file #49638:  gcxgc.m added by None (7KiB - text/x-objcsrc - version 3)
file #49584:  gcxgc.m added by None (6KiB - text/x-objcsrc - version2)
file #49560:  gcxgc2.m added by philipnienhuis (6KiB - application/octet-stream)
file #49517:  gcxgc.m added by None (5KiB - text/x-objcsrc)

 

Depends on the following items: None found

Items that depend on this one: None found

 

CC list is empty

 

Do you think this task is very important?
If so, you can add your encouragement to it.
This task has 0 encouragements so far.

Only logged-in users can vote.

 

 

 

Follow 13 latest changes.

Date Changed by Updated Field Previous Value => Replaced by
2020-08-21 philipnienhuis Attached File- => Added gcxgc.m, #49694
2020-08-16 philipnienhuis StatusIn Progress => Done
    Open/ClosedOpen => Closed
2020-08-12 philipnienhuis Attached File- => Added gcxgc3.m, #49662
2020-08-07 None Attached File- => Added gcxgc.m, #49642
2020-08-06 None Attached File- => Added gcxgc.m, #49638
2020-07-30 philipnienhuis StatusNeed Info => In Progress
2020-07-30 None Attached File- => Added gcxgc.m, #49584
2020-07-28 philipnienhuis Attached File- => Added gcxgc2.m, #49560
2020-07-26 philipnienhuis CategoryNone => Forge : new function
    StatusNone => Need Info
    Assigned toNone => philipnienhuis
2020-07-20 None Attached File- => Added gcxgc.m, #49517

Back to the top


Powered by Savane 3.5