`orthodrome`¶

exception Slow[source]¶

class Loc(lat, lon)[source]¶

Simple location representation.

Attrib lat: Latitude in [deg].
Attrib lon: Longitude in [deg].

clip(x, mi, ma)[source]¶

Clip values of an array.

Parameters

x (numpy.ndarray) – Continunous data to be clipped.
mi (float) – Clip minimum.
ma (float) – Clip maximum.

Returns

Clipped data.

Return type

numpy.ndarray

wrap(x, mi, ma)[source]¶

Wrapping continuous data to fundamental phase values.

$x_{\mathrm{wrapped}} = x_{\mathrm{cont},i} - \frac{ x_{\mathrm{cont},i} - r_{\mathrm{min}} } { r_{\mathrm{max}} - r_{\mathrm{min}}} \cdot ( r_{\mathrm{max}} - r_{\mathrm{min}}),\quad x_{\mathrm{wrapped}}\; \in \;[ r_{\mathrm{min}},\, r_{\mathrm{max}}].$

Parameters

x (numpy.ndarray) – Continunous data to be wrapped.
mi (float) – Minimum value of wrapped data.
ma (float) – Maximum value of wrapped data.

Returns

Wrapped data.

Return type

numpy.ndarray

cosdelta(*args)[source]¶

Cosine of the angular distance between two points a and b on a sphere.

This function (find implementation below) returns the cosine of the distance angle ‘delta’ between two points a and b, coordinates of which are expected to be given in geographical coordinates and in degrees. For numerical stability a maximum of 1.0 is enforced.

