Re: [gpsd-dev] CEP(50) 0.009 meters

[Fred Wright]

On Mon, 25 Jul 2016, Gary E. Miller wrote:
> On Mon, 25 Jul 2016 12:30:07 -0700
> Hal Murray <address@hidden> wrote:
>
> > address@hidden said:
> > > Looks like I have a bit more work to do, but close.  Changes so far
> > > are in git head.  Next I'll add WGS84 ellipsoid correction for
> > > latitude.
> >
> > Why?  Is that going to be significant?
>
> It could be.  Skytraq claims, not yet verified, that they can get
> cEnti-meter accuracy off a 1km baseline.  At that distance the
> slightly less than 1% error would be noticeable.  At least to someone
> very picky.

At your latitude it's about 0.27%.  See below.

> Even if not important, if the last guy that worked on gpsprof had
> worked to rediculous overkill in the math then I would not have had to
> suffer to fix it.  So I'm paying it forward to the next guy pushing the
> envelope.
>
> > Assuming you start with doubles, after you do the subtract, how many
> > bits of precision do you have left?
>
> Right now I am subtracting these two:
>
>   121.31405412
> - 121.31405411
> --------------
>
>     0.00000001
>
> I figure the latitude is 35 binary digits of precision.  Plus a sign bit.
>
> An IEEE 754 double float is 53 bits of binary precision.  Plus a sign bit.
>
> So that leaves 18 bits binary precision to play with.  About 5.5 decimal
> digits of precision.
>
> If I only know the meters inn a degree (lat or lon) to about 1% that
> seems gross compared to the rest.
>
> The next step, that I'll not take, would be altitude correction.  I've been
> reading PhD papers that do indeed go that far.

Strictly speaking, for this sort of calculation you really want
straight-line distance, not great-circle distance, though the difference
for short distances is miniscule.  OTOH, that calculation is
mathematically exact and reasonably numerically accurate, without needing
iteration (beyond the internal iteration in the trig functions, that is).
Basically, you just convert the points to ECEF coordinates and then do the
usual Pythagorean distance calculation (3D version).

I'm attaching some code I started working on a while back to do this sort
of thing.  It's more general in that it works for arbitrary ellipsoids,
and also allows you to specify the value of PI (since the GPS specs call
out a specific PI value that one is supposed to use for consistency).

Note that everything calculated in Ellipsoid.__new__() is constant for any
given ellipsoid (and PI value).

Sample:

MacPro:Util fw$ python
Python 2.7.12 (default, Jun 29 2016, 12:46:54)
[GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import geo_utils as gu
>>> e = gu.Ellipsoid.ByName('wgs84')
>>> e.DistElev(44.058188, 121.31405412, 3629 / .3048, 44.058188, 121.31405411, 
>>> 3629 / .3048)
(0.000802770017091834, 0.0)
>>> gu.math.sqrt(e._GetLocParams(44.058188, 121.31405412, 3629 / 
>>> .3048).norm_sqr)
6388485.687314996
>>> gu.math.sqrt(e._GetLocParams(44.058188, 121.31405412, 3629 / 
>>> .3048).norm_sqr) / 6371000
1.0027445749984296
>>>

You didn't post your latitude (or altitude), so I used what Google Earth
gave me for central Bend.

The first result is the distance in meters, and the elevation angle (0.0
in this case).

The second result is the length of the ellipsoidal normal in meters at
your latitude.

The third result divides this by the usually accepted mean radius value to
show the discrepancy between the ellipsoidal and spherical approaches.

Fred Wright

geo_utils.py
Description: Some geographic calculation utilities