geovectorslib.geod.direct() - Code Metrics - omdv/geovectors - Measure and Improve Code Quality continuously with Scrutinizer

geovectorslib.geod.direct() B
last analyzed 2021-06-04 00:46 UTC

↳ Parent: geovectorslib.geod

Complexity

Conditions

Size

Total Lines	98
Code Lines	59

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	3
eloc	59
nop	4
dl	0
loc	98
rs	8.3417
c	0
b	0
f	0

How to fix Long Method

import numpy as np
from geographiclib.geodesic import Geodesic

from geovectorslib.utils import wrap90deg, wrap180deg, wrap360deg


def inverse(lats1: 'list', lons1: 'list', lats2: 'list', lons2: 'list') -> 'dict':
    """
    Inverse geodesic using Vincenty approach.

    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees)

    All args must be of equal length.

    Ref:
    https://www.movable-type.co.uk/scripts/latlong-vincenty.html
    """

    eps = np.finfo(float).eps

    lon1, lon2 = map(wrap180deg, [np.array(lons1), np.array(lons2)])
    lat1, lat2 = map(wrap90deg, [np.array(lats1), np.array(lats2)])
    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])

    # WGS84
    a = 6378137.0
    b = 6356752.314245
    f = 1.0 / 298.257223563

    # L = difference in longitude
    L = lon2 - lon1

    # U = reduced latitude, defined by tan U = (1-f)·tanφ.
    tan_U1 = (1 - f) * np.tan(lat1)
    cos_U1 = 1 / np.sqrt((1 + tan_U1 * tan_U1))
    sin_U1 = tan_U1 * cos_U1
    tan_U2 = (1 - f) * np.tan(lat2)
    cos_U2 = 1 / np.sqrt((1 + tan_U2 * tan_U2))
    sin_U2 = tan_U2 * cos_U2

    # checks for antipodal points
    antipodal = (np.abs(L) > np.pi / 2) | (np.abs(lat2 - lat1) > np.pi / 2)

    # delta = difference in longitude on an auxiliary sphere
    delta = L
    sin_delta = np.sin(delta)
    cos_delta = np.cos(delta)

    # sigma = angular distance P₁ P₂ on the sphere
    sigma = np.zeros(lat1.shape)
    sigma[antipodal] = np.pi

    cos_sigma = np.ones(lat1.shape)
    cos_sigma[antipodal] = -1

    # sigmam = angular distance on the sphere from the equator
    # to the midpoint of the line
    # azi = azimuth of the geodesic at the equator
    # init delta_prime at very high value to not miss iterations when delta=0
    delta_prime = 1 / eps * np.ones(lat1.shape)
    iterations = 0

    # init before loop to allow mask indexing
    sin_Sq_sigma = np.zeros(lat1.shape)
    sin_sigma = np.zeros(lat1.shape)
    cos_sigma = np.zeros(lat1.shape)
    sin_azi = np.zeros(lat1.shape)
    cos_Sq_azi = np.ones(lat1.shape)
    cos_2_sigma_m = np.ones(lat1.shape)
    C = np.ones(lat1.shape)

    # init mask
    m = np.ones(lat1.shape, dtype=bool)
    conv_mask = np.zeros(lat1.shape, dtype=bool)

    while (np.abs(delta[m] - delta_prime[m]) > 1e-12).any():
        sin_delta = np.sin(delta)
        cos_delta = np.cos(delta)
        sin_Sq_sigma = (cos_U2 * sin_delta) * (cos_U2 * sin_delta) + (
            cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_delta
        ) * (cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_delta)

        # co-incident/antipodal points mask - exclude from the rest of the loop
        # the value has to be about 1.e-4 experimentally
        # otherwise there are issues with convergence for near antipodal points
        m = (np.abs(sin_Sq_sigma) > eps) & (sin_Sq_sigma != np.nan)

        sin_sigma[m] = np.sqrt(sin_Sq_sigma[m])
        cos_sigma[m] = sin_U1[m] * sin_U2[m] + cos_U1[m] * cos_U2[m] * cos_delta[m]
        sigma[m] = np.arctan2(sin_sigma[m], cos_sigma[m])
        sin_azi[m] = cos_U1[m] * cos_U2[m] * sin_delta[m] / sin_sigma[m]
        cos_Sq_azi[m] = 1 - sin_azi[m] * sin_azi[m]

        # on equatorial line cos²azi = 0
        cos_2_sigma_m[m] = cos_sigma[m] - 2 * sin_U1[m] * sin_U2[m] / cos_Sq_azi[m]
        cos_2_sigma_m[cos_Sq_azi == 0] = 0

        C[m] = f / 16 * cos_Sq_azi[m] * (4 + f * (4 - 3 * cos_Sq_azi[m]))

        delta_prime = delta
        delta = L + (1 - C) * f * sin_azi * (
            sigma
            + (C * sin_sigma)
            * (cos_2_sigma_m + C * cos_sigma * (-1 + 2 * cos_2_sigma_m * cos_2_sigma_m))
        )

        # Exceptions
        iteration_check = np.abs(delta)
        iteration_check[antipodal] = iteration_check[antipodal] - np.pi
        if (iteration_check > np.pi).any():
            raise Exception('delta > np.pi')

        iterations += 1
        if iterations >= 50:
            conv_mask = np.abs(delta - delta_prime) > 1e-12
            break

    uSq = cos_Sq_azi * (a * a - b * b) / (b * b)
    A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
    B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
    d_sigma = (B * sin_sigma) * (
        cos_2_sigma_m
        + (B / 4)
        * (
            cos_sigma * (-1 + 2 * cos_2_sigma_m * cos_2_sigma_m)
            - (B / 6 * cos_2_sigma_m)
            * (-3 + 4 * sin_sigma * sin_sigma)
            * (-3 + 4 * cos_2_sigma_m * cos_2_sigma_m)
        )
    )
    # length of the geodesic
    s12 = b * A * (sigma - d_sigma)

    # note special handling of exactly antipodal points where sin²sigma = 0
    # (due to discontinuity atan2(0, 0) = 0 but atan2(eps, 0) = np.pi/2 / 90°)
    # in which case bearing is always meridional, due north (or due south!)

    azi1 = np.arctan2(cos_U2 * sin_delta, cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_delta)
    azi1[np.abs(sin_Sq_sigma) < eps] = 0

    azi2 = np.arctan2(
        cos_U1 * sin_delta, -sin_U1 * cos_U2 + cos_U1 * sin_U2 * cos_delta
    )
    azi2[np.abs(sin_Sq_sigma) < eps] = np.pi

    azi1 = wrap360deg(azi1 * 180 / np.pi)
    azi2 = wrap360deg(azi2 * 180 / np.pi)

    # if distance is too small - return 0 instead of `None`
    # to be consistent with geographiclib
    mask = s12 < eps
    azi1[mask] = 0
    azi2[mask] = 0

    # use geographiclib for points which didn't converge
    if conv_mask[conv_mask].shape[0] > 0:
        vInverse = np.vectorize(Geodesic.WGS84.Inverse)
        _ps = vInverse(
            np.array(lats1)[conv_mask],
            np.array(lons1)[conv_mask],
            np.array(lats2)[conv_mask],
            np.array(lons2)[conv_mask],
        )
        s12[conv_mask] = [_p['s12'] for _p in _ps]
        azi1[conv_mask] = [_p['azi1'] for _p in _ps]
        azi2[conv_mask] = [_p['azi2'] for _p in _ps]

    return {'s12': s12, 'azi1': azi1, 'azi2': azi2, 'iterations': iterations}


def direct(lats1: 'list', lons1: 'list', brgs: 'list', dists: 'list') -> 'dict':
    """
    Direct geodesic using Vincenty approach.

    Given a start point, initial bearing (0° is North, clockwise)
    and distance in meters, this will calculate the destination point
    and final bearing travelling along a (shortest distance) great circle arc.

    Bearing in 0-360deg starting from North clockwise.

    all variables must be of same length

    Ref:
    https://www.movable-type.co.uk/scripts/latlong.html
    https://www.movable-type.co.uk/scripts/latlong-vincenty.html
    """

    lon1, lat1, brg = map(np.radians, [lons1, lats1, brgs])

    # WGS84
    a = 6378137.0
    b = 6356752.314245
    f = 1.0 / 298.257223563

    azi1 = brg
    s = dists

    sin_azi1 = np.sin(azi1)
    cos_azi1 = np.cos(azi1)

    tan_U1 = (1 - f) * np.tan(lat1)
    cos_U1 = 1 / np.sqrt((1 + tan_U1 * tan_U1))
    sin_U1 = tan_U1 * cos_U1

    # sigma1 = angular distance on the sphere from the equator to P1
    sigma1 = np.arctan2(tan_U1, cos_azi1)

    sin_azi = cos_U1 * sin_azi1  # azi = azimuth of the geodesic at the equator
    cos_Sq_azi = 1 - sin_azi * sin_azi
    uSq = cos_Sq_azi * (a * a - b * b) / (b * b)
    A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
    B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))

    sigma = s / (b * A)
    # sigma = angular distance P₁ P₂ on the sphere
    # sigmam = angular distance on the sphere from the equator to
    # the midpoint of the line
    sin_sigma = np.sin(sigma)
    cos_sigma = np.cos(sigma)
    cos_2_sigma_m = np.cos(2 * sigma1 + sigma)

    iterations = 0
    sigma_prime = 0

    while (np.abs(sigma - sigma_prime) > 1e-12).any():
        cos_2_sigma_m = np.cos(2 * sigma1 + sigma)
        sin_sigma = np.sin(sigma)
        cos_sigma = np.cos(sigma)
        delta_sigma = (B * sin_sigma) * (
            cos_2_sigma_m
            + B
            / 4
            * (
                cos_sigma * (-1 + 2 * cos_2_sigma_m * cos_2_sigma_m)
                - (B / 6 * cos_2_sigma_m)
                * (-3 + 4 * sin_sigma * sin_sigma)
                * (-3 + 4 * cos_2_sigma_m * cos_2_sigma_m)
            )
        )
        sigma_prime = sigma
        sigma = s / (b * A) + delta_sigma
        iterations += 1
        if iterations >= 50:
            raise Exception('Vincenty formula failed to converge')

    x = sin_U1 * sin_sigma - cos_U1 * cos_sigma * cos_azi1
    lat2 = np.arctan2(
        sin_U1 * cos_sigma + cos_U1 * sin_sigma * cos_azi1,
        (1 - f) * np.sqrt(sin_azi * sin_azi + x * x),
    )
    lambda_ = np.arctan2(
        sin_sigma * sin_azi1, cos_U1 * cos_sigma - sin_U1 * sin_sigma * cos_azi1
    )
    C = f / 16 * cos_Sq_azi * (4 + f * (4 - 3 * cos_Sq_azi))
    L = lambda_ - (1 - C) * f * sin_azi * (
        sigma
        + C
        * sin_sigma
        * (cos_2_sigma_m + C * cos_sigma * (-1 + 2 * cos_2_sigma_m * cos_2_sigma_m))
    )
    lon2 = lon1 + L
    azi2 = np.arctan2(sin_azi, -x)

    lat2 = wrap90deg(lat2 * 180.0 / np.pi)
    lon2 = wrap180deg(lon2 * 180.0 / np.pi)
    azi2 = wrap360deg(azi2 * 180.0 / np.pi)

    return {'lat2': lat2, 'lon2': lon2, 'azi2': azi2, 'iterations': iterations}