$A_{\mathrm{lat'}} = \frac{ \pi}{180} \cdot A_{lat}, \quad A_{\mathrm{lon'}} = \frac{ \pi}{180} \cdot A_{lon}, \quad B_{\mathrm{lat'}} = \frac{ \pi}{180} \cdot B_{lat}, \quad B_{\mathrm{lon'}} = \frac{ \pi}{180} \cdot B_{lon}\\[0.5cm] \cos(\Delta) = \min( 1.0, \quad \sin( A_{\mathrm{lat'}}) \sin( B_{\mathrm{lat'}} ) + \cos(A_{\mathrm{lat'}}) \cos( B_{\mathrm{lat'}} ) \cos( B_{\mathrm{lon'}} - A_{\mathrm{lon'}} )$

Parameters

a (pyrocko.orthodrome.Loc) – Location point A.
b (pyrocko.orthodrome.Loc) – Location point B.

Returns

Cosdelta.

Return type

float

cosdelta_numpy(a_lats, a_lons, b_lats, b_lons)[source]¶

Cosine of the angular distance between two points a and b on a sphere.

This function returns the cosines of the distance angles delta between two points a and b given as numpy.ndarray. The coordinates are expected to be given in geographical coordinates and in degrees. For numerical stability a maximum of 1.0 is enforced.

Please find the details of the implementation in the documentation of the function pyrocko.orthodrome.cosdelta() above.

Parameters

a_lats (numpy.ndarray) – Latitudes in [deg] point A.
a_lons (numpy.ndarray) – Longitudes in [deg] point A.
b_lats (numpy.ndarray) – Latitudes in [deg] point B.
b_lons (numpy.ndarray, (N)) – Longitudes in [deg] point B.

Returns

Cosdelta.

azimuth(*args)[source]¶

Azimuth calculation.

This function (find implementation below) returns azimuth … between points a and b, coordinates of which are expected to be given in geographical coordinates and in degrees.

$A_{\mathrm{lat'}} = \frac{ \pi}{180} \cdot A_{lat}, \quad A_{\mathrm{lon'}} = \frac{ \pi}{180} \cdot A_{lon}, \quad B_{\mathrm{lat'}} = \frac{ \pi}{180} \cdot B_{lat}, \quad B_{\mathrm{lon'}} = \frac{ \pi}{180} \cdot B_{lon}\\ \varphi_{\mathrm{azi},AB} = \frac{180}{\pi} \arctan \left[ \frac{ \cos( A_{\mathrm{lat'}}) \cos( B_{\mathrm{lat'}} ) \sin(B_{\mathrm{lon'}} - A_{\mathrm{lon'}} )} {\sin ( B_{\mathrm{lat'}} ) - \sin( A_{\mathrm{lat'}} cosdelta) } \right]$

Parameters

a (pyrocko.orthodrome.Loc) – Location point A.
b (pyrocko.orthodrome.Loc) – Location point B.

Returns

Azimuth in degree

azimuth_numpy(a_lats, a_lons, b_lats, b_lons, _cosdelta=None)[source]¶

Calculation of the azimuth (track angle) from a location A towards B.

This function returns azimuths (track angles) from locations A towards B given in numpy.ndarray. Coordinates are expected to be given in geographical coordinates and in degrees.

Please find the details of the implementation in the documentation of the function pyrocko.orthodrome.azimuth().

Parameters

a_lats (numpy.ndarray, (N)) – Latitudes in [deg] point A.
a_lons (numpy.ndarray, (N)) – Longitudes in [deg] point A.
b_lats (numpy.ndarray, (N)) – Latitudes in [deg] point B.
b_lons (numpy.ndarray, (N)) – Longitudes in [deg] point B.

Returns

Azimuths in [deg].

Return type

numpy.ndarray, (N)

azibazi(*args, **kwargs)[source]¶

Azimuth and backazimuth from location A towards B and back.

Returns: Azimuth in [deg] from A to B, back azimuth in [deg] from B to A.
Return type: tuple[float, float]

azibazi_numpy(a_lats, a_lons, b_lats, b_lons, implementation='c')[source]¶

Azimuth and backazimuth from location A towards B and back.

Arguments are given as numpy.ndarray.

Parameters

a_lats (numpy.ndarray) – Latitude(s) in [deg] of point A.
a_lons (numpy.ndarray) – Longitude(s) in [deg] of point A.
b_lats (numpy.ndarray) – Latitude(s) in [deg] of point B.
b_lons (numpy.ndarray) – Longitude(s) in [deg] of point B.

Returns

Azimuth(s) in [deg] from A to B, back azimuth(s) in [deg] from B to A.

Return type

numpy.ndarray, numpy.ndarray

azidist_numpy(*args)[source]¶

Calculation of the azimuth (track angle) and the distance from locations A towards B on a sphere.

The assisting functions used are pyrocko.orthodrome.cosdelta() and pyrocko.orthodrome.azimuth()

Parameters

a_lats (numpy.ndarray, (N)) – Latitudes in [deg] point A.
a_lons (numpy.ndarray, (N)) – Longitudes in [deg] point A.
b_lats (numpy.ndarray, (N)) – Latitudes in [deg] point B.
b_lons (numpy.ndarray, (N)) – Longitudes in [deg] point B.

Returns

Azimuths in [deg], distances in [deg].

Return type

numpy.ndarray, (2xN)

distance_accurate50m(*args, **kwargs)[source]¶

Accurate distance calculation based on a spheroid of rotation.

Function returns distance in meter between points A and B, coordinates of which must be given in geographical coordinates and in degrees. The returned distance should be accurate to 50 m using WGS84. Values for the Earth’s equator radius and the Earth’s oblateness (f_oblate) are defined in the pyrocko configuration file pyrocko.config.

From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:

Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell, Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1

$F = \frac{\pi}{180} \frac{(A_{lat} + B_{lat})}{2}, \quad G = \frac{\pi}{180} \frac{(A_{lat} - B_{lat})}{2}, \quad l = \frac{\pi}{180} \frac{(A_{lon} - B_{lon})}{2} \quad \\[0.5cm] S = \sin^2(G) \cdot \cos^2(l) + \cos^2(F) \cdot \sin^2(l), \quad \quad C = \cos^2(G) \cdot \cos^2(l) + \sin^2(F) \cdot \sin^2(l)$

$w = \arctan \left( \sqrt{ \frac{S}{C}} \right) , \quad r = \sqrt{\frac{S}{C} }$

The spherical-earth distance D between A and B, can be given with:

$D_{sphere} = 2w \cdot R_{equator}$

The oblateness of the Earth requires some correction with correction factors h1 and h2:

$h_1 = \frac{3r - 1}{2C}, \quad h_2 = \frac{3r +1 }{2S}\\[0.5cm] D = D_{\mathrm{sphere}} \cdot [ 1 + h_1 \,f_{\mathrm{oblate}} \cdot \sin^2(F) \cos^2(G) - h_2\, f_{\mathrm{oblate}} \cdot \cos^2(F) \sin^2(G)]$

Parameters

a (pyrocko.orthodrome.Loc) – Location point A.
b (pyrocko.orthodrome.Loc) – Location point B.

Returns

Distance in [m].

Return type

float

distance_accurate50m_numpy(a_lats, a_lons, b_lats, b_lons, implementation='c')[source]¶

Accurate distance calculation based on a spheroid of rotation.

Function returns distance in meter between points a and b, coordinates of which must be given in geographical coordinates and in degrees. The returned distance should be accurate to 50 m using WGS84. Values for the Earth’s equator radius and the Earth’s oblateness (f_oblate) are defined in the pyrocko configuration file pyrocko.config.

From wikipedia (http://de.wikipedia.org/wiki/Orthodrome), based on:

Meeus, J.: Astronomical Algorithms, S 85, Willmann-Bell, Richmond 2000 (2nd ed., 2nd printing), ISBN 0-943396-61-1

$F_i = \frac{\pi}{180} \frac{(a_{lat,i} + a_{lat,i})}{2}, \quad G_i = \frac{\pi}{180} \frac{(a_{lat,i} - b_{lat,i})}{2}, \quad l_i= \frac{\pi}{180} \frac{(a_{lon,i} - b_{lon,i})}{2} \\[0.5cm] S_i = \sin^2(G_i) \cdot \cos^2(l_i) + \cos^2(F_i) \cdot \sin^2(l_i), \quad \quad C_i = \cos^2(G_i) \cdot \cos^2(l_i) + \sin^2(F_i) \cdot \sin^2(l_i)$

$w_i = \arctan \left( \sqrt{\frac{S_i}{C_i}} \right), \quad r_i = \sqrt{\frac{S_i}{C_i} }$

The spherical-earth distance D between a and b, can be given with:

$D_{\mathrm{sphere},i} = 2w_i \cdot R_{\mathrm{equator}}$

The oblateness of the Earth requires some correction with correction factors h1 and h2:

$h_{1.i} = \frac{3r - 1}{2C_i}, \quad h_{2,i} = \frac{3r +1 }{2S_i}\\[0.5cm] D_{AB,i} = D_{\mathrm{sphere},i} \cdot [1 + h_{1,i} \,f_{\mathrm{oblate}} \cdot \sin^2(F_i) \cos^2(G_i) - h_{2,i}\, f_{\mathrm{oblate}} \cdot \cos^2(F_i) \sin^2(G_i)]$

Parameters

a_lats (numpy.ndarray, (N)) – Latitudes in [deg] point A.
a_lons (numpy.ndarray, (N)) – Longitudes in [deg] point A.
b_lats (numpy.ndarray, (N)) – Latitudes in [deg] point B.
b_lons (numpy.ndarray, (N)) – Longitudes in [deg] point B.

Returns

Distances in [m].

Return type

numpy.ndarray, (N)

ne_to_latlon(lat0, lon0, north_m, east_m)[source]¶

Transform local cartesian coordinates to latitude and longitude.

From east and north coordinates (x and y coordinate numpy.ndarray) relative to a reference differences in longitude and latitude are calculated, which are effectively changes in azimuth and distance, respectively:

$\text{distance change:}\; \Delta {\bf{a}} &= \sqrt{{\bf{y}}^2 + {\bf{x}}^2 }/ \mathrm{R_E}, \text{azimuth change:}\; \Delta \bf{\gamma} &= \arctan( \bf{x} / \bf{y}).$

The projection used preserves the azimuths of the input points.

Parameters

lat0 (float) – Latitude origin of the cartesian coordinate system in [deg].
lon0 (float) – Longitude origin of the cartesian coordinate system in [deg].
north_m (numpy.ndarray, (N)) – Northing distances from origin in [m].
east_m (numpy.ndarray, (N)) – Easting distances from origin in [m].

Returns

Array with latitudes and longitudes in [deg].

Return type

numpy.ndarray, (2xN)

azidist_to_latlon(lat0, lon0, azimuth_deg, distance_deg)[source]¶

Absolute latitudes and longitudes are calculated from relative changes.

Convenience wrapper to azidist_to_latlon_rad() with azimuth and distance given in degrees.

Parameters

lat0 (float) – Latitude origin of the cartesian coordinate system in [deg].
lon0 (float) – Longitude origin of the cartesian coordinate system in [deg].
azimuth_deg (numpy.ndarray, (N)) – Azimuth from origin in [deg].
distance_deg (numpy.ndarray, (N)) – Distances from origin in [deg].

Returns

Array with latitudes and longitudes in [deg].

Return type

numpy.ndarray, (2xN)

azidist_to_latlon_rad(lat0, lon0, azimuth_rad, distance_rad)[source]¶

Absolute latitudes and longitudes are calculated from relative changes.

For numerical stability a range between of -1.0 and 1.0 is enforced for c and alpha.

$\Delta {\bf a}_i \; \text{and} \; \Delta \gamma_i \; \text{are relative distances and azimuths from lat0 and lon0 for \textit{i} source points of a finite source.}$

$\mathrm{b} &= \frac{\pi}{2} -\frac{\pi}{180}\;\mathrm{lat_0}\\ {\bf c}_i &=\arccos[\; \cos(\Delta {\bf{a}}_i) \cos(\mathrm{b}) + |\Delta \gamma_i| \, \sin(\Delta {\bf a}_i) \sin(\mathrm{b})\; ] \\ \mathrm{lat}_i &= \frac{180}{\pi} \left(\frac{\pi}{2} - {\bf c}_i \right)$

$\alpha_i &= \arcsin \left[ \; \frac{ \sin(\Delta {\bf a}_i ) \sin(|\Delta \gamma_i|)}{\sin({\bf c}_i)}\; \right] \\ \alpha_i &= \begin{cases} \alpha_i, &\text{if} \; \cos(\Delta {\bf a}_i) - \cos(\mathrm{b}) \cos({\bf{c}}_i) > 0, \; \text{else} \\ \pi - \alpha_i, & \text{if} \; \alpha_i > 0,\; \text{else}\\ -\pi - \alpha_i, & \text{if} \; \alpha_i < 0. \end{cases} \\ \mathrm{lon}_i &= \mathrm{lon_0} + \frac{180}{\pi} \, \frac{\Delta \gamma_i }{|\Delta \gamma_i|} \cdot \alpha_i \text{, with $\alpha_i \in [-\pi,\pi]$}$

Parameters

lat0 (float) – Latitude origin of the cartesian coordinate system in [deg].
lon0 (float) – Longitude origin of the cartesian coordinate system in [deg].
distance_rad (numpy.ndarray, (N)) – Distances from origin in [rad].
azimuth_rad (numpy.ndarray, (N)) – Azimuth from origin in [rad].

Returns

Array with latitudes and longitudes in [deg].

Return type

numpy.ndarray, (2xN)

ne_to_latlon_alternative_method(lat0, lon0, north_m, east_m)[source]¶

Transform local cartesian coordinates to latitude and longitude.

Like pyrocko.orthodrome.ne_to_latlon(), but this method (implementation below), although it should be numerically more stable, suffers problems at points which are across the pole as seen from the cartesian origin.

$\text{distance change:}\; \Delta {{\bf a}_i} &= \sqrt{{\bf{y}}^2_i + {\bf{x}}^2_i }/ \mathrm{R_E},\\ \text{azimuth change:}\; \Delta {\bf \gamma}_i &= \arctan( {\bf x}_i {\bf y}_i). \\ \mathrm{b} &= \frac{\pi}{2} -\frac{\pi}{180} \;\mathrm{lat_0}\\$

${{\bf z}_1}_i &= \cos{\left( \frac{\Delta {\bf a}_i - \mathrm{b}}{2} \right)} \cos {\left( \frac{|\gamma_i|}{2} \right) }\\ {{\bf n}_1}_i &= \cos{\left( \frac{\Delta {\bf a}_i + \mathrm{b}}{2} \right)} \sin {\left( \frac{|\gamma_i|}{2} \right) }\\ {{\bf z}_2}_i &= \sin{\left( \frac{\Delta {\bf a}_i - \mathrm{b}}{2} \right)} \cos {\left( \frac{|\gamma_i|}{2} \right) }\\ {{\bf n}_2}_i &= \sin{\left( \frac{\Delta {\bf a}_i + \mathrm{b}}{2} \right)} \sin {\left( \frac{|\gamma_i|}{2} \right) }\\ {{\bf t}_1}_i &= \arctan{\left( \frac{{{\bf z}_1}_i} {{{\bf n}_1}_i} \right) }\\ {{\bf t}_2}_i &= \arctan{\left( \frac{{{\bf z}_2}_i} {{{\bf n}_2}_i} \right) } \\[0.5cm] c &= \begin{cases} 2 \cdot \arccos \left( {{\bf z}_1}_i / \sin({{\bf t}_1}_i) \right),\; \text{if } |\sin({{\bf t}_1}_i)| > |\sin({{\bf t}_2}_i)|,\; \text{else} \\ 2 \cdot \arcsin{\left( {{\bf z}_2}_i / \sin({{\bf t}_2}_i) \right)}. \end{cases}\\$

${\bf {lat}}_i &= \frac{180}{ \pi } \left( \frac{\pi}{2} - {\bf {c}}_i \right) \\ {\bf {lon}}_i &= {\bf {lon}}_0 + \frac{180}{ \pi } \frac{\gamma_i}{|\gamma_i|}, \text{ with}\; \gamma \in [-\pi,\pi]$

Parameters

lat0 (float) – Latitude origin of the cartesian coordinate system in [deg].
lon0 (float) – Longitude origin of the cartesian coordinate system in [deg].
north_m (numpy.ndarray, (N)) – Northing distances from origin in [m].
east_m (numpy.ndarray, (N)) – Easting distances from origin in [m].

Returns

Array with latitudes and longitudes in [deg].

Return type

numpy.ndarray, (2xN)

latlon_to_ne(*args)[source]¶

Relative cartesian coordinates with respect to a reference location.

For two locations, a reference location A and another location B, given in geographical coordinates in degrees, the corresponding cartesian coordinates are calculated. Assisting functions are pyrocko.orthodrome.azimuth() and pyrocko.orthodrome.distance_accurate50m().

$D_{AB} &= \mathrm{distance\_accurate50m(}A, B \mathrm{)}, \quad \varphi_{\mathrm{azi},AB} = \mathrm{azimuth(}A,B \mathrm{)}\\[0.3cm] n &= D_{AB} \cdot \cos( \frac{\pi }{180} \varphi_{\mathrm{azi},AB} )\\ e &= D_{AB} \cdot \sin( \frac{\pi }{180} \varphi_{\mathrm{azi},AB})$

Parameters

refloc (pyrocko.orthodrome.Loc) – Location reference point.
loc (pyrocko.orthodrome.Loc) – Location of interest.

Returns

Northing and easting from refloc to location in [m].

Return type

tuple[float, float]

latlon_to_ne_numpy(lat0, lon0, lat, lon)[source]¶

Relative cartesian coordinates with respect to a reference location.

For two locations, a reference location (lat0, lon0) and another location B, given in geographical coordinates in degrees, the corresponding cartesian coordinates are calculated. Assisting functions are azimuth() and distance_accurate50m().

Parameters

lat0 – Latitude of the reference location in [deg].
lon0 – Longitude of the reference location in [deg].
lat – Latitude of the absolute location in [deg].
lon – Longitude of the absolute location in [deg].

Returns

(n, e): relative north and east positions in [m].

Return type

numpy.ndarray, (2xN)

Implemented formulations:

$D_{AB} &= \mathrm{distance\_accurate50m(}A, B \mathrm{)}, \quad \varphi_{\mathrm{azi},AB} = \mathrm{azimuth(}A,B \mathrm{)}\\[0.3cm] n &= D_{AB} \cdot \cos( \frac{\pi }{180} \varphi_{ \mathrm{azi},AB} )\\ e &= D_{AB} \cdot \sin( \frac{\pi }{180} \varphi_{ \mathrm{azi},AB} )$

positive_region(region)[source]¶

Normalize parameterization of a rectangular geographical region.

Parameters: region (tuple of float) – (west, east, south, north) in [deg].
Returns: (west, east, south, north), where west <= east and where west and east are in the range [-180., 180.+360.].
Return type: tuple of float

points_in_region(p, region)[source]¶

Check what points are contained in a rectangular geographical region.

Parameters

p (numpy.ndarray (N, 2)) – (lat, lon) pairs in [deg].
region (tuple of float) – (west, east, south, north) region boundaries in [deg].

Returns

Mask, returning True for each point within the region.

Return type

numpy.ndarray of bool, shape (N)

point_in_region(p, region)[source]¶

Check if a point is contained in a rectangular geographical region.

Parameters

p (tuple of float) – (lat, lon) in [deg].
region (tuple of float) – (west, east, south, north) region boundaries in [deg].

Returns

True, if point is in region, else False.

Return type

bool

radius_to_region(lat, lon, radius)[source]¶

Get a rectangular region which fully contains a given circular region.

Parameters

lat (float) – Latitude of the center point of circular region in [deg].
lon (float) – Longitude of the center point of circular region in [deg].
radius (float) – Radius of circular region in [m].

Returns

Rectangular region as (east, west, south, north) in [deg] or None.

Return type

tuple[float, float, float, float]

geographic_midpoint(lats, lons, weights=None, depths=None, earthradius=6371000.0)[source]¶

Calculate geographic midpoints by finding the center of gravity.

This method suffers from instabilities if points are centered around the poles.

Parameters

lats (numpy.ndarray, (N)) – Latitudes in [deg].
lons (numpy.ndarray, (N)) – Longitudes in [deg].
weights (optional, numpy.ndarray, (N)) – Weighting factors.
depths (optional, numpy.ndarray, (N)) – Depths in [m].

Returns

Latitudes and longitudes of the midpoint in [deg] (and depth [m] if depths are given).

Return type

(lat, lon) or (lat, lon, depth)

geodetic_to_ecef(lat, lon, alt)[source]¶

Convert geodetic coordinates to Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates. 3 4

Parameters

lat (float) – Geodetic latitude in [deg].
lon (float) – Geodetic longitude in [deg].
alt (float) – Geodetic altitude (height) in [m] (positive for points outside the geoid).

Returns

ECEF Cartesian coordinates (X, Y, Z) in [m].

Return type

tuple[float, float, float]

3: https://en.wikipedia.org/wiki/ECEF
4: https://en.wikipedia.org/wiki/Geographic_coordinate_conversion #From_geodetic_to_ECEF_coordinates

ecef_to_geodetic(X, Y, Z)[source]¶

Convert Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates to geodetic coordinates (Ferrari’s solution).

Parameters: Z (X, Y,) – Cartesian coordinates in ECEF system in [m].
Returns: Geodetic coordinates (lat, lon, alt). Latitude and longitude are in [deg] and altitude is in [m] (positive for points outside the geoid).
Return type: tuple, float

Navigation

orthodrome¶

Navigation

`orthodrome`¶