nvector package

Submodules

nvector.info module

Introduction to Nvector

Nvector is a suite of tools written in Python to solve geographical position calculations like:

  • Calculate the surface distance between two geographical positions:*
  • Convert positions given in one reference frame into another reference frame*
  • Find the destination point given start position, azimuth/bearing and distance*
  • Find the mean position (center/midpoint) of several geographical positions.*
  • Find the intersection between two paths.*
  • Find the cross track distance between a path and a position.*

In this library, we represent position with an “n-vector”, which is the normal vector to the Earth model (the same reference ellipsoid that is used for latitude and longitude). When using n-vector, all Earth-positions are treated equally, and there is no need to worry about singularities or discontinuities. An additional benefit with using n-vector is that many position calculations can be solved with simple vector algebra (e.g. dot product and cross product).

Converting between n-vector and latitude/longitude is unambiguous and easy using the provided functions.

n_E is n-vector in the program code, while in documents we use nE. E denotes an Earth-fixed coordinate frame, and it indicates that the three components of n-vector are along the three axes of E. More details about the notation used are found here: http://www.navlab.net/nvector/

The core functions provided are:

lat_lon2n_E:
Converts latitude and longitude to n-vector.
n_E2lat_lon:
Converts n-vector to latitude and longitude.
n_EB_E2p_EB_E:
Converts n-vector to Cartesian position vector in meters.
p_EB_E2n_EB_E:
Converts Cartesian position vector in meters to n-vector.
n_EA_E_and_n_EB_E2p_AB_E:
From two positions A and B, finds the delta position.

Nvector also provide an object oriented interface.

FrameE:
frame of reference rotates and moves with the Earth. Origo = Earth’s centre. z-axis->North, x-axis->Latitude=Longitude=0
FrameB:
frame of reference rotates and moves with Body. Origo = Body’s centre. x-axis -> forward, y-axis -> starboard, z-axis -> down
FrameN:
frame of reference moves with Body and rotates with Earth. Origo = Beneath/above Body at Earth’s surface. x-axis -> North, y-axis -> East, z-axis -> down
FrameL:
frame of reference moves with Body, but does not rotate with Earth. Origo = Beneath/above Body at Earth’s surface.
ECEFvector:
Geographical position given as Cartesian position vector in frame E
GeoPoint:
Geographical position given as latitude, longitude, depth in frame E
Nvector:
Geographical position given as N-vector and depth in frame E
GeoPath:
Geodesic path between two points in Frame E

Documentation is at: http://www.navlab.net/nvector/

Code and issue tracker is at https://github.com/pbrod/nvector.

Latest stable release is at http://pypi.python.org/pypi/Nvector.

To test if the toolbox is working paste the following in an interactive python session:

import nvector as nv
nv.test(coverage=True, doctests=True)

Getting Started

Example 1: “A and B to delta”

Given two positions, A and B as latitudes, longitudes and depths relative to Earth, E.

Find the exact vector between the two positions, given in meters north, east, and down, and find the direction (azimuth) to B, relative to north. Assume WGS-84 ellipsoid. The given depths are from the ellipsoid surface. Use position A to define north, east, and down directions. (Due to the curvature of Earth and different directions to the North Pole, the north, east, and down directions will change (relative to Earth) for different places. A must be outside the poles for the north and east directions to be defined.)

Solution:
>>> import numpy as np
>>> import nvector as nv
>>> wgs84 = nv.FrameE(name='WGS84')
>>> pointA = wgs84.GeoPoint(latitude=1, longitude=2, z=3, degrees=True)
>>> pointB = wgs84.GeoPoint(latitude=4, longitude=5, z=6, degrees=True)
Step 1: Convert to ECEF vectors
>>> p_EA_E = pointA.to_ecef_vector()
>>> p_EB_E = pointB.to_ecef_vector()
Step 2: Find p_AB_E (delta decomposed in E).
>>> p_AB_E = p_EB_E - p_EA_E  # (delta decomposed in E).
Step 3: Find p_AB_N (delta decomposed in N).
>>> frame_N = nv.FrameN(pointA)
>>> p_AB_N = p_AB_E.change_frame(frame_N)
>>> p_AB_N = p_AB_N.pvector.ravel()
>>> valtxt = '{0:8.2f}, {1:8.2f}, {2:8.2f}'.format(*p_AB_N)
>>> 'delta north, east, down = {}'.format(valtxt)
'delta north, east, down = 331730.23, 332997.87, 17404.27'
Step4: Also find the direction (azimuth) to B, relative to north:
>>> azimuth = np.arctan2(p_AB_N[1], p_AB_N[0])
>>> 'azimuth = {0:4.2f} deg'.format(np.rad2deg(azimuth))
'azimuth = 45.11 deg'

Example 2: “B and delta to C”

A radar or sonar attached to a vehicle B (Body coordinate frame) measures the distance and direction to an object C. We assume that the distance and two angles (typically bearing and elevation relative to B) are already combined to the vector p_BC_B (i.e. the vector from B to C, decomposed in B). The position of B is given as n_EB_E and z_EB, and the orientation (attitude) of B is given as R_NB (this rotation matrix can be found from roll/pitch/yaw by using zyx2R).

Find the exact position of object C as n-vector and depth ( n_EC_E and z_EC ), assuming Earth ellipsoid with semi-major axis a and flattening f. For WGS-72, use a = 6 378 135 m and f = 1/298.26.

Solution:
>>> import nvector as nv
>>> wgs72 = nv.FrameE(name='WGS72')
>>> wgs72 = nv.FrameE(a=6378135, f=1.0/298.26)
Step 1: Position and orientation of B is given 400m above E:
>>> n_EB_E = wgs72.Nvector(nv.unit([[1], [2], [3]]), z=-400)
Step 2: Delta BC decomposed in B
>>> frame_B = nv.FrameB(n_EB_E, yaw=10, pitch=20, roll=30, degrees=True)
>>> p_BC_B = frame_B.Pvector(np.r_[3000, 2000, 100].reshape((-1, 1)))
Step 3: Decompose delta BC in E
>>> p_BC_E = p_BC_B.to_ecef_vector()
Step 4: Find point C by adding delta BC to EB
>>> p_EB_E = n_EB_E.to_ecef_vector()
>>> p_EC_E = p_EB_E + p_BC_E
>>> pointC = p_EC_E.to_geo_point()
>>> lat, lon, z = pointC.latitude_deg, pointC.longitude_deg, pointC.z
>>> msg = 'Pos C: lat, lon = {:4.2f}, {:4.2f} deg,  height = {:4.2f} m'
>>> msg.format(lat[0], lon[0], -z[0])
'Pos C: lat, lon = 53.33, 63.47 deg,  height = 406.01 m'

Example 3: “ECEF-vector to geodetic latitude”

Position B is given as an “ECEF-vector” p_EB_E (i.e. a vector from E, the center of the Earth, to B, decomposed in E). Find the geodetic latitude, longitude and height (latEB, lonEB and hEB), assuming WGS-84 ellipsoid.

Solution:
>>> import nvector as nv
>>> wgs84 = nv.FrameE(name='WGS84')
>>> position_B = 6371e3 * np.vstack((0.9, -1, 1.1))  # m
>>> p_EB_E = wgs84.ECEFvector(position_B)
Step 1: Find position B as geodetic latitude, longitude and height
>>> pointB = p_EB_E.to_geo_point()
Step 2: Extract latitude and longitude in degrees
>>> lat, lon, h = pointB.latitude_deg, pointB.longitude_deg, -pointB.z
>>> msg = 'Pos B: lat, lon = {:4.2f}, {:4.2f} deg, height = {:9.2f} m'
>>> msg.format(lat[0], lon[0], h[0])
'Pos B: lat, lon = 39.38, -48.01 deg, height = 4702059.83 m'

Example 4: “Geodetic latitude to ECEF-vector”

>>> wgs84 = nv.FrameE(name='WGS84')
>>> pointB = wgs84.GeoPoint(latitude=1, longitude=2, z=-3, degrees=True)
>>> p_EB_E = pointB.to_ecef_vector()
>>> 'Ex4: p_EB_E = {} m'.format(p_EB_E.pvector.ravel())
'Ex4: p_EB_E = [ 6373290.27721828   222560.20067474   110568.82718179] m'

Example 5: “Surface distance”

Find the surface distance sAB (i.e. great circle distance) between two positions A and B. The heights of A and B are ignored, i.e. if they don’t have zero height, we seek the distance between the points that are at the surface of the Earth, directly above/below A and B. The Euclidean distance (chord length) dAB should also be found. Use Earth radius 6371e3 m.

Solution:
>>> frame_E = nv.FrameE(a=6371e3, f=0)
>>> positionA = frame_E.GeoPoint(latitude=88, longitude=0, degrees=True)
>>> positionB = frame_E.GeoPoint(latitude=89, longitude=-170, degrees=True)
>>> s_AB, _azia, _azib = positionA.distance_and_azimuth(positionB)
>>> p_AB_E = positionB.to_ecef_vector() - positionA.to_ecef_vector()
>>> d_AB = np.linalg.norm(p_AB_E.pvector, axis=0)[0]
>>> msg = 'Great circle and Euclidean distance = {:5.2f} km, {:5.2f} km'
>>> msg.format(s_AB / 1000, d_AB / 1000)
'Great circle and Euclidean distance = 332.46 km, 332.42 km'
Alternative solution:
>>> path = nv.GeoPath(positionA, positionB)
>>> s_AB2 = path.track_distance(method='greatcircle').ravel()
>>> d_AB2 = path.track_distance(method='euclidean').ravel()
>>> msg.format(s_AB2[0] / 1000, d_AB2[0] / 1000)
'Great circle and Euclidean distance = 332.46 km, 332.42 km'
Exact solution for the WGS84 ellipsoid:
>>> wgs84 = nv.FrameE(name='WGS84')
>>> point1 = wgs84.GeoPoint(latitude=88, longitude=0, degrees=True)
>>> point2 = wgs84.GeoPoint(latitude=89, longitude=-170, degrees=True)
>>> s_12, _azi1, _azi2 = point1.distance_and_azimuth(point2)
>>> p_12_E = point2.to_ecef_vector() - point1.to_ecef_vector()
>>> d_12 = np.linalg.norm(p_12_E.pvector, axis=0)[0]
>>> msg.format(s_12 / 1000, d_12 / 1000)
'Great circle and Euclidean distance = 333.95 km, 333.91 km'

Example 7: “Mean position”

Three positions A, B, and C are given as n-vectors n_EA_E, n_EB_E, and n_EC_E. Find the mean position, M, given as n_EM_E. Note that the calculation is independent of the depths of the positions.

Solution:
>>> points = nv.GeoPoint(latitude=[90, 60, 50],
...                      longitude=[0, 10, -20], degrees=True)
>>> nvectors = points.to_nvector()
>>> n_EM_E = nvectors.mean_horizontal_position()
>>> g_EM_E = n_EM_E.to_geo_point()
>>> lat, lon = g_EM_E.latitude_deg, g_EM_E.longitude_deg
>>> msg = 'Pos M: lat, lon = {:4.2f}, {:4.2f} deg'
>>> msg.format(lat[0], lon[0])
'Pos M: lat, lon = 67.24, -6.92 deg'

Example 8: “A and azimuth/distance to B”

We have an initial position A, direction of travel given as an azimuth (bearing) relative to north (clockwise), and finally the distance to travel along a great circle given as sAB. Use Earth radius 6371e3 m to find the destination point B.

In geodesy this is known as “The first geodetic problem” or “The direct geodetic problem” for a sphere, and we see that this is similar to Example 2, but now the delta is given as an azimuth and a great circle distance. (“The second/inverse geodetic problem” for a sphere is already solved in Examples 1 and 5.)

Solution:
>>> frame = nv.FrameE(a=6371e3, f=0)
>>> pointA = frame.GeoPoint(latitude=80, longitude=-90, degrees=True)
>>> pointB, _azimuthb = pointA.geo_point(distance=1000, azimuth=200,
...                                      degrees=True)
>>> latB, lonB = pointB.latitude_deg, pointB.longitude_deg
>>> 'Ex8, Destination: lat, lon = {:4.2f}, {:4.2f} deg'.format(latB, lonB)
'Ex8, Destination: lat, lon = 79.99, -90.02 deg'

Example 9: “Intersection of two paths”

Define a path from two given positions (at the surface of a spherical Earth), as the great circle that goes through the two points.

Path A is given by A1 and A2, while path B is given by B1 and B2.

Find the position C where the two paths intersect.

Solution 9:

>>> pointA1 = nv.GeoPoint(10, 20, degrees=True)
>>> pointA2 = nv.GeoPoint(30, 40, degrees=True)
>>> pointB1 = nv.GeoPoint(50, 60, degrees=True)
>>> pointB2 = nv.GeoPoint(70, 80, degrees=True)
>>> pathA = nv.GeoPath(pointA1, pointA2)
>>> pathB = nv.GeoPath(pointB1, pointB2)
>>> pointC = pathA.intersection(pathB)
>>> lat, lon = pointC.latitude_deg, pointC.longitude_deg
>>> msg = 'Ex9, Intersection: lat, long = {:4.2f}, {:4.2f} deg'
>>> msg.format(lat[0], lon[0])
'Ex9, Intersection: lat, long = 40.32, 55.90 deg'

Example 10: “Cross track distance”

Path A is given by the two positions A1 and A2 (similar to the previous example).

Find the cross track distance sxt between the path A (i.e. the great circle through A1 and A2) and the position B (i.e. the shortest distance at the surface, between the great circle and B).

Also find the Euclidean distance dxt between B and the plane defined by the great circle. Use Earth radius 6371e3.

Solution 10:
>>> frame = nv.FrameE(a=6371e3, f=0)
>>> pointA1 = frame.GeoPoint(0, 0, degrees=True)
>>> pointA2 = frame.GeoPoint(10, 0, degrees=True)
>>> pointB = frame.GeoPoint(1, 0.1, degrees=True)
>>> pathA = nv.GeoPath(pointA1, pointA2)
>>> s_xt = pathA.cross_track_distance(pointB, method='greatcircle').ravel()
>>> d_xt = pathA.cross_track_distance(pointB, method='euclidean').ravel()
>>> val_txt = '{:4.2f} km, {:4.2f} km'.format(s_xt[0]/1000, d_xt[0]/1000)
>>> msg = 'cross track distance from path A to position B'
>>> '{}, s_xt, d_xt = {}'.format(msg, val_txt)
'cross track distance from path A to position B, s_xt, d_xt = 11.12 km, 11.12 km'

Below we also give the functional solutions to example 1.

Example 1: Find the exact vector between the two positions, given in meters
north, east, and down, i.e. find p_AB_N:
>>> import numpy as np
>>> import nvector as nv
>>> from nvector import rad, deg
>>> lat_EA, lon_EA, z_EA = rad(1), rad(2), 3
>>> lat_EB, lon_EB, z_EB = rad(4), rad(5), 6

SOLUTION: Step1: Convert to n-vectors:

>>> n_EA_E = nv.lat_lon2n_E(lat_EA, lon_EA)
>>> n_EB_E = nv.lat_lon2n_E(lat_EB, lon_EB)

Step2: Find p_AB_E (delta decomposed in E). WGS-84 ellipsoid is default:

>>> p_AB_E = nv.n_EA_E_and_n_EB_E2p_AB_E(n_EA_E, n_EB_E, z_EA, z_EB)
Step3: Find R_EN for position A:
>>> R_EN = nv.n_E2R_EN(n_EA_E)
Step4: Find p_AB_N (delta decomposed in N).
>>> p_AB_N = np.dot(R_EN.T, p_AB_E).ravel()
>>> 'delta north, east, down = {0:8.2f}, {1:8.2f}, {2:8.2f}'.format(*p_AB_N)
'delta north, east, down = 331730.23, 332997.87, 17404.27'
Step5: Also find the direction (azimuth) to B, relative to north:
>>> azimuth = np.arctan2(p_AB_N[1], p_AB_N[0]) # positive angle about down-axis
>>> 'azimuth = {0:4.2f} deg'.format(deg(azimuth))
'azimuth = 45.11 deg'

See also

geographiclib

nvector.info.test_docstrings()[source]

nvector.navigator module

Module providing function and classes for calculating distances and bearing from ship to sensor as well as the position of the sensor relative to the ship.

class nvector.navigator.Navigator(source=None)[source]

Bases: nvector.navigator.Subject

Reads ship and sensor positions from file, and calculates distance and bearing from ship to sensor as well as the position of the sensor relative to the ship.

source : BufferSource

distance : real scalar
distance from gps to sensor
bearing : real scalar
angle formed by North-pole, gps-position and sensor-position
true_heading : real scalar
of ship
speed: real scalar
of ship in knots

relative_position :

handle_sample()[source]

Handles the next sentences and processes them

set_sample_time(sample_time)[source]
update_source(source)[source]

Updates source and initializes sensor location

class nvector.navigator.NmeaFileReader(file_path)[source]

Bases: object

file_path : string
file name with path
get_data()[source]

Return all data from the NMEA file

iget_data()[source]

Iterates over data from NMEA file

nvector.ngs module

Vincenty’s formulae are two related iterative methods used in geodesy to calculate the distance between two points on the surface of a spheroid, developed by Thaddeus Vincenty (1975a) They are based on the assumption that the figure of the Earth is an oblate spheroid, and hence are more accurate than methods such as great-circle distance which assume a spherical Earth. The first (direct) method computes the location of a point which is a given distance and azimuth (direction) from another point. The second (inverse) method computes the geographical distance and azimuth between two given points. They have been widely used in geodesy because they are accurate to within 0.5 mm (0.020’‘) on the Earth ellipsoid.

Reference

name: inverse version: 201105.xx author: stephen j. frakes last mod: dr. dennis milbert purpose: to compute a geodetic inverse

and then display output information

a semimajor axis equatorial (in meters) f flattening b semiminor axis polar (in meters) baz azimuth back (in radians)

dlon temporary value for difference in longitude (radians)

edist ellipsoid distance (in meters) elips ellipsoid choice

faz azimuth forward (in radians)

finv reciprocal flattening

option user prompt response

name1 name of station one glat1,glon1 station one - (lat & lon in radians )

name2 name of station two glat2,glon2 station two - (lat & lon in radians )

class nvector.ngs.Geodesic(a=6378137, f=0.0033528106647474805, name='')[source]

Bases: object

Solve geodesic problems.

The following illustrates its use

>>> import numpy as np
>>> from ngs import Geodesic
>>> wgs84 = Geodesic(name='WGS84')

# The geodesic inverse problem >>> lat1, lon1 = np.deg2rad((-41.32, 174.81)) >>> lat2, lon2 = np.deg2rad((40.96, -5.50)) >>> s12, az1, az2 = wgs84.inverse(lat1, lon1, lat2, lon2)[:3]

# The geodesic direct problem >>> lat1, lon1, az1 = np.deg2rad((40.6, -73.8, 45)) >>> lat2, lon2, az2 = wgs84.direct(lat1, lon1, az1, 10000e3)

All angles (latitudes, longitudes, azimuths, spherical arc lengths) are measured in radians. Latitudes must lie in [-pi/2,pi/2]. All lengths (distance, reduced length) are measured in meters.

direct(lat1, lon1, faz, S)[source]

SOLUTION OF THE GEODETIC DIRECT PROBLEM AFTER T.VINCENTY MODIFIED RAINSFORD’S METHOD WITH HELMERT’S ELLIPTICAL TERMS EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL

A IS THE SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID F IS THE FLATTENING OF THE REFERENCE ELLIPSOID LATITUDES AND LONGITUDES IN RADIANS POSITIVE NORTH AND EAST AZIMUTHS IN RADIANS CLOCKWISE FROM NORTH GEODESIC DISTANCE S ASSUMED IN UNITS OF SEMI-MAJOR AXIS A

PROGRAMMED FOR CDC-6600 BY LCDR L.PFEIFER NGS ROCKVILLE MD 20FEB75 MODIFIED FOR SYSTEM 360 BY JOHN G GERGEN NGS ROCKVILLE MD 750608

inverse(lat1, lon1, lat2, lon2)[source]

Return ellipsoidal distance between points inverse for long-line and antipodal cases. latitudes may be 90 degrees exactly. latitude positive north, longitude positive east, radians. azimuth clockwise from north, radians. original programmed by thaddeus vincenty, 1975, 1976 removed back side solution option, debugged, revised – 2011may01 – dgm this version of code is interim – antipodal boundary needs work

s12: real scalar
ellopsoidal distance between point 1 and 2
faz, baz: real scalars
forward and backward azimuth

sig, spherical distance on auxiliary sphere lam, longitude difference on auxiliary sphere kind, solution flag: kind=1, long-line; kind=2, antipodal it, iteration count

nvector.ngs.main()[source]
nvector.ngs.select_ellipsoid(name)[source]

nvector.objects module

Created on 29. des. 2015

@author: pab

class nvector.objects.FrameE(a=None, f=None, name='WGS84', north='z')[source]

Bases: nvector.objects._BaseFrame

Name:
Earth
Position:
The origin coincides with Earth’s centre (geometrical centre of ellipsoid model).
Orientation:
The x-axis is along the Earth’s rotation axis, pointing north (the yz-plane coincides with the equatorial plane), the y-axis points towards longitude +90x (east).
Comments:
The frame is Earth-fixed (rotates and moves with the Earth). The choice of axis directions ensures that at zero latitude and longitude, N (described below) has the same orientation as E. If roll/pitch/yaw are zero, also B (described below) has this orientation. Note that these properties are not valid for another common choice of the axis directions, denoted e (lower case), which has z pointing north and x pointing to latitude=longitude=0.
ECEFvector(*args, **kwds)[source]
GeoPoint(*args, **kwds)[source]
Nvector(*args, **kwds)[source]
direct(lat_a, lon_a, azimuth, distance, z=0, long_unroll=False, degrees=False)[source]
inverse(lat_a, lon_a, lat_b, lon_b, z=0, long_unroll=False, degrees=False)[source]
class nvector.objects.FrameB(position, yaw=0, pitch=0, roll=0, degrees=False)[source]

Bases: nvector.objects.FrameN

Name:
Body (typically of a vehicle)
Position:
The origin is in the vehicle’s reference point.
Orientation:
The x-axis points forward, the y-axis to the right (starboard) and the z-axis in the vehicle’s down direction.
Comments:
The frame is fixed to the vehicle.
R_EN
class nvector.objects.FrameL(position, wander_azimuth=0)[source]

Bases: nvector.objects.FrameN

Name:
Local level, Wander azimuth
Position:
The origin is directly beneath or above the vehicle (B), at Earth’s surface (surface of ellipsoid model).
Orientation:
The z-axis is pointing down. Initially, the x-axis points towards north, and the y-axis points towards east, but as the vehicle moves they are not rotating about the z-axis (their angular velocity relative to the Earth has zero component along the z-axis). (Note: Any initial horizontal direction of the x- and y-axes is valid for L, but if the initial position is outside the poles, north and east are usually chosen for convenience.)
Comments:
The L-frame is equal to the N-frame except for the rotation about the z-axis, which is always zero for this frame (relative to E). Hence, at a given time, the only difference between the frames is an angle between the x-axis of L and the north direction; this angle is called the wander azimuth angle. The L-frame is well suited for general calculations, as it is non-singular.
class nvector.objects.FrameN(position)[source]

Bases: nvector.objects._BaseFrame

Name:
North-East-Down (local level)
Position:
The origin is directly beneath or above the vehicle (B), at Earth’s surface (surface of ellipsoid model).
Orientation:
The x-axis points towards north, the y-axis points towards east (both are horizontal), and the z-axis is pointing down.
Comments:
When moving relative to the Earth, the frame rotates about its z-axis to allow the x-axis to always point towards north. When getting close to the poles this rotation rate will increase, being infinite at the poles. The poles are thus singularities and the direction of the x- and y-axes are not defined here. Hence, this coordinate frame is NOT SUITABLE for general calculations.
Pvector(pvector)[source]
class nvector.objects.GeoPoint(latitude, longitude, z=0, frame=None, degrees=False)[source]

Bases: object

Geographical position given as latitude, longitude, depth in frame E

latitude, longitude: real scalars or vectors of length n.
Geodetic latitude and longitude given in [rad or deg]
z: real scalar or vector of length n.
Depth(s) [m] relative to the ellipsoid (depth = -height)
frame: FrameE object
reference ellipsoid. The default ellipsoid model used is WGS84, but other ellipsoids/spheres might be specified.
degrees: bool
True if input are given in degrees otherwise radians are assumed.

Solve geodesic problems.

The following illustrates its use

>>> import nvector as nv
>>> wgs84 = nv.FrameE(name='WGS84')

# The geodesic inverse problem >>> point1 = wgs84.GeoPoint(-41.32, 174.81, degrees=True)) >>> point2 = wgs84.GeoPoint(40.96, -5.50, degrees=True) >>> s12, az1, az2 = point1.distance_and_azimuth(point2, degrees=True) >>> ‘s12 = {:5.2f}, az1 = {:5.2f}, az2 = {:5.2f}’.format(s12, az1, az2) ‘s12 = 19959679.27, az1 = 161.07, az2 = 18.83’

# The geodesic direct problem >>> point1 = wgs84.GeoPoint(40.6, -73.8, degrees=True) >>> az1, distance = 45, 10000e3 >>> point2, az2 = point1.geo_point(distance, az1, degrees=True) >>> lat2, lon2 = point2.latitude_deg, point2.longitude_deg >>> ‘lat2 = {:5.2f}, lon2 = {:5.2f}, az2 = {:5.2f}’.format(lat2, lon2, az2) ‘lat2 = 32.64, lon2 = 49.01, az2 = 140.37’

distance_and_azimuth(point, long_unroll=False, degrees=False)[source]

Return ellipsoidal distance between positions as well as the direction.

point: GeoPoint object
Latitude and longitude of position B.
degrees: bool
azimuths are returned in degrees if True otherwise in radians.
s_ab: real scalar
ellipsoidal distance [m] between position A and B.
azimuth_a, azimuth_b
direction [rad or deg] of line at position A and B relative to North, respectively.
geo_point(distance, azimuth, long_unroll=False, degrees=False)[source]

Return position B computed from current position, distance and azimuth.

distance: real scalar
ellipsoidal distance [m] between position A and B.
azimuth_a:
azimuth [rad or deg] of line at position A.
degrees: bool
azimuths are given in degrees if True otherwise in radians.
point_b: GeoPoint object
latitude and longitude of position B.
azimuth_b
azimuth [rad or deg] of line at position B.
latitude_deg
longitude_deg
ravel(degrees=False)[source]
to_ecef_vector()[source]
to_nvector()[source]

Converts latitude and longitude to n-vector.

latitude, longitude: real scalars or vectors of length n.
Geodetic latitude and longitude given in [rad]
n_E: 3 x n array
n-vector(s) [no unit] decomposed in E.

n_E2lat_lon.

class nvector.objects.GeoPath(point1, point2)[source]

Bases: object

Geographical path between two points in Frame E

cross_track_distance(point, method='greatcircle', radius=None)[source]

Return cross track distance from the path to a point.

point: GeoPoint, Nvector or ECEFvector object
position to measure the cross track distance to.
radius: real scalar
radius of sphere in [m]. Default mean Earth radius
method: string
defining distance calculated. Options are: ‘greatcircle’ or ‘euclidean’
distance: real scalar
distance in [m]
intersection(path)[source]

Return the intersection between the paths

path: GeoPath object
path to intersect
point: GeoPoint
point of intersection between paths
nvectors()[source]
track_distance(method='greatcircle', radius=None)[source]

Return the distance of the path.

class nvector.objects.Nvector(normal, z=0, frame=None)[source]

Bases: object

Geographical position given as N-vector and depth in frame E

normal: 3 x n array
n-vector(s) [no unit] decomposed in E.
z: real scalar or vector of length n.
Depth(s) [m] relative to the ellipsoid (depth = -height)
frame: FrameE object
reference ellipsoid. The default ellipsoid model used is WGS84, but other ellipsoids/spheres might be specified.

The position of B (typically body) relative to E (typically Earth) is given into this function as n-vector, n_EB_E and a depth, z relative to the ellipsiod.

GeoPoint, ECEFvector, Pvector

mean_horizontal_position()[source]

Return the n-vector of the horizontal mean position.

p_EM_E: 3 x 1 array
n-vector [no unit] of the mean position, decomposed in E.
to_ecef_vector()[source]

Converts n-vector to Cartesian position vector (“ECEF-vector”)

p_EB_E: ECEFvector object
Cartesian position vector(s) from E to B, decomposed in E.

The calculation is excact, taking the ellipsity of the Earth into account. It is also non-singular as both n-vector and p-vector are non-singular (except for the center of the Earth).

n_EB_E2p_EB_E, ECEFvector, Pvector, GeoPoint

to_geo_point()[source]

Converts n-vector to geo-point.

n_E2lat_lon, GeoPoint, ECEFvector, Pvector

to_nvector()[source]
class nvector.objects.Pvector(pvector, frame)[source]

Bases: object

to_ecef_vector()[source]
to_geo_point()[source]
to_nvector()[source]
class nvector.objects.ECEFvector(pvector, frame=None)[source]

Bases: object

Geographical position given as Cartesian position vector in frame E

pvector: 3 x n array
Cartesian position vector(s) [m] from E to B, decomposed in E.
frame: FrameE object
reference ellipsoid. The default ellipsoid model used is WGS84, but other ellipsoids/spheres might be specified.

The position of B (typically body) relative to E (typically Earth) is given into this function as p-vector, p_EB_E relative to the center of the frame.

GeoPoint, ECEFvector, Pvector

change_frame(frame)[source]

Converts to Cartesian position vector in another frame

frame: FrameB, FrameN or frameL object
Frame N used to convert p_AB_E (position vector from A to B, decomposed in E) to p_AB_N.
p_AB_N: Pvector object
position vector from A to B, decomposed in frame N.

n_EB_E2p_EB_E, n_EA_E_and_p_AB_E2n_EB_E, n_EA_E_and_n_EB_E2p_AB_E.

to_geo_point()[source]

Converts ECEF-vector to geo-point.

point: GeoPoint object
containing geodetic latitude and longitude given in [rad or deg] and depth, z, relative to the ellipsoid (depth = -height).

n_E2lat_lon, n_EB_E2p_EB_E, GeoPoint, Nvector, ECEFvector, Pvector

to_nvector()[source]

Converts Cartesian position vector to n-vector.

n_EB_E: Nvector object
n-vector(s) [no unit] of position B, decomposed in E.

The calculation is excact, taking the ellipsity of the Earth into account. It is also non-singular as both n-vector and p-vector are non-singular (except for the center of the Earth).

n_EB_E2p_EB_E, Nvector

nvector.plot module

Created on 9. des. 2015

@author: pab

nvector.plot.plot_mean_position()[source]
nvector.plot.plot_world(map1)[source]
map1: Basemap object
map1 to plot.

nvector.skeleton module

This is a skeleton file that can serve as a starting point for a Python console script. To run this script uncomment the following line in the entry_points section in setup.cfg:

console_scripts =
hello_world = nvector.module:function

Then run python setup.py install which will install the command hello_world inside your current environment. Besides console scripts, the header (i.e. until _logger...) of this file can also be used as template for Python modules.

Note: This skeleton file can be safely removed if not needed!

nvector.skeleton.main(args)[source]
nvector.skeleton.parse_args(args)[source]

Parse command line parameters

Parameters:args – command line parameters as list of strings
Returns:command line parameters as argparse.Namespace
nvector.skeleton.run()[source]

Module contents

Introduction to Nvector

Nvector is a suite of tools written in Python to solve geographical position calculations like:

  • Calculate the surface distance between two geographical positions:*
  • Convert positions given in one reference frame into another reference frame*
  • Find the destination point given start position, azimuth/bearing and distance*
  • Find the mean position (center/midpoint) of several geographical positions.*
  • Find the intersection between two paths.*
  • Find the cross track distance between a path and a position.*

In this library, we represent position with an “n-vector”, which is the normal vector to the Earth model (the same reference ellipsoid that is used for latitude and longitude). When using n-vector, all Earth-positions are treated equally, and there is no need to worry about singularities or discontinuities. An additional benefit with using n-vector is that many position calculations can be solved with simple vector algebra (e.g. dot product and cross product).

Converting between n-vector and latitude/longitude is unambiguous and easy using the provided functions.

n_E is n-vector in the program code, while in documents we use nE. E denotes an Earth-fixed coordinate frame, and it indicates that the three components of n-vector are along the three axes of E. More details about the notation used are found here: http://www.navlab.net/nvector/

The core functions provided are:

lat_lon2n_E:
Converts latitude and longitude to n-vector.
n_E2lat_lon:
Converts n-vector to latitude and longitude.
n_EB_E2p_EB_E:
Converts n-vector to Cartesian position vector in meters.
p_EB_E2n_EB_E:
Converts Cartesian position vector in meters to n-vector.
n_EA_E_and_n_EB_E2p_AB_E:
From two positions A and B, finds the delta position.

Nvector also provide an object oriented interface.

FrameE:
frame of reference rotates and moves with the Earth. Origo = Earth’s centre. z-axis->North, x-axis->Latitude=Longitude=0
FrameB:
frame of reference rotates and moves with Body. Origo = Body’s centre. x-axis -> forward, y-axis -> starboard, z-axis -> down
FrameN:
frame of reference moves with Body and rotates with Earth. Origo = Beneath/above Body at Earth’s surface. x-axis -> North, y-axis -> East, z-axis -> down
FrameL:
frame of reference moves with Body, but does not rotate with Earth. Origo = Beneath/above Body at Earth’s surface.
ECEFvector:
Geographical position given as Cartesian position vector in frame E
GeoPoint:
Geographical position given as latitude, longitude, depth in frame E
Nvector:
Geographical position given as N-vector and depth in frame E
GeoPath:
Geodesic path between two points in Frame E

Documentation is at: http://www.navlab.net/nvector/

Code and issue tracker is at https://github.com/pbrod/nvector.

Latest stable release is at http://pypi.python.org/pypi/Nvector.

To test if the toolbox is working paste the following in an interactive python session:

import nvector as nv
nv.test(coverage=True, doctests=True)

Getting Started

Example 1: “A and B to delta”

Given two positions, A and B as latitudes, longitudes and depths relative to Earth, E.

Find the exact vector between the two positions, given in meters north, east, and down, and find the direction (azimuth) to B, relative to north. Assume WGS-84 ellipsoid. The given depths are from the ellipsoid surface. Use position A to define north, east, and down directions. (Due to the curvature of Earth and different directions to the North Pole, the north, east, and down directions will change (relative to Earth) for different places. A must be outside the poles for the north and east directions to be defined.)

Solution:
>>> import numpy as np
>>> import nvector as nv
>>> wgs84 = nv.FrameE(name='WGS84')
>>> pointA = wgs84.GeoPoint(latitude=1, longitude=2, z=3, degrees=True)
>>> pointB = wgs84.GeoPoint(latitude=4, longitude=5, z=6, degrees=True)
Step 1: Convert to ECEF vectors
>>> p_EA_E = pointA.to_ecef_vector()
>>> p_EB_E = pointB.to_ecef_vector()
Step 2: Find p_AB_E (delta decomposed in E).
>>> p_AB_E = p_EB_E - p_EA_E  # (delta decomposed in E).
Step 3: Find p_AB_N (delta decomposed in N).
>>> frame_N = nv.FrameN(pointA)
>>> p_AB_N = p_AB_E.change_frame(frame_N)
>>> p_AB_N = p_AB_N.pvector.ravel()
>>> valtxt = '{0:8.2f}, {1:8.2f}, {2:8.2f}'.format(*p_AB_N)
>>> 'delta north, east, down = {}'.format(valtxt)
'delta north, east, down = 331730.23, 332997.87, 17404.27'
Step4: Also find the direction (azimuth) to B, relative to north:
>>> azimuth = np.arctan2(p_AB_N[1], p_AB_N[0])
>>> 'azimuth = {0:4.2f} deg'.format(np.rad2deg(azimuth))
'azimuth = 45.11 deg'