1			import numpy as np
2			from geographiclib.geodesic import Geodesic
3
4			from geovectorslib.utils import wrap90deg, wrap180deg, wrap360deg
5
6
7			def inverse(lats1: 'list', lons1: 'list', lats2: 'list', lons2: 'list') -> 'dict':
8			"""
9			Inverse geodesic using Vincenty approach.
10
11			Calculate the great circle distance between two points
12			on the earth (specified in decimal degrees)
13
14			All args must be of equal length.
15
16			Ref:
17			https://www.movable-type.co.uk/scripts/latlong-vincenty.html
18			"""
19
20			eps = np.finfo(float).eps
21
22			lon1, lon2 = map(wrap180deg, [np.array(lons1), np.array(lons2)])
23			lat1, lat2 = map(wrap90deg, [np.array(lats1), np.array(lats2)])
24			lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
25
26			# WGS84
27			a = 6378137.0
28			b = 6356752.314245
29			f = 1.0 / 298.257223563
30
31			# L = difference in longitude
32			L = lon2 - lon1
33
34			# U = reduced latitude, defined by tan U = (1-f)·tanφ.
35			tan_U1 = (1 - f) * np.tan(lat1)
36			cos_U1 = 1 / np.sqrt((1 + tan_U1 * tan_U1))
37			sin_U1 = tan_U1 * cos_U1
38			tan_U2 = (1 - f) * np.tan(lat2)
39			cos_U2 = 1 / np.sqrt((1 + tan_U2 * tan_U2))
40			sin_U2 = tan_U2 * cos_U2
41
42			# checks for antipodal points
43			antipodal = (np.abs(L) > np.pi / 2) \| (np.abs(lat2 - lat1) > np.pi / 2)
44
45			# delta = difference in longitude on an auxiliary sphere
46			delta = L
47			sin_delta = np.sin(delta)
48			cos_delta = np.cos(delta)
49
50			# sigma = angular distance P₁ P₂ on the sphere
51			sigma = np.zeros(lat1.shape)
52			sigma[antipodal] = np.pi
53
54			cos_sigma = np.ones(lat1.shape)
55			cos_sigma[antipodal] = -1
56
57			# sigmam = angular distance on the sphere from the equator
58			# to the midpoint of the line
59			# azi = azimuth of the geodesic at the equator
60			# init delta_prime at very high value to not miss iterations when delta=0
61			delta_prime = 1 / eps * np.ones(lat1.shape)
62			iterations = 0
63
64			# init before loop to allow mask indexing
65			sin_Sq_sigma = np.zeros(lat1.shape)
66			sin_sigma = np.zeros(lat1.shape)
67			cos_sigma = np.zeros(lat1.shape)
68			sin_azi = np.zeros(lat1.shape)
69			cos_Sq_azi = np.ones(lat1.shape)
70			cos_2_sigma_m = np.ones(lat1.shape)
71			C = np.ones(lat1.shape)
72
73			# init mask
74			m = np.ones(lat1.shape, dtype=bool)
75			conv_mask = np.zeros(lat1.shape, dtype=bool)
76
77			while (np.abs(delta[m] - delta_prime[m]) > 1e-12).any():
78			sin_delta = np.sin(delta)
79			cos_delta = np.cos(delta)
80			sin_Sq_sigma = (cos_U2 * sin_delta) * (cos_U2 * sin_delta) + (
81			cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_delta
82			) * (cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_delta)
83
84			# co-incident/antipodal points mask - exclude from the rest of the loop
85			# the value has to be about 1.e-4 experimentally
86			# otherwise there are issues with convergence for near antipodal points
87			m = (np.abs(sin_Sq_sigma) > eps) & (sin_Sq_sigma != np.nan)
88
89			sin_sigma[m] = np.sqrt(sin_Sq_sigma[m])
90			cos_sigma[m] = sin_U1[m] * sin_U2[m] + cos_U1[m] * cos_U2[m] * cos_delta[m]
91			sigma[m] = np.arctan2(sin_sigma[m], cos_sigma[m])
92			sin_azi[m] = cos_U1[m] * cos_U2[m] * sin_delta[m] / sin_sigma[m]
93			cos_Sq_azi[m] = 1 - sin_azi[m] * sin_azi[m]
94
95			# on equatorial line cos²azi = 0
96			cos_2_sigma_m[m] = cos_sigma[m] - 2 * sin_U1[m] * sin_U2[m] / cos_Sq_azi[m]
97			cos_2_sigma_m[cos_Sq_azi == 0] = 0
98
99			C[m] = f / 16 * cos_Sq_azi[m] * (4 + f * (4 - 3 * cos_Sq_azi[m]))
100
101			delta_prime = delta
102			delta = L + (1 - C) * f * sin_azi * (
103			sigma
104			+ (C * sin_sigma)
105			* (cos_2_sigma_m + C * cos_sigma * (-1 + 2 * cos_2_sigma_m * cos_2_sigma_m))
106			)
107
108			# Exceptions
109			iteration_check = np.abs(delta)
110			iteration_check[antipodal] = iteration_check[antipodal] - np.pi
111			if (iteration_check > np.pi).any():
112			raise Exception('delta > np.pi')
113
114			iterations += 1
115			if iterations >= 50:
116			conv_mask = np.abs(delta - delta_prime) > 1e-12
117			break
118
119			uSq = cos_Sq_azi * (a * a - b * b) / (b * b)
120			A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
121			B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
122			d_sigma = (B * sin_sigma) * (
123			cos_2_sigma_m
124			+ (B / 4)
125			* (
126			cos_sigma * (-1 + 2 * cos_2_sigma_m * cos_2_sigma_m)
127			- (B / 6 * cos_2_sigma_m)
128			* (-3 + 4 * sin_sigma * sin_sigma)
129			* (-3 + 4 * cos_2_sigma_m * cos_2_sigma_m)
130			)
131			)
132			# length of the geodesic
133			s12 = b * A * (sigma - d_sigma)
134
135			# note special handling of exactly antipodal points where sin²sigma = 0
136			# (due to discontinuity atan2(0, 0) = 0 but atan2(eps, 0) = np.pi/2 / 90°)
137			# in which case bearing is always meridional, due north (or due south!)
138
139			azi1 = np.arctan2(cos_U2 * sin_delta, cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_delta)
140			azi1[np.abs(sin_Sq_sigma) < eps] = 0
141
142			azi2 = np.arctan2(
143			cos_U1 * sin_delta, -sin_U1 * cos_U2 + cos_U1 * sin_U2 * cos_delta
144			)
145			azi2[np.abs(sin_Sq_sigma) < eps] = np.pi
146
147			azi1 = wrap360deg(azi1 * 180 / np.pi)
148			azi2 = wrap360deg(azi2 * 180 / np.pi)
149
150			# if distance is too small - return 0 instead of `None`
151			# to be consistent with geographiclib
152			mask = s12 < eps
153			azi1[mask] = 0
154			azi2[mask] = 0
155
156			# use geographiclib for points which didn't converge
157			if conv_mask[conv_mask].shape[0] > 0:
158			vInverse = np.vectorize(Geodesic.WGS84.Inverse)
159			_ps = vInverse(
160			np.array(lats1)[conv_mask],
161			np.array(lons1)[conv_mask],
162			np.array(lats2)[conv_mask],
163			np.array(lons2)[conv_mask],
164			)
165			s12[conv_mask] = [_p['s12'] for _p in _ps]
166			azi1[conv_mask] = [_p['azi1'] for _p in _ps]
167			azi2[conv_mask] = [_p['azi2'] for _p in _ps]
168
169			return {'s12': s12, 'azi1': azi1, 'azi2': azi2, 'iterations': iterations}
170
171
172			def direct(lats1: 'list', lons1: 'list', brgs: 'list', dists: 'list') -> 'dict':
173			"""
174			Direct geodesic using Vincenty approach.
175
176			Given a start point, initial bearing (0° is North, clockwise)
177			and distance in meters, this will calculate the destination point
178			and final bearing travelling along a (shortest distance) great circle arc.
179
180			Bearing in 0-360deg starting from North clockwise.
181
182			all variables must be of same length
183
184			Ref:
185			https://www.movable-type.co.uk/scripts/latlong.html
186			https://www.movable-type.co.uk/scripts/latlong-vincenty.html
187			"""
188
189			lon1, lat1, brg = map(np.radians, [lons1, lats1, brgs])
190
191			# WGS84
192			a = 6378137.0
193			b = 6356752.314245
194			f = 1.0 / 298.257223563
195
196			azi1 = brg
197			s = dists
198
199			sin_azi1 = np.sin(azi1)
200			cos_azi1 = np.cos(azi1)
201
202			tan_U1 = (1 - f) * np.tan(lat1)
203			cos_U1 = 1 / np.sqrt((1 + tan_U1 * tan_U1))
204			sin_U1 = tan_U1 * cos_U1
205
206			# sigma1 = angular distance on the sphere from the equator to P1
207			sigma1 = np.arctan2(tan_U1, cos_azi1)
208
209			sin_azi = cos_U1 * sin_azi1 # azi = azimuth of the geodesic at the equator
210			cos_Sq_azi = 1 - sin_azi * sin_azi
211			uSq = cos_Sq_azi * (a * a - b * b) / (b * b)
212			A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
213			B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
214
215			sigma = s / (b * A)
216			# sigma = angular distance P₁ P₂ on the sphere
217			# sigmam = angular distance on the sphere from the equator to
218			# the midpoint of the line
219			sin_sigma = np.sin(sigma)
220			cos_sigma = np.cos(sigma)
221			cos_2_sigma_m = np.cos(2 * sigma1 + sigma)
222
223			iterations = 0
224			sigma_prime = 0
225
226			while (np.abs(sigma - sigma_prime) > 1e-12).any():
227			cos_2_sigma_m = np.cos(2 * sigma1 + sigma)
228			sin_sigma = np.sin(sigma)
229			cos_sigma = np.cos(sigma)
230			delta_sigma = (B * sin_sigma) * (
231			cos_2_sigma_m
232			+ B
233			/ 4
234			* (
235			cos_sigma * (-1 + 2 * cos_2_sigma_m * cos_2_sigma_m)
236			- (B / 6 * cos_2_sigma_m)
237			* (-3 + 4 * sin_sigma * sin_sigma)
238			* (-3 + 4 * cos_2_sigma_m * cos_2_sigma_m)
239			)
240			)
241			sigma_prime = sigma
242			sigma = s / (b * A) + delta_sigma
243			iterations += 1
244			if iterations >= 50:
245			raise Exception('Vincenty formula failed to converge')
246
247			x = sin_U1 * sin_sigma - cos_U1 * cos_sigma * cos_azi1
248			lat2 = np.arctan2(
249			sin_U1 * cos_sigma + cos_U1 * sin_sigma * cos_azi1,
250			(1 - f) * np.sqrt(sin_azi * sin_azi + x * x),
251			)
252			lambda_ = np.arctan2(
253			sin_sigma * sin_azi1, cos_U1 * cos_sigma - sin_U1 * sin_sigma * cos_azi1
254			)
255			C = f / 16 * cos_Sq_azi * (4 + f * (4 - 3 * cos_Sq_azi))
256			L = lambda_ - (1 - C) * f * sin_azi * (
257			sigma
258			+ C
259			* sin_sigma
260			* (cos_2_sigma_m + C * cos_sigma * (-1 + 2 * cos_2_sigma_m * cos_2_sigma_m))
261			)
262			lon2 = lon1 + L
263			azi2 = np.arctan2(sin_azi, -x)
264
265			lat2 = wrap90deg(lat2 * 180.0 / np.pi)
266			lon2 = wrap180deg(lon2 * 180.0 / np.pi)
267			azi2 = wrap360deg(azi2 * 180.0 / np.pi)
268
269			return {'lat2': lat2, 'lon2': lon2, 'azi2': azi2, 'iterations': iterations}
270

omdv / geovectors

geovectorslib.geod.direct() B last analyzed 2021-06-04 00:46 UTC

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like

geovectorslib.geod.direct() B
last analyzed 2021-06-04 00:46 UTC