New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix vincenty inverse for equatorial line #2282

Merged
merged 3 commits into from Jan 11, 2019

Conversation

Projects
None yet
3 participants
@Brtle
Copy link
Contributor

Brtle commented Jan 10, 2019

I looked in detail the calculation of vincenty inverse to fix the unit test on Alpine (see #1825). I think that a case is not correctly handled when cos²α = 0. On this website, a note describes this case:

Vincenty observes that eqn. (18) becomes indeterminate over equatorial lines (since cos²α → 0); in this case, set cos2σm to 0 and the result is computed correctly.

The problematic unit test calculates the distance between (0, 0) and (0, 17) which is precisely an equatorial line. Handling the cos²α = 0 fixes this issue. Another test on the equatorial line was expected to fail (distance between (0, 0) and (0, 13)) and is converging now.

Moreover, this case is also handled in geopy in the same way.

I also did some minor changes:

  • Some cos/sin was already calculated, I reused directly the value instead of recalculating the cos/sin.
  • α was calculated but only cos²α was useful, I calculated it directly without using asin.
  • This code block was inside the loop for no reason:
u2 = pow(math.cos(alpha), 2) * (a * a - b * b) / (b * b)
_a = 1 + (u2 / 16384) * (4096 + u2 * (-768 + u2 *
                                        (320 - 175 * u2)))
_b = (u2 / 1024) * (256 + u2 * (-128 + u2 * (74 - 47 * u2)))
delta_sigma = _b * sin_sigma * \
    (cos2sigma_m + (_b / 4) *
        (cos_sigma * (-1 + 2 * pow(cos2sigma_m, 2)) - (_b / 6) *
            cos2sigma_m * (-3 + 4 * sqr_sin_sigma) *
            (-3 + 4 * pow(cos2sigma_m, 2))))

dist = b * _a * (sigma - delta_sigma)
alpha12 = math.atan2(
    (math.cos(u_2) * math.sin(dlon)),
    (math.cos(u_1) * math.sin(u_2) -
        math.sin(u_1) * math.cos(u_2) * math.cos(dlon)))
alpha21 = math.atan2(
    (math.cos(u_1) * math.sin(dlon)),
    (-math.sin(u_1) * math.cos(u_2) +
        math.cos(u_1) * math.sin(u_2) * math.cos(dlon)))

It was useless to calculate these variables at each iteration, I moved the block at the end of the function.

@QuLogic
Copy link
Member

QuLogic left a comment

Makes sense to me; not sure what the state of the maintenance branch is though.

Show resolved Hide resolved obspy/geodetics/base.py
Show resolved Hide resolved obspy/geodetics/base.py Outdated
@megies

This comment has been minimized.

Copy link
Member

megies commented Jan 11, 2019

Looks good, thanks! Can you add a one-liner to the changelog as well?

@megies megies added this to the 1.1.1 milestone Jan 11, 2019

@megies

megies approved these changes Jan 11, 2019

@Brtle

This comment has been minimized.

Copy link
Contributor

Brtle commented Jan 11, 2019

Looks good, thanks! Can you add a one-liner to the changelog as well?

Oops, I forgot that. Done!

@megies

This comment has been minimized.

Copy link
Member

megies commented Jan 11, 2019

Thanks @Brtle! 🎉

@megies megies merged commit 67df2e1 into maintenance_1.1.x Jan 11, 2019

0 of 3 checks passed

ci/circleci Your tests failed on CircleCI
Details
continuous-integration/appveyor/pr AppVeyor build failed
Details
continuous-integration/travis-ci/pr The Travis CI build failed
Details

@megies megies deleted the fix_vincenty_inverse branch Jan 11, 2019

@Brtle Brtle referenced this pull request Jan 11, 2019

Closed

Tests fail on Alpine Linux #1825

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment