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

 
 

patch #10021: [octave forge] (mapping) meridianfwd

Submitter:  None
Submitted:  Sun 14 Feb 2021 09:45:00 PM 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 27 Jul 2021 08:45:46 PM UTC, comment #28: 

Maybe azimuth.m can be just a wrapper around geodeticarc.m?

I have no opinion on what's the best thing to do - just absorb (& polish) azimuth.m, or write a simple wrapper.

Philip Nienhuis <philipnienhuis>
Group Member
Tue 27 Jul 2021 08:42:05 PM UTC, comment #27: 

Yeah perhaps private/ is a better place for sph_chk.m to live in.
I've only pushed it yet - a next step is to invoke it in the various functions that use reference ellipsoids and -spheres.
Something for later this summer, I still don't have that much time for Octave.

Both vincenty.m and geodeticarc.m have as 2nd output parameter the azimuths at starting & end points; AFAICS that works quite well. I suppose that the azimuth of the starting point is what "azimuth.m" computes - can you confirm?
My thinking is that both functions compute the geodesics = shortest connections between two points, so the azimuths at the starting point should point in the optimal direction.
From azimuth.m at Alfredo Foltran's site:

:
## @example
## >> azimuth([10, 10], [10, 40])
## ans = 87.336
## >> azimuth([0, 10], [0, 40])
## ans = 90
## >> azimuth(pi / 4, 0, pi / 4, -pi / 2, "radians")
## ans = 5.3279
:


and current dev mapping package's geodeticarc.m:

>> [a, b] = geodeticarc ([10, 10], [10, 40])
a = 3.2880e+06
b =
   87.336   92.664

>> [a, b] = geodeticarc ([0, 10], [0, 40])
a = 3.3396e+06
b =
   90   90

>> rad2deg ([pi / 4, 0])
ans =
   45    0
>> rad2deg ([pi / 4, -pi / 2])
ans =
   45  -90
>> [~, b] = geodeticarc ([45, 0], [45, -90])
b =
   -54.726  -125.274
>> mean (b)
ans = -90
## makes sense, the entire displacement is due west = -90 degrees
>> deg2rad (b)
ans =
  -0.9552  -2.1864

>> 2 * pi + deg2rad (b)
5.3280   4.0967


so all basic functionality is there.
geodeticarc.m as it stands needs a few improvements:

  • allow radians as input
  • allow separate lat/lon values as inputs besides [lat lon] vectors/matrices
  • apply wrapTo360 / wrapTo2Pi to the azimuth outputs.


Well, something for this summer. Or would you have a go at it?

Philip Nienhuis <philipnienhuis>
Group Member
Tue 27 Jul 2021 07:03:36 PM UTC, comment #26: 

Perfect.

For my understanding should sph_chk be in inst/private?
This should be useful in having cleaner code.


Also A. Foltran's azimuth function allows a reference ellipsoid.
Should this be looked into?

Anonymous
Sun 25 Jul 2021 08:49:24 PM UTC, comment #25: 

A last bit of this patch report (sph_chk.m) has been pushed here:
http://hg.code.sf.net/p/octave/mapping/rev/c4af94e14e20

I've adapted it so that it can check for specific struct fields so that each individual function calling it just has to specify only the fields it requires. Proper class of the struct field values isn't performed though. Maybe later.

Philip Nienhuis <philipnienhuis>
Group Member
Fri 16 Jul 2021 09:00:07 PM UTC, comment #24: 

I've pushed A. Foltran's vincenty.m here:
http://hg.code.sf.net/p/octave/mapping/rev/e9c566d17498
and later applied some style fixes.

Next I wrote a new geodeticarc.m (ellipsoid "extension" of distance.m) and pushed that here:
http://hg.code.sf.net/p/octave/mapping/rev/578e79ca9dba

Closing report.

Philip Nienhuis <philipnienhuis>
Group Member
Mon 12 Jul 2021 08:21:02 PM UTC, comment #23: 

Sorted it all out - vincentyDirect and derived geodeticfwd are simply much more accurate than what is described in Vermeer's notes mentioned in comment #8.

I pushed some final changes to geodeticfwd.m here:
http://hg.code.sf.net/p/octave/mapping/rev/cd1b0bbc5f7f
Later on I'll check out if the accuracy can be improved for small values of the sin() functions in Vincenty's formulas, given the notorious accuracy degradation for small values (see also patch #9974).

Next I applied some style fixes to meridianfwd.m and pushed it:
http://hg.code.sf.net/p/octave/mapping/rev/ffb0b654559e

Formally this patch report is finished, hence Status "Ready for test", but there's some spin-off: I'll proceed with vincenty.m and geodeticarc.m as announced in comment #21.

Philip Nienhuis <philipnienhuis>
Group Member
Sun 11 Jul 2021 09:58:17 PM UTC, comment #22: 

Progress report:

After some (much) more experimenting with vincentyDirect.m, vincenty.m, and after trying various examples I got a bit confused.
The test examples in the original geodeticfwd.m (from the comment #0) can't be exactly reproduced by vincentyDirect.m as it stands.  The difference is about 10 km. But it looks like the original geodeticfwd.m was based on too simple calculation methods. Anyway this needs to be sorted out before further work on meridianfwd.m can continue.
That said, for angle displacement inputs vincentyDirect.m turns out to return quite exact answers AFAICS. Inputting vincentyDirect.m's output into vicenty.m and vice versa fairly exactly reproduces the input displacement and end coordinates, resp.

In the mean time I improved the "new" geodeticfwd.m:

  • absorbed the actual vincentyDirect.m code into it; geodeticfwd.m isn't a wrapper for vincentyDirect.m anymore but self-contained;
  • vectorized the code so multiple displacements can be computed relatively rapidly in one call;
  • geodeticfwd.m accepts scalars, vectors and 2D/ND arrays, provided all input vectors and arrays match regarding sizes and dimensions. Scalars on the one hand and vectors/arrays OTOH can be mixed;
  • also check the ellipsoid LengthUnit to match that of the "length" input (to avoid a potentially deceiving issue with length units).

These last changes are here:
http://hg.code.sf.net/p/octave/mapping/rev/347c2ac91480
There are still pending improvements: more BIST (esp. input validation), remove some debug statements; maybe ellipsoid validity checks. That'll happen later this week.

I did this to just leave A. Foltran's version as much as possible his work; I intend to keep vincentyDirect.m in the mapping package for at least some time.

Philip Nienhuis <philipnienhuis>
Group Member
Thu 08 Jul 2021 09:21:52 PM UTC, comment #21: 

I've pushed Alfredo Foltran's vincentyDirect.m function here:
http://hg.code.sf.net/p/octave/mapping/rev/e13cb68d8423
(fixed a minor typo along the way)

and then pushed a new geodeticfwd.m here:
http://hg.code.sf.net/p/octave/mapping/rev/7f9640d8956c
which is merely a wrapper around vincentyDirect.m that offers more flexible input.

I've added some BIST but will add more in the coming days/weeks.
meridianfwd is next to be added then.

Later on I'd also like to incorporate Alfedo's vincenty.m and a wrapper around it (probably called geodeticarc).

Philip Nienhuis <philipnienhuis>
Group Member
Sat 05 Jun 2021 08:37:40 PM UTC, comment #20: 

(coming back to this patch)
Looks like I'm able to spend time on this the next weeks.
I compared vincentyDirect.m's results to those of grodeticfwd.m but I see some different results; I'm just trying to find out where that comes from.
I'd first like to prepare one of those for including in the mapping pkg, after that I'll look into the reference ellipsoids checks.

Philip Nienhuis <philipnienhuis>
Group Member
Tue 06 Apr 2021 09:48:25 PM UTC, comment #19: 

If a user supplies
The two fields
SemimajorAxis
InverseFlattening

It would work similar to the vector version and calculate the other fields.
The ones with strings would be

    Code = ""
    Name = ""
    LengthUnit = meter (default)


I would think most users would do

E = wgs84Ellipsoid;
function (foo,bar,E,"degrees")


Anonymous
Tue 06 Apr 2021 07:39:21 PM UTC, comment #18: 


> Check all the fields: No, I would like a better solution.
>If a user gave 2 things semi major axis and inverse flatting  all the maths would work, and be unitless. The rest are just strings.

An ellipsoid struct by definition has several fields, so some checks are needed but not for all fields. Maybe just the "defining properties" and the rest can be corrected/augmented by paramcalc() (with a warning in that case).
Supplying a vector with just SMA and (inverse) flattening is the other option in __ell_chk__.m , so that's covered.

I have to dig up the exact reckon issues I hit with patch #9974 comment #14 and comment #15 - the very same issues apply here but in somewhat different disguise.
I d/led 'vincentyDirect.m' from Alfredo Foltran's site, AFAICS it does the same as geodeticfwd.m but maybe more solidly. Eve then it needs tests and vectorization however.

Stay tuned.

Philip Nienhuis <philipnienhuis>
Group Member
Tue 06 Apr 2021 03:18:18 PM UTC, comment #17: 

Rename: no issue.
I thought if it was in the private folder was all that is needed.

Error message handling: Are there other packages which have this?

Check all the fields: No, I would like a better solution.
If a user gave 2 things semi major axis and inverse flatting  all the maths would work, and be unitless. The rest are just strings.
So making this a private function is fine.

Focus efforts on the patch instead.
Not trying to create distractions.


Anonymous
Mon 05 Apr 2021 02:17:18 PM UTC, comment #16: 

Sorry for a bit of delay, other patches + real life crept in.

IMO sph_chk.m had better be called "ell_chk.m" as it checks ellipsoids rather than spheroids (they seem to have different fields). Moreover as it isn't meant to be a user-callable function  "__ell_chk__.m' is even better.
I've attached a first draft here, actually a renamed and morphed sph_chk.m

TBD: find a way to get the name of the calling function (for the error messages). TTBOMK there's no way for that in .m-file functions, so maybe the caller needs to specify its name as 2nd arg.

I'm wondering if we need to rigidly check all fields. Cf. Matlab, ellipsoids have defining properties "Code", "Name", "SemimajorAxis", "SemimiorAxis", "LengthUnit", "InverseFlattening" and "Eccentricity", while "Flattening, ThirdFlattening", "MeanRadius", "SurfaceArea" and "Volume" are "additional" (derived?) properties. The latter can be derived / computed using paramcalc (now a subfunc in referenceEllipsoid.m). We can separate paramcalc and make it a private function.

What is your opinion on all these points?

The next days/weeks I'll try to work on geodeticfwd and meridianfwd, sph_chk/__ell_chk__ is actually an auxilliary thing (maybe even a distraction :-) ).

(file #51193)

Philip Nienhuis <philipnienhuis>
Group Member
Fri 19 Mar 2021 09:05:08 AM UTC, comment #15: 

Improved the checks for the structure.

Your last statement should be caught by the referenceEllipsoid

(file #51091)

Anonymous
Fri 12 Mar 2021 11:17:44 AM UTC, comment #14: 

Sorry I have little time to spend on it.
A quick scan w.r.t. sph_chk.m:

  • I like the idea vert much, could help reduce a lot of code duplication in the referenceEllipsoid/-Spheroid/whatever functions.
  • The 1st elseif clause (isstruct...): not just check struct but also if it is a valid spheroid/ellipsoidstruct as well, by checking field names and class of field contents. My thought is that anything which passes sph_chk.m will be fed blindly to the functions calling it and those must rely on valid contents of output 'E'; otherwise users will get errors they might be able not relate to their input.
  • Shouldn't the last 'else' be an elseif checking for char input (name of ellipsoid)? Then there should be another 'else' withan error message


Philip Nienhuis <philipnienhuis>
Group Member
Wed 10 Mar 2021 11:15:50 PM UTC, comment #13: 

Yes, will work on the meridianfwd wrapper after.
What are your thoughts about the private function in comment #2?

Anonymous
Sat 06 Mar 2021 10:33:10 PM UTC, comment #12: 

Alfredo Foltran confirmed so I'll try to his vicentyDirect.m function to see how much better it works (accuracy, performance) and if OK, vectorize his function, get it into the mapping pkg.

If it all works out OK, would you see how it can be invoked behind meridianfwd instead of geodeticfwd?

BTW Somewhere I read some critical comments on the Vicenty formulas but first let's see how it behaves.

Alfredo Foltran also has another reckon.m version. TBH I'd first need to rediscover where the reckon.m currently in the mapping package is inaccurate.

Philip Nienhuis <philipnienhuis>
Group Member
Fri 05 Mar 2021 12:53:15 PM UTC, comment #11: 


> Should this patch be closed?

Not yet I think. You supplied meridianfwd.m which is useful in itself.

I'll ask Alfredo Foltran if he agrees with absorbing some of his functions in the mapping package. It's all GPL so no strict need but it's the minimal politeness for these sort of cases.
(ISTR his name from some other occasion)
Also vicentyDirect.m needs some work.

I'll keep you informed here

Philip Nienhuis <philipnienhuis>
Group Member
Fri 05 Mar 2021 12:25:51 PM UTC, comment #10: 

Agreed no sense in redoing work.
Since you state this appears to work with the current cases.

Should this patch be closed?


Anonymous
Thu 04 Mar 2021 10:41:09 PM UTC, comment #9: 

Hmmm, I wonder if Alfredo Foltran has already sorted it out, why bother to reinvent the wheel?

His code needs to be adapted to accept vectors (or even matrices) but otherwise it works good on a Unit Sphere.

Philip Nienhuis <philipnienhuis>
Group Member
Thu 04 Mar 2021 06:12:16 PM UTC, comment #8: 

Look at Vermeer's notes pg 8.


A good start for Vincenty's method would be here:

https://github.com/alfoltranteam/octave-map/blob/master/src/vincentyDirect.m




Anonymous
Thu 04 Mar 2021 04:53:42 PM UTC, comment #7: 

I read (IIRC on wikipedia) that the formulas generally used for this operation (due to Vincenty) tend to break down on antipodal points [*].
Yet reckon.m seems to do fairly well for that case. Surprising as I found that is is not very accurate in some other cases (not quite corner cases).
Where did you get the used equations from? couldn't find a reference in the code.

I already had inserted 'wrapToPi (az)' as L.142 (see comment #4) but it didn't make much difference here, still imaginary longitudes.

Indeed, I figured the erroneous longitudes would be due to a quadrant / principal value issue.

I really hope the issues can be solved. geodeticfwd.m seems a very useful function to me.

[*] https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid

Philip Nienhuis <philipnienhuis>
Group Member
Thu 04 Mar 2021 01:04:13 PM UTC, comment #6: 

Thank you for the test cases

Adding on L142

az = wrapToPi (az);


Takes care of the imaginary portions.

The other is a quadrant issue
Doing an asin (0) but the domain is (-90, 90)

reckon is using atan2 to avoid this issue.

So as in comment #4
There needs be more for L143-144

Anonymous
Sat 27 Feb 2021 11:32:33 AM UTC, comment #5: 

Also,

>> [lat, lon] = geodeticfwd (0, 0, [0 90], 2 * pi() / 2, "Unit Sphere", "degrees")
lat =
   7.0165e-15   4.2962e-31

lon =
            0   7.0165e-15


where I'd expect lon to be [180 180] (going from equator to antipode on equator, over the North Pole and along the equator).
In fact any azimuth value should lead to this result.
E.g.,

>> [latr, lonr] = reckon (0, 0, 180, [-90 : 10: 90])
latr =
   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0

lonr =
 Columns 1 through 14:
  -180  -180  -180  -180  -180  -180  -180  -180  -180  -180  -180  -180  -180  -180
 Columns 15 through 19:
  -180  -180  -180  -180  -180


Moreover,

>> [lat, lon] = geodeticfwd (0, 0, [-90 : 30 : 90], pi()/2, "Unit Sphere", "degrees")
lat =
 Columns 1 through 6:
   3.5082e-15   3.0000e+01   6.0000e+01   9.0000e+01   6.0000e+01   3.0000e+01
 Column 7:
   3.5082e-15

lon =
 Columns 1 through 4:
  -90.0000 +       0i  -90.0000 +       0i  -90.0000 +  0.0000i         0 +       0i
 Columns 5 through 7:
   90.0000 -  0.0000i   90.0000 +       0i   90.0000 +       0i

>> [latr, lonr] = reckon (0, 0, 90, [-90 : 30: 90])
latr =
         0   30.0000   60.0000   90.0000   60.0000   30.0000         0

lonr =
  -90  -90  -90    0   90   90   90


I'd first try to let geodeticfwd.m match reckon.m's result for a Unit Sphere (note: radius =1.000).
Latitude seem to match, longitudes don't.

Philip Nienhuis <philipnienhuis>
Group Member
Mon 22 Feb 2021 09:48:17 PM UTC, comment #4: 

Just playing around a bit.

>> [lat, lon] = geodeticfwd (0, 0, [-3600:180:3600], earthRadius * 2 * pi() / 360, "Sphere", "degrees")
lat =
 Columns 1 through 22:
   1  -1   1  -1   1  -1   1  -1   1  -1   1  -1   1  -1   1  -1   1  -1   1  -1   1  -1
 Columns 23 through 41:
   1  -1   1  -1   1  -1   1  -1   1  -1   1  -1   1  -1   1  -1   1  -1   1

lon =
 Columns 1 through 6:
   2.4495e-15   1.2261e-15   2.2045e-15   1.4710e-15   1.9596e-15   1.7160e-15
 Columns 7 through 12:
   1.7146e-15   1.9609e-15   1.4697e-15   2.2059e-15   1.2247e-15  -1.1023e-15
 Columns 13 through 18:
   9.7978e-16  -8.5731e-16   7.3484e-16  -6.1237e-16   4.8989e-16  -3.6742e-16
 Columns 19 through 24:
   2.4495e-16  -1.2247e-16            0   1.2247e-16  -2.4495e-16   3.6742e-16
 Columns 25 through 30:
  -4.8989e-16   6.1237e-16  -7.3484e-16   8.5731e-16  -9.7978e-16   1.1023e-15
 Columns 31 through 36:
  -1.2247e-15  -2.2059e-15  -1.4697e-15  -1.9609e-15  -1.7146e-15  -1.7160e-15
 Columns 37 through 41:
  -1.9596e-15  -1.4710e-15  -2.2045e-15  -1.2261e-15  -2.4495e-15


==> looks worthwile to wrapTo180 / wrapToPi the azimuth argument behind the scenes in order to limit the max error.

Then,

>> [lat, lon] = geodeticfwd (0, 0, [-360:180:360], earthRadius * 2 * pi() / 4, "Sphere", "degrees")
lat =
   90  -90   90  -90   90

lon =
 Columns 1 through 3:
    90.0000 - 118.2262i   -90.0000 +  75.4561i          0 +        0i
 Columns 4 and 5:
    90.0000 -  75.4561i   -90.0000 + 118.2262i


... nasty, complex output. I suspect this corner case (moving from equator to the poles) needs better preprocessing around Ls. 143-144.

Moving from the poles to the equator goes better:

>> [lat, lon] = geodeticfwd (90, 0, [-360:180:360], earthRadius * 2 * pi() / 4, "Sphere", "degrees")
lat =
   7.0165e-15            0   7.0165e-15            0   7.0165e-15

lon =
   1.4033e-14  -7.0165e-15            0   7.0165e-15  -1.4033e-14


This is just with a spherical Earth model, haven't tried ellipsoids yet.

Philip Nienhuis <philipnienhuis>
Group Member
Wed 17 Feb 2021 08:49:48 PM UTC, comment #3: 

Maybe geodeticfwd.m could be more intimately linked to reckon.m? that is, it does for an ellipsoid what reckonm does for a spheroid.

I'll check if geodeticfwd.m is more accurate - I think we've discussed in the patches for scxsc.m that reckon.m as it stands could/should be more accurate.

Philip Nienhuis <philipnienhuis>
Group Member
Tue 16 Feb 2021 01:17:25 PM UTC, comment #2: 

Correct there is no geodeticfwd.
It looks like this would supersede meridianfwd.
When writing meridianfwd it seemed better to write the stronger version and then create the wrapper.

Also there are many functions which check the spheroid, and it has been improved on.
Would it be better to have this be a private function?
Where you would only need to update in one place.
Something like the attached.



(file #50866)

Anonymous
Mon 15 Feb 2021 08:59:47 PM UTC, comment #1: 

AFAICS Matlab has no geodeticfwd.

meridianfwd is the inverse of meridianarc.
geodeticfwd is extension over reference ellipsoids of reckon.
Right?

Anyway both look useful to me.

Philip Nienhuis <philipnienhuis>
Group Member
Sun 14 Feb 2021 09:45:00 PM UTC, original submission:  

I think this is a simplified version of the geodetic forward.

Implemented the geodeticfwd and created a wrapper for
meridianfwd.


Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #51193:  __ell_chk__.m added by philipnienhuis (3KiB - text/plain)
file #51091:  sph_chk.m added by None (2KiB - text/x-objcsrc - version 2)
file #50866:  sph_chk.m added by None (1KiB - application/octet-stream)
file #50863:  geodeticfwd.m added by None (5KiB - application/octet-stream)
file #50864:  meridianfwd.m added by None (2KiB - 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
  •  

    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 11 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-07-16 philipnienhuis StatusReady For Test Done
        Open/ClosedOpen Closed
    2021-07-12 philipnienhuis StatusIn Progress Ready For Test
    2021-04-05 philipnienhuis Attached File- Added _ell_chk_.m, #51193
    2021-03-19 None Attached File- Added sph_chk.m, #51091
    2021-03-10 philipnienhuis StatusNone In Progress
    2021-02-22 philipnienhuis Assigned toNone philipnienhuis
        Carbon-Copy- Added philipnienhuis
    2021-02-16 None Attached File- Added sph_chk.m, #50866
    2021-02-14 None Attached File- Added geodeticfwd.m, #50863
        Attached File- Added meridianfwd.m, #50864

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code