Example 2: “B and delta to C”

A radar or sonar attached to a vehicle B (Body coordinate frame) measures the distance and direction to an object C. We assume that the distance and two angles (typically bearing and elevation relative to B) are already combined to the vector p_BC_B (i.e. the vector from B to C, decomposed in B). The position of B is given as n_EB_E and z_EB, and the orientation (attitude) of B is given as R_NB (this rotation matrix can be found from roll/pitch/yaw by using zyx2R).

Find the exact position of object C as n-vector and depth ( n_EC_E and z_EC ), assuming Earth ellipsoid with semi-major axis a and flattening f. For WGS-72, use a = 6 378 135 m and f = 1/298.26.

Solution:
>>> import nvector as nv
>>> wgs72 = nv.FrameE(name='WGS72')
>>> wgs72 = nv.FrameE(a=6378135, f=1.0/298.26)
Step 1: Position and orientation of B is given 400m above E:
>>> n_EB_E = wgs72.Nvector(nv.unit([[1], [2], [3]]), z=-400)
Step 2: Delta BC decomposed in B
>>> frame_B = nv.FrameB(n_EB_E, yaw=10, pitch=20, roll=30, degrees=True)
>>> p_BC_B = frame_B.Pvector(np.r_[3000, 2000, 100].reshape((-1, 1)))
Step 3: Decompose delta BC in E
>>> p_BC_E = p_BC_B.to_ecef_vector()
Step 4: Find point C by adding delta BC to EB
>>> p_EB_E = n_EB_E.to_ecef_vector()
>>> p_EC_E = p_EB_E + p_BC_E
>>> pointC = p_EC_E.to_geo_point()
>>> lat, lon, z = pointC.latitude_deg, pointC.longitude_deg, pointC.z
>>> msg = 'Pos C: lat, lon = {:4.2f}, {:4.2f} deg,  height = {:4.2f} m'
>>> msg.format(lat[0], lon[0], -z[0])
'Pos C: lat, lon = 53.33, 63.47 deg,  height = 406.01 m'

Example 3: “ECEF-vector to geodetic latitude”

Position B is given as an “ECEF-vector” p_EB_E (i.e. a vector from E, the center of the Earth, to B, decomposed in E). Find the geodetic latitude, longitude and height (latEB, lonEB and hEB), assuming WGS-84 ellipsoid.

Solution:
>>> import nvector as nv
>>> wgs84 = nv.FrameE(name='WGS84')
>>> position_B = 6371e3 * np.vstack((0.9, -1, 1.1))  # m
>>> p_EB_E = wgs84.ECEFvector(position_B)
Step 1: Find position B as geodetic latitude, longitude and height
>>> pointB = p_EB_E.to_geo_point()
Step 2: Extract latitude and longitude in degrees
>>> lat, lon, h = pointB.latitude_deg, pointB.longitude_deg, -pointB.z
>>> msg = 'Pos B: lat, lon = {:4.2f}, {:4.2f} deg, height = {:9.2f} m'
>>> msg.format(lat[0], lon[0], h[0])
'Pos B: lat, lon = 39.38, -48.01 deg, height = 4702059.83 m'

Example 4: “Geodetic latitude to ECEF-vector”

>>> wgs84 = nv.FrameE(name='WGS84')
>>> pointB = wgs84.GeoPoint(latitude=1, longitude=2, z=-3, degrees=True)
>>> p_EB_E = pointB.to_ecef_vector()
>>> 'Ex4: p_EB_E = {} m'.format(p_EB_E.pvector.ravel())
'Ex4: p_EB_E = [ 6373290.27721828   222560.20067474   110568.82718179] m'

Example 5: “Surface distance”

Find the surface distance sAB (i.e. great circle distance) between two positions A and B. The heights of A and B are ignored, i.e. if they don’t have zero height, we seek the distance between the points that are at the surface of the Earth, directly above/below A and B. The Euclidean distance (chord length) dAB should also be found. Use Earth radius 6371e3 m.

Solution:
>>> frame_E = nv.FrameE(a=6371e3, f=0)
>>> positionA = frame_E.GeoPoint(latitude=88, longitude=0, degrees=True)
>>> positionB = frame_E.GeoPoint(latitude=89, longitude=-170, degrees=True)
>>> s_AB, _azia, _azib = positionA.distance_and_azimuth(positionB)
>>> p_AB_E = positionB.to_ecef_vector() - positionA.to_ecef_vector()
>>> d_AB = np.linalg.norm(p_AB_E.pvector, axis=0)[0]
>>> msg = 'Great circle and Euclidean distance = {:5.2f} km, {:5.2f} km'
>>> msg.format(s_AB / 1000, d_AB / 1000)
'Great circle and Euclidean distance = 332.46 km, 332.42 km'
Alternative solution:
>>> path = nv.GeoPath(positionA, positionB)
>>> s_AB2 = path.track_distance(method='greatcircle').ravel()
>>> d_AB2 = path.track_distance(method='euclidean').ravel()
>>> msg.format(s_AB2[0] / 1000, d_AB2[0] / 1000)
'Great circle and Euclidean distance = 332.46 km, 332.42 km'
Exact solution for the WGS84 ellipsoid:
>>> wgs84 = nv.FrameE(name='WGS84')
>>> point1 = wgs84.GeoPoint(latitude=88, longitude=0, degrees=True)
>>> point2 = wgs84.GeoPoint(latitude=89, longitude=-170, degrees=True)
>>> s_12, _azi1, _azi2 = point1.distance_and_azimuth(point2)
>>> p_12_E = point2.to_ecef_vector() - point1.to_ecef_vector()
>>> d_12 = np.linalg.norm(p_12_E.pvector, axis=0)[0]
>>> msg.format(s_12 / 1000, d_12 / 1000)
'Great circle and Euclidean distance = 333.95 km, 333.91 km'

