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

 
 

patch #9974: [octave forge](mapping) gc2sc scxsc gcxsc

Submitted by:  None
Submitted on:  Sat 19 Sep 2020 04:06:34 AM 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)

Mon 12 Oct 2020 03:37:34 PM UTC, comment #17: 

There is still rhxrh, but I do not see much priority for it. (But I can be convinced.)

Anyway, I'm grateful for the many contributions you've made to date.
I don't mind having to do a bit of work, after all this is hobby; the issue rather is about how to fit it into limited available time. Luckily (or not) my company's network has been down since Wednesday last week so this time I could allocate some time.
Circle intersections was an instructive subject for me. It reminded me of my geography project in high school (or is it called college?) about navigation on the globe, one of my other hobbies is sea-sailing; both relate to that same subject. Nice.

If you want to switch to and do work on some particular subject it may be wise to contact me first, to avoid doing stuff I've already done some work on or I know someone else has worked on.
Privately by email (my address is about everywhere), or if you feel uneasy there (you're still anonymous to me) I can also open a task on the task tracker to discuss next steps. But I'll have to find out first if anonymous comments can be made there.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Mon 12 Oct 2020 01:38:49 PM UTC, comment #16: 

Sorry you had to do a lot of work.

In the wiki grouped list is the Coordinates, Geodesy, and Projections section. This is the area I have been working on.
Modelling the Earth subsection is almost complete. I will need to learn more for the other subsections.

Alexander Barth was the last author for both distance and reckon.
Maybe the author could provide some insight.

Anonymous
Sun 11 Oct 2020 08:54:54 PM UTC, comment #15: 

Pushed scxsc.m here:
http://hg.code.sf.net/p/octave/mapping/rev/d878d7828ae1

Pushed gcxsc.m here:
http://hg.code.sf.net/p/octave/mapping/rev/1d59b083067f

... and I fixed a few other things in gcxgc.m.

All in all this patch was much, much more work than I expected :-)
While testing all these 3 functions + gcxgc.m together, ever more corner cases turned up. Let's say I've earned a lot.

You were right about the haversine formula being more accurate than the spherical cosine formula for small values. After some experimenting I found that below 25-40 degrees the haversine formula is more accurate and even much more accurate once the spherical distance gets below 5 degrees. So the Wikipedia article about the haversine formula (comment #10) could use some improvement.
For spherical distances over 179 degrees the haversine formula degenerates quickly. In scxsc.m I had to make a split to get suffcient accuracy for small and large distances between circle centres.
All in all scxsc.m is now very accurate AFAICS.

The function distance.m that also invokes the spherical cosine formula is similarly affected for small point distances. I looked into it but found it has been written in such a convoluted way that fixing it is unduly hard for now. Maybe later.

The inaccuracy of gc2sc.m I reported in comment #14 is caused by reckon.m being slightly inaccurate. I suppose there are also formulas for reckoning based on haversine and coversine functions but I couldn't find one, nor a useful reference. Can you?

Anyway, closing report.

I'll play around with the new functions a little the coming days, if I find no more bugs I'll make a new mapping package release.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Sat 10 Oct 2020 11:46:41 AM UTC, comment #14: 

Well, maybe not so fast ;-)

While fixing the scxsc.m tests I now hit the issue that gc2sc.m is a bit inaccurate.
As touched in passing in comment #11, in test #2 the great circles from test #4 in gcxgc.m are morphed into small circle notation, but along the way enough precision is lost to screw up the test in scxsc.m.
Maybe I could use your assistance in improving gc2sc.m's precision.

About precision and haversines, indeed for small radii scxsc.m gets inaccurate, but trying with haversines only gave me more inaccuracy. I want to figure that out more thoroughly.
You can try for yourself with tangent circles, e.g.:

[la, lo] = scxsc (0, 0, 1, 0, 2, 1)

and/or maybe smaller radii and thus smaller distances between circle centers. Or even a degenerate case with zero radius:

[la, lo] = scxsc (0, 0, 0, 0, 2, 2)

For large circles (r > 30 degrees) scxsc.m turns out to be quite accurate.

In the wiki there's a long list of functions needing implementation.

Would you have any experience and/or interest in projections?
I'd very much like to have projection functions in the mapping package based on the proj package. I've already tried a few things but very basic. Once projections are implemented a next step would be implementing geoshow etc.
One difficult step is to parse projection strings into something more manageable.

