geovectorslib.geod.direct() - Code Metrics - Inspection of "Inverse and Direct functions release" - omdv/geovectors - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Push — master ( 9b4d36...ed671d )

by Oleg

created 2020-01-04 04:48 UTC

geovectorslib.geod.direct() B

↳ Parent: geovectorslib.geod

Complexity

Conditions

Size

Total Lines	101
Code Lines	60

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	3
eloc	60
nop	4
dl	0
loc	101
rs	8.309
c	0
b	0
f	0

How to fix Long Method

import numpy as np

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φ.
    tanU1 = (1 - f) * np.tan(lat1)
    cosU1 = 1 / np.sqrt((1 + tanU1 * tanU1))
    sinU1 = tanU1 * cosU1
    tanU2 = (1 - f) * np.tan(lat2)
    cosU2 = 1 / np.sqrt((1 + tanU2 * tanU2))
    sinU2 = tanU2 * cosU2

    # 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

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

    cossigma = np.ones(lat1.shape)
    cossigma[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
    deltaprime = np.zeros(lat1.shape)
    iterations = 0

    # init before loop to allow mask indexing
    sinSqsigma = np.zeros(lat1.shape)
    sinsigma = np.zeros(lat1.shape)
    cossigma = np.zeros(lat1.shape)
    sinazi = np.zeros(lat1.shape)
    cosSqazi = np.ones(lat1.shape)
    cos2sigmam = np.ones(lat1.shape)
    C = np.ones(lat1.shape)

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

    while (np.abs(delta[m] - deltaprime[m]) > 1e-12).any():
        sindelta = np.sin(delta)
        cosdelta = np.cos(delta)
        sinSqsigma = (cosU2 * sindelta) * (cosU2 * sindelta) + (
            cosU1 * sinU2 - sinU1 * cosU2 * cosdelta
        ) * (cosU1 * sinU2 - sinU1 * cosU2 * cosdelta)

        # co-incident/antipodal points mask - exclude from the rest of the loop
        m = (np.abs(sinSqsigma) > eps) & (sinSqsigma != np.nan)

        sinsigma[m] = np.sqrt(sinSqsigma[m])
        cossigma[m] = sinU1[m] * sinU2[m] + cosU1[m] * cosU2[m] * cosdelta[m]
        sigma[m] = np.arctan2(sinsigma[m], cossigma[m])
        sinazi[m] = cosU1[m] * cosU2[m] * sindelta[m] / sinsigma[m]
        cosSqazi[m] = 1 - sinazi[m] * sinazi[m]

        # on equatorial line cos²azi = 0
        cos2sigmam[m] = cossigma[m] - 2 * sinU1[m] * sinU2[m] / cosSqazi[m]
        cos2sigmam[cosSqazi == 0] = 0

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

        deltaprime = delta
        delta = L + (1 - C) * f * sinazi * (
            sigma
            + C
            * sinsigma
            * (cos2sigmam + C * cossigma * (-1 + 2 * cos2sigmam * cos2sigmam))
        )

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

        iterations += 1
        if iterations >= 20:
            raise Exception('Vincenty formula failed to converge')

    uSq = cosSqazi * (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)))
    dsigma = (
        B
        * sinsigma
        * (
            cos2sigmam
            + B
            / 4
            * (
                cossigma * (-1 + 2 * cos2sigmam * cos2sigmam)
                - B
                / 6
                * cos2sigmam
                * (-3 + 4 * sinsigma * sinsigma)
                * (-3 + 4 * cos2sigmam * cos2sigmam)
            )
        )
    )
    # length of the geodesic
    s12 = b * A * (sigma - dsigma)

    # 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(cosU2 * sindelta, cosU1 * sinU2 - sinU1 * cosU2 * cosdelta)

    azi1[np.abs(sinSqsigma) < eps] = 0

    azi2 = np.arctan2(cosU1 * sindelta, -sinU1 * cosU2 + cosU1 * sinU2 * cosdelta)
    azi2[np.abs(sinSqsigma) < eps] = np.pi

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

    # if distance is too small
    mask = s12 < eps
    azi1[mask] = None
    azi2[mask] = None

    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

    sinazi1 = np.sin(azi1)
    cosazi1 = np.cos(azi1)

    tanU1 = (1 - f) * np.tan(lat1)
    cosU1 = 1 / np.sqrt((1 + tanU1 * tanU1))
    sinU1 = tanU1 * cosU1

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

    sinazi = cosU1 * sinazi1  # azi = azimuth of the geodesic at the equator
    cosSqazi = 1 - sinazi * sinazi
    uSq = cosSqazi * (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

    iterations = 0
    sigmaprime = 0

    while (np.abs(sigma - sigmaprime) > 1e-12).any():
        cos2sigmam = np.cos(2 * sigma1 + sigma)
        sinsigma = np.sin(sigma)
        cossigma = np.cos(sigma)
        deltasigma = (
            B
            * sinsigma
            * (
                cos2sigmam
                + B
                / 4
                * (
                    cossigma * (-1 + 2 * cos2sigmam * cos2sigmam)
                    - B
                    / 6
                    * cos2sigmam
                    * (-3 + 4 * sinsigma * sinsigma)
                    * (-3 + 4 * cos2sigmam * cos2sigmam)
                )
            )
        )
        sigmaprime = sigma
        sigma = s / (b * A) + deltasigma
        iterations += 1
        if iterations >= 100:
            raise Exception('Vincenty formula failed to converge')

    x = sinU1 * sinsigma - cosU1 * cossigma * cosazi1

    lat2 = np.arctan2(
        sinU1 * cossigma + cosU1 * sinsigma * cosazi1,
        (1 - f) * np.sqrt(sinazi * sinazi + x * x),
    )
    lambd = np.arctan2(
        sinsigma * sinazi1, cosU1 * cossigma - sinU1 * sinsigma * cosazi1
    )
    C = f / 16 * cosSqazi * (4 + f * (4 - 3 * cosSqazi))
    L = lambd - (1 - C) * f * sinazi * (
        sigma
        + C
        * sinsigma
        * (cos2sigmam + C * cossigma * (-1 + 2 * cos2sigmam * cos2sigmam))

    )
    lon2 = lon1 + L
    azi2 = np.arctan2(sinazi, -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
3			from geovectorslib.utils import wrap90deg, wrap180deg, wrap360deg
4
5
6			def inverse(lats1: 'list', lons1: 'list', lats2: 'list', lons2: 'list') -> 'dict':
7			"""
8			Inverse geodesic using Vincenty approach.
9
10			Calculate the great circle distance between two points
11			on the earth (specified in decimal degrees)
12
13			All args must be of equal length.
14
15			Ref:
16			https://www.movable-type.co.uk/scripts/latlong-vincenty.html
17			"""
18
19			eps = np.finfo(float).eps
20
21			lon1, lon2 = map(wrap180deg, [np.array(lons1), np.array(lons2)])
22			lat1, lat2 = map(wrap90deg, [np.array(lats1), np.array(lats2)])
23			lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
24
25			# WGS84
26			a = 6378137.0
27			b = 6356752.314245
28			f = 1.0 / 298.257223563
29
30			# L = difference in longitude
31			L = lon2 - lon1
32
33			# U = reduced latitude, defined by tan U = (1-f)·tanφ.
34			tanU1 = (1 - f) * np.tan(lat1)
35			cosU1 = 1 / np.sqrt((1 + tanU1 * tanU1))
36			sinU1 = tanU1 * cosU1
37			tanU2 = (1 - f) * np.tan(lat2)
38			cosU2 = 1 / np.sqrt((1 + tanU2 * tanU2))
39			sinU2 = tanU2 * cosU2
40
41			# checks for antipodal points
42			antipodal = (np.abs(L) > np.pi / 2) \| (np.abs(lat2 - lat1) > np.pi / 2)
43
44			# delta = difference in longitude on an auxiliary sphere
45			delta = L
46
47			# sigma = angular distance P₁ P₂ on the sphere
48			sigma = np.zeros(lat1.shape)
49			sigma[antipodal] = np.pi
50
51			cossigma = np.ones(lat1.shape)
52			cossigma[antipodal] = -1
53
54			# sigmam = angular distance on the sphere from the equator
55			# to the midpoint of the line
56			# azi = azimuth of the geodesic at the equator
57			deltaprime = np.zeros(lat1.shape)
58			iterations = 0
59
60			# init before loop to allow mask indexing
61			sinSqsigma = np.zeros(lat1.shape)
62			sinsigma = np.zeros(lat1.shape)
63			cossigma = np.zeros(lat1.shape)
64			sinazi = np.zeros(lat1.shape)
65			cosSqazi = np.ones(lat1.shape)
66			cos2sigmam = np.ones(lat1.shape)
67			C = np.ones(lat1.shape)
68
69			# init mask
70			m = np.ones(lat1.shape, dtype=bool)
71
72			while (np.abs(delta[m] - deltaprime[m]) > 1e-12).any():
73			sindelta = np.sin(delta)
74			cosdelta = np.cos(delta)
75			sinSqsigma = (cosU2 * sindelta) * (cosU2 * sindelta) + (
76			cosU1 * sinU2 - sinU1 * cosU2 * cosdelta
77			) * (cosU1 * sinU2 - sinU1 * cosU2 * cosdelta)
78
79			# co-incident/antipodal points mask - exclude from the rest of the loop
80			m = (np.abs(sinSqsigma) > eps) & (sinSqsigma != np.nan)
81
82			sinsigma[m] = np.sqrt(sinSqsigma[m])
83			cossigma[m] = sinU1[m] * sinU2[m] + cosU1[m] * cosU2[m] * cosdelta[m]
84			sigma[m] = np.arctan2(sinsigma[m], cossigma[m])
85			sinazi[m] = cosU1[m] * cosU2[m] * sindelta[m] / sinsigma[m]
86			cosSqazi[m] = 1 - sinazi[m] * sinazi[m]
87
88			# on equatorial line cos²azi = 0
89			cos2sigmam[m] = cossigma[m] - 2 * sinU1[m] * sinU2[m] / cosSqazi[m]
90			cos2sigmam[cosSqazi == 0] = 0
91
92			C[m] = f / 16 * cosSqazi[m] * (4 + f * (4 - 3 * cosSqazi[m]))
93
94			deltaprime = delta
95			delta = L + (1 - C) * f * sinazi * (
96			sigma
97			+ C
98			* sinsigma
99			* (cos2sigmam + C * cossigma * (-1 + 2 * cos2sigmam * cos2sigmam))
100			)
101
102			# Exceptions
103			iterationCheck = np.abs(delta)
104			iterationCheck[antipodal] = iterationCheck[antipodal] - np.pi
105			if (iterationCheck > np.pi).any():
106			raise Exception('delta > np.pi')
107
108			iterations += 1
109			if iterations >= 20:
110			raise Exception('Vincenty formula failed to converge')
111
112			uSq = cosSqazi * (a * a - b * b) / (b * b)
113			A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
114			B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
115			dsigma = (
116			B
117			* sinsigma
118			* (
119			cos2sigmam
120			+ B
121			/ 4
122			* (
123			cossigma * (-1 + 2 * cos2sigmam * cos2sigmam)
124			- B
125			/ 6
126			* cos2sigmam
127			* (-3 + 4 * sinsigma * sinsigma)
128			* (-3 + 4 * cos2sigmam * cos2sigmam)
129			)
130			)
131			)
132			# length of the geodesic
133			s12 = b * A * (sigma - dsigma)
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(cosU2 * sindelta, cosU1 * sinU2 - sinU1 * cosU2 * cosdelta)
			0 ignored issues – show introduced 2020-01-04 04:49 UTC by Report Bug Copy Issue Report The variable `sindelta` does not seem to be defined in case the `while` loop on line `72` is not entered. Are you sure this can never be the case? Loading history... introduced 2020-01-04 04:49 UTC by Report Bug Copy Issue Report The variable `cosdelta` does not seem to be defined in case the `while` loop on line `72` is not entered. Are you sure this can never be the case? Loading history...
140			azi1[np.abs(sinSqsigma) < eps] = 0
141
142			azi2 = np.arctan2(cosU1 * sindelta, -sinU1 * cosU2 + cosU1 * sinU2 * cosdelta)
143			azi2[np.abs(sinSqsigma) < eps] = np.pi
144
145			azi1 = wrap360deg(azi1 * 180 / np.pi)
146			azi2 = wrap360deg(azi2 * 180 / np.pi)
147
148			# if distance is too small
149			mask = s12 < eps
150			azi1[mask] = None
151			azi2[mask] = None
152
153			return {'s12': s12, 'azi1': azi1, 'azi2': azi2, 'iterations': iterations}
154
155
156			def direct(lats1: 'list', lons1: 'list', brgs: 'list', dists: 'list') -> 'dict':
157			"""
158			Direct geodesic using Vincenty approach.
159
160			Given a start point, initial bearing (0° is North, clockwise)
161			and distance in meters, this will calculate the destination point
162			and final bearing travelling along a (shortest distance) great circle arc.
163
164			Bearing in 0-360deg starting from North clockwise.
165
166			all variables must be of same length
167
168			Ref:
169			https://www.movable-type.co.uk/scripts/latlong.html
170			https://www.movable-type.co.uk/scripts/latlong-vincenty.html
171			"""
172
173			lon1, lat1, brg = map(np.radians, [lons1, lats1, brgs])
174
175			# WGS84
176			a = 6378137.0
177			b = 6356752.314245
178			f = 1.0 / 298.257223563
179
180			azi1 = brg
181			s = dists
182
183			sinazi1 = np.sin(azi1)
184			cosazi1 = np.cos(azi1)
185
186			tanU1 = (1 - f) * np.tan(lat1)
187			cosU1 = 1 / np.sqrt((1 + tanU1 * tanU1))
188			sinU1 = tanU1 * cosU1
189
190			# sigma1 = angular distance on the sphere from the equator to P1
191			sigma1 = np.arctan2(tanU1, cosazi1)
192
193			sinazi = cosU1 * sinazi1 # azi = azimuth of the geodesic at the equator
194			cosSqazi = 1 - sinazi * sinazi
195			uSq = cosSqazi * (a * a - b * b) / (b * b)
196			A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
197			B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
198
199			sigma = s / (b * A)
200			# sigma = angular distance P₁ P₂ on the sphere
201			# sigmam = angular distance on the sphere from the equator to
202			# the midpoint of the line
203
204			iterations = 0
205			sigmaprime = 0
206
207			while (np.abs(sigma - sigmaprime) > 1e-12).any():
208			cos2sigmam = np.cos(2 * sigma1 + sigma)
209			sinsigma = np.sin(sigma)
210			cossigma = np.cos(sigma)
211			deltasigma = (
212			B
213			* sinsigma
214			* (
215			cos2sigmam
216			+ B
217			/ 4
218			* (
219			cossigma * (-1 + 2 * cos2sigmam * cos2sigmam)
220			- B
221			/ 6
222			* cos2sigmam
223			* (-3 + 4 * sinsigma * sinsigma)
224			* (-3 + 4 * cos2sigmam * cos2sigmam)
225			)
226			)
227			)
228			sigmaprime = sigma
229			sigma = s / (b * A) + deltasigma
230			iterations += 1
231			if iterations >= 100:
232			raise Exception('Vincenty formula failed to converge')
233
234			x = sinU1 * sinsigma - cosU1 * cossigma * cosazi1
			0 ignored issues – show introduced 2020-01-04 04:49 UTC by Report Bug Copy Issue Report The variable `sinsigma` does not seem to be defined in case the `while` loop on line `207` is not entered. Are you sure this can never be the case? Loading history... introduced 2020-01-04 04:49 UTC by Report Bug Copy Issue Report The variable `cossigma` does not seem to be defined in case the `while` loop on line `207` is not entered. Are you sure this can never be the case? Loading history...
235			lat2 = np.arctan2(
236			sinU1 * cossigma + cosU1 * sinsigma * cosazi1,
237			(1 - f) * np.sqrt(sinazi * sinazi + x * x),
238			)
239			lambd = np.arctan2(
240			sinsigma * sinazi1, cosU1 * cossigma - sinU1 * sinsigma * cosazi1
241			)
242			C = f / 16 * cosSqazi * (4 + f * (4 - 3 * cosSqazi))
243			L = lambd - (1 - C) * f * sinazi * (
244			sigma
245			+ C
246			* sinsigma
247			* (cos2sigmam + C * cossigma * (-1 + 2 * cos2sigmam * cos2sigmam))
			0 ignored issues – show introduced 2020-01-04 04:49 UTC by Report Bug Copy Issue Report The variable `cos2sigmam` does not seem to be defined in case the `while` loop on line `207` is not entered. Are you sure this can never be the case? Loading history...
248			)
249			lon2 = lon1 + L
250			azi2 = np.arctan2(sinazi, -x)
251
252			lat2 = wrap90deg(lat2 * 180.0 / np.pi)
253			lon2 = wrap180deg(lon2 * 180.0 / np.pi)
254			azi2 = wrap360deg(azi2 * 180.0 / np.pi)
255
256			return {'lat2': lat2, 'lon2': lon2, 'azi2': azi2, 'iterations': iterations}
257

omdv / geovectors

Push — master ( 9b4d36...ed671d )

geovectorslib.geod.direct() B

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like