Example 7: “Mean position”

Three positions A, B, and C are given as n-vectors n_EA_E, n_EB_E, and n_EC_E. Find the mean position, M, given as n_EM_E. Note that the calculation is independent of the depths of the positions.

Solution:
>>> points = nv.GeoPoint(latitude=[90, 60, 50],
...                      longitude=[0, 10, -20], degrees=True)
>>> nvectors = points.to_nvector()
>>> n_EM_E = nvectors.mean_horizontal_position()
>>> g_EM_E = n_EM_E.to_geo_point()
>>> lat, lon = g_EM_E.latitude_deg, g_EM_E.longitude_deg
>>> msg = 'Pos M: lat, lon = {:4.2f}, {:4.2f} deg'
>>> msg.format(lat[0], lon[0])
'Pos M: lat, lon = 67.24, -6.92 deg'

Example 8: “A and azimuth/distance to B”

We have an initial position A, direction of travel given as an azimuth (bearing) relative to north (clockwise), and finally the distance to travel along a great circle given as sAB. Use Earth radius 6371e3 m to find the destination point B.

In geodesy this is known as “The first geodetic problem” or “The direct geodetic problem” for a sphere, and we see that this is similar to Example 2, but now the delta is given as an azimuth and a great circle distance. (“The second/inverse geodetic problem” for a sphere is already solved in Examples 1 and 5.)

Solution:
>>> frame = nv.FrameE(a=6371e3, f=0)
>>> pointA = frame.GeoPoint(latitude=80, longitude=-90, degrees=True)
>>> pointB, _azimuthb = pointA.geo_point(distance=1000, azimuth=200,
...                                      degrees=True)
>>> latB, lonB = pointB.latitude_deg, pointB.longitude_deg
>>> 'Ex8, Destination: lat, lon = {:4.2f}, {:4.2f} deg'.format(latB, lonB)
'Ex8, Destination: lat, lon = 79.99, -90.02 deg'

Example 9: “Intersection of two paths”

Define a path from two given positions (at the surface of a spherical Earth), as the great circle that goes through the two points.

Path A is given by A1 and A2, while path B is given by B1 and B2.

Find the position C where the two paths intersect.

Solution 9:

>>> pointA1 = nv.GeoPoint(10, 20, degrees=True)
>>> pointA2 = nv.GeoPoint(30, 40, degrees=True)
>>> pointB1 = nv.GeoPoint(50, 60, degrees=True)
>>> pointB2 = nv.GeoPoint(70, 80, degrees=True)
>>> pathA = nv.GeoPath(pointA1, pointA2)
>>> pathB = nv.GeoPath(pointB1, pointB2)
>>> pointC = pathA.intersection(pathB)
>>> lat, lon = pointC.latitude_deg, pointC.longitude_deg
>>> msg = 'Ex9, Intersection: lat, long = {:4.2f}, {:4.2f} deg'
>>> msg.format(lat[0], lon[0])
'Ex9, Intersection: lat, long = 40.32, 55.90 deg'

Example 10: “Cross track distance”

Path A is given by the two positions A1 and A2 (similar to the previous example).

Find the cross track distance sxt between the path A (i.e. the great circle through A1 and A2) and the position B (i.e. the shortest distance at the surface, between the great circle and B).

Also find the Euclidean distance dxt between B and the plane defined by the great circle. Use Earth radius 6371e3.

Solution 10:
>>> frame = nv.FrameE(a=6371e3, f=0)
>>> pointA1 = frame.GeoPoint(0, 0, degrees=True)
>>> pointA2 = frame.GeoPoint(10, 0, degrees=True)
>>> pointB = frame.GeoPoint(1, 0.1, degrees=True)
>>> pathA = nv.GeoPath(pointA1, pointA2)
>>> s_xt = pathA.cross_track_distance(pointB, method='greatcircle').ravel()
>>> d_xt = pathA.cross_track_distance(pointB, method='euclidean').ravel()
>>> val_txt = '{:4.2f} km, {:4.2f} km'.format(s_xt[0]/1000, d_xt[0]/1000)
>>> msg = 'cross track distance from path A to position B'
>>> '{}, s_xt, d_xt = {}'.format(msg, val_txt)
'cross track distance from path A to position B, s_xt, d_xt = 11.12 km, 11.12 km'

Below we also give the functional solutions to example 1.

Example 1: Find the exact vector between the two positions, given in meters
north, east, and down, i.e. find p_AB_N:
>>> import numpy as np
>>> import nvector as nv
>>> from nvector import rad, deg
>>> lat_EA, lon_EA, z_EA = rad(1), rad(2), 3
>>> lat_EB, lon_EB, z_EB = rad(4), rad(5), 6

SOLUTION: Step1: Convert to n-vectors:

>>> n_EA_E = nv.lat_lon2n_E(lat_EA, lon_EA)
>>> n_EB_E = nv.lat_lon2n_E(lat_EB, lon_EB)

Step2: Find p_AB_E (delta decomposed in E). WGS-84 ellipsoid is default:

>>> p_AB_E = nv.n_EA_E_and_n_EB_E2p_AB_E(n_EA_E, n_EB_E, z_EA, z_EB)
Step3: Find R_EN for position A:
>>> R_EN = nv.n_E2R_EN(n_EA_E)
Step4: Find p_AB_N (delta decomposed in N).
>>> p_AB_N = np.dot(R_EN.T, p_AB_E).ravel()
>>> 'delta north, east, down = {0:8.2f}, {1:8.2f}, {2:8.2f}'.format(*p_AB_N)
'delta north, east, down = 331730.23, 332997.87, 17404.27'
Step5: Also find the direction (azimuth) to B, relative to north:
>>> azimuth = np.arctan2(p_AB_N[1], p_AB_N[0]) # positive angle about down-axis
>>> 'azimuth = {0:4.2f} deg'.format(deg(azimuth))
'azimuth = 45.11 deg'

See also

geographiclib