So you see there s always work to be done ;-)

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Sat 10 Oct 2020 10:32:02 AM UTC, comment #13: 

Perfect.

I will look for other functions to work on.
There are not many left for me.

Anonymous
Thu 08 Oct 2020 08:26:11 PM UTC, comment #12: 

Yeah I know the Wikipedia article you refer to, it explicitly says:

On computer systems with low floating-point precision ...

(emphasis by me)
where I guess they refer to 32-bit floats ("singles").
On 64-bit HW (almost all computers these days) it is of no more relevance.
But you are right in that it's always worth a check :-)

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Thu 08 Oct 2020 08:19:46 PM UTC, comment #11: 

As my company's network was down all day (and largely unusable yesterday) I made some progress. In fact I have scxsc.m almost finished.

Indeed my assumption that once all corner cases (actually implicit assumptions in the math) are out of the way, the math & strategy will work w/o any checks, turned out to be a lucky guess.

I revamped the sin/cos into sind/cosd (easier debugged when checking values in the variable editor :-) ). It now yields the exact same outcomes as Matlab.

Please ignore comment #8, that was written at a time I got lost myself.
BTW, in comment #9 where I write "_all tests in your code_" I didn't mean the BIST tests! rather the if blocks checking for e.g., identical axes etc.

Currenty I'm trying to fix the BIST tests. Using the "new" scxsc.m I discovered that the 4th test in gcxgc.m which got inherited as test #2 in scxsc, is in error - only the 2nd great circles there coincide. So gcxgc.m needs fixing as well (and scxsc.m proves to work very well).

Hopefully I can get all this finished soon - I'd like to make a new mapping pkg (+ io package) release next week.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Thu 08 Oct 2020 05:38:58 PM UTC, comment #10: 

Thank you for looking into this.

I have read and understand how I could have used only vectors.
I hope to have some time to look and understand comment #8 and comment #9

It looks like you are seeing some good corner cases to test.
As you said s.o. was a good starting point.
If you get to it before me there would be no issues.

https://en.wikipedia.org/wiki/Great-circle_distance#Computational_formulas

States specifically use haversine for small distances.
If your results are robust then the code can be cleaner.

Anonymous
Wed 07 Oct 2020 10:18:21 PM UTC, comment #9: 

I had a little time to look into it further. Next week I might do more.
But I think I have nailed the issues.

My ideas (to be validated !) are as follows:

The strategy & math in the S.O. article mentioned in scxsc.m is absolutely valid, but fragile. I think before it can be safely invoked various cases have to be catched to make sure the implicit assumptions are valid.

Case 1 is non-intersecting small circles because they're too far apart. Easily checked, just see if the smallest spherical distance between the circle centers is greater than the sum of their radii.

Case 2 is coinciding axes. Easily checked by comparing the circle centers or antipodes to differ less than -say- 2*eps. Once there:
Case 2A is different circle radii ==> no intersections.
Case 2B is equal radii ==> coinciding circles, infinite intersections.

All three lead to the same outcome, but different warnings.

All these cases can be explored by first computing the spherical distance between the circle centers for all pairs.
In my spherical trig goniometry textbook (from 1940, a classic) I found the spherical cosine formula that after experimenting seems to be sufficiently robust and accurate (no haversines etc. needed with 64-bit floats).
Template / example function attached.

After all these cases have been neutralized, I think all tests in your code can be dropped as there are no more mathematical issues to watch out for, not even sqrt of small negative values leading to atan2 choking on complex numbers.

Will you have a go at it or shall I do it next week or the week after?

(file #49935)

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Fri 02 Oct 2020 06:29:17 AM UTC, comment #8: 

Just to keep you informed:
Yesterday I've been playing a bit with the logic in the Stackoverflow post mentioned in scxsc.m. Using sind/cosd/atan2d I got it to work perfectly for small circles on perpendicular axes (lat1/lon1 = (0, 0); lat2/lon2 = (0, 90) with radii of 45 degrees, i.e. output of 2 X lat/lon = (0, 45), the one intersection/tangential point.
But any other values (e.g., radii of 45 and 46 degrees) gave wrong results (like 134.some for lon and/or lat => wrong hemisphere). Could be an atan2 input issue (IOW output is not the principal value). But also smaller values failed in other ways.
The SO procedure seems mathematically correct to me but I'm a bit wary of the corner cases.

In any case there should be checks on x02 and n2 to be max. 1, to avoid complex intermediate outcome when later on sqrt() is fed a very small negative value.

Which makes me think that the scxsc.m code at large may need to become more robust against internal rounding errors / inaccuracies of the goniometric functions. Also, as stated in the S.O. post, the cross product's accuracy may be a factor.

To be continued :-)

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Wed 30 Sep 2020 09:43:32 PM UTC, comment #7: 

Attached is my version so far, based on yours in comment #6.
I had a bit of trouble vectorizing the two stanzas, but I finally realized that I tried it with three pair of circles all the time which led cross() to take the wrong dimension. cross(), like dot() needs unambiguous indication of which dimensions to work on. From then on it was not too hard.
I've also added a third output indicating the intersection status of each pair of circles.

Comparing with Matlab, I find some consistent differences.

A perhaps small issue:
Experimenting to get two small circles with just one intersection point gives for Octave:

>> [u, v, w] = scxsc (0, 0, 1, 0, 2, [1.0000000002340808, 1.0000000002340809]
## (where the first radius lead to disjoint circles, and the second leads to one intersection point)

but Matlab needs:

>> [u, v] = scxsc ([0 0],[0 0], [1 1], [0 0], [2 2], [0.999998 .999999])
## ditto

Maybe Matlab seems to work with the "flat" circle radius (of the circle in the plane that intersects the surface) where the Octave version seems to take the radius over the surface of the sphere? Anyway it's just small differences.

A more profound issue is this:
Consider two circles with centers 90 degrees apart and then check with radii of 44 and 46 degrees. Only the second should have intersection points.
Octave does:

>> [u, v] = scxsc ([0 0], [0 0], [45 45], [0 0], [90 90], [44  46])
:
u =
   NaN   NaN
   NaN   NaN

v =
   NaN   NaN
   NaN   NaN

and complains about coinciding circles (which they aren't, their axes are perpendicular).
Matlab does (correctly IMO):

>> [u, v] = scxsc ([0 0], [0 0], [45 45], [0 0], [90 90], [44  46])
:
u =
        NaN       NaN
     7.5908   -7.5908
v =
    44.4912   44.4912

Your first and second scxsc versions and mine yield identical results, so I suppose there's something wrong in the applied math somewhere.

(file #49897)

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Sun 27 Sep 2020 03:15:23 PM UTC, comment #6: 

Ok now it returns vectors

[g, h] = scxsc (a, b, c, d, e, f)
warning: Octave: coinciding-small-circles
          scxsc: no solution circle centers are either the same or antipodal.
warning: Octave: coinciding-small-circles
          scxsc: no solution circle centers are either the same or antipodal.
g =

  -57.997936587566194    1.012254952815915
                  NaN                  NaN
   57.997936587566194   -1.012254952815915
   50.768479516407744   -0.886077123792614
   33.537497989074723   -0.585339762790331
                  NaN                  NaN
  -33.537497989074758    0.585339762790332
  -50.768479516407737    0.886077123792614
  -57.997936587566194    1.012254952815915

h =

  -102.764389682754597     1.348015699823234
                   NaN                   NaN
   122.235610317245346    -1.008178790369112
   144.735610317245374    -0.615479708670388
   167.235610317245374    -0.222780626971663
                   NaN                   NaN
  -147.764389682754626     0.562617536425785
  -125.264389682754654     0.955316618124509
  -102.764389682754597     1.348015699823234

There are some places I did not know how to do with out looping.
I marked these with TODO.

Added this as a test

test scxsc
PASSES 20 out of 20 tests

(file #49868)

Anonymous
Thu 24 Sep 2020 11:18:28 AM UTC, comment #5: 

Thank you for looking into this.

From comment #2 issues

Putting the code block at L111

>> [g, h] = scxsc (a, b, c, d, e, f)
error: horizontal dimensions mismatch (9x3 vs 1x1)
error: called from
    scxsc at line 77 column 10

Which is different from your error.

I then followed a gcxgc and placed the code block after the isreal section L77.

 [g, h] = scxsc (a, b, c, d, e, f)
error: Octave: coinciding-small-circles
          scxsc: no solution circle centers are either the same or antipodal.

However following your idx example and doing one at a time values are returned for all 6.

Will look into these issues, then look into comment #3
I think it is a good idea

Anonymous
Wed 23 Sep 2020 05:20:38 PM UTC, comment #4: 
Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Wed 23 Sep 2020 01:50:43 PM UTC, comment #3: 

Maybe we could return an extra output arg, similar to gcxgc.m, indicating conditions of the intersection points, a la:

  • 0: intersecting circles with two distinct intersection points
  • 1: coinciding circles (NaNs returned)
  • 2: tangent circles, just one intersection point returned as duplicates
  • -1 or 3: small circles too far apart, no intersections, NaNs returned.

What do you think?
If it's too difficult I could give it a try once you've solved the other issues in scxsc.m mentioned in comment #2.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Tue 22 Sep 2020 07:28:39 PM UTC, comment #2: 

Playing around with gc2sc.m and scxsc.m yields a few issues:

  • gc2sc with vector arguments always returns a scalar for radius.

Cause:
  range = pi / 2;
  radius = range;
...should read:
  range = pi / 2;
  radius = repmat (range, size (vect(:, 1)));
and a little lower:
    radius = rad2deg (range);
...should read:
    radius = rad2deg (radius);
and then it's fixed.

I'll fix gc2sc.m and polish some style things and push it soon.

Then,

  • scxsc.m doesn't really accept vectors.

Copying over L.94-112 from gcxgc.m before L.111 fixes that but then still the code doesn't work properly.
What I did:

## These are copied from test #3 in gcxgc:
[a, b, c] = gc2sc (45, [0:45:360], 45)
[d, e, f] = gc2sc (-45, -135, -45)
## and then just:
>> [g, h] = scxsc (a, b, c, d, e, f)
g =
   NaN   NaN
h =
   NaN   NaN

However, feeding entries from a, b, and c one by one yields correct results when compared with the results of test #3 in gcxgc run as:
[lats, lons, idl] = gcxgc (45, [0:45:360], 45, -45, -135, -45)
e.g.,:

>> id = 3; [g, h] = scxsc (a(id), b(id), c(id), d,  e, f)
g =
   57.998  -57.998
h =
   122.236   -57.764

  • I didn't see any scxsc warnings for coinciding circles (id = 2 and 6), NaNs should be returned there.

Could you please fix the above issue in scxsc?

In the mean time I'll next look at gcxsc.m.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Sun 20 Sep 2020 06:31:44 PM UTC, comment #1: 

Thanks for the contribution (that you announced in patch #9958).

I'll have a look in the coming weeks.

Philip Nienhuis <philipnienhuis>
Project MemberIn charge of this item.
Sat 19 Sep 2020 04:06:34 AM UTC, original submission:  

Modified Alfredo Foltran's gc2sc.
Then found a walk through example for scxsc.

gcxsc is wrapper of the two.

Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #49935:  sphdist.m added by philipnienhuis (188B - application/octet-stream)
file #49897:  scxsc.m added by philipnienhuis (9KiB - application/octet-stream)
file #49868:  scxsc.m added by None (7KiB - text/x-objcsrc - version 2)
file #49816:  gc2sc.m added by None (5KiB - text/x-objcsrc)
file #49817:  scxsc.m added by None (6KiB - text/x-objcsrc)
file #49818:  gcxsc.m added by None (4KiB - 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 10 latest changes.

Date Changed by Updated Field Previous Value => Replaced by
2020-10-11 philipnienhuis Open/ClosedOpen => Closed
2020-10-11 philipnienhuis StatusIn Progress => Done
2020-10-07 philipnienhuis Attached File- => Added sphdist.m, #49935
2020-09-30 philipnienhuis Attached File- => Added scxsc.m, #49897
2020-09-27 None Attached File- => Added scxsc.m, #49868
2020-09-23 philipnienhuis StatusNone => In Progress
2020-09-20 philipnienhuis Assigned toNone => philipnienhuis
2020-09-19 None Attached File- => Added gc2sc.m, #49816
    Attached File- => Added scxsc.m, #49817
    Attached File- => Added gcxsc.m, #49818

Back to the top


Powered by Savane 3.5