Source code for obspy.geodetics.base

# -*- coding: utf-8 -*-
"""
Various geodetic utilities for ObsPy.

:copyright:
    The ObsPy Development Team (devs@obspy.org)
:license:
    GNU Lesser General Public License, Version 3
    (https://www.gnu.org/copyleft/lesser.html)
"""
import math
import warnings

import numpy as np

from obspy.core.util.misc import to_int_or_zero


# checking for geographiclib
try:
    import geographiclib  # @UnusedImport # NOQA
    from geographiclib.geodesic import Geodesic
    HAS_GEOGRAPHICLIB = True
    try:
        GEOGRAPHICLIB_VERSION_AT_LEAST_1_34 = [1, 34] <= list(map(
            to_int_or_zero, geographiclib.__version__.split(".")))
    except AttributeError:
        GEOGRAPHICLIB_VERSION_AT_LEAST_1_34 = False
except ImportError:
    HAS_GEOGRAPHICLIB = False
    GEOGRAPHICLIB_VERSION_AT_LEAST_1_34 = False


WGS84_A = 6378137.0
WGS84_F = 1 / 298.257223563


[docs] def _check_latitude(latitude, variable_name='latitude'): """ Check whether latitude is in the -90 to +90 range. """ if latitude is None: return if latitude > 90 or latitude < -90: msg = '{} out of bounds! (-90 <= {} <=90)'.format( variable_name, variable_name) raise ValueError(msg)
[docs] def _normalize_longitude(longitude): """ Normalize longitude in the -180 to +180 range. """ if longitude is None: return while longitude > 180: longitude -= 360 while longitude < -180: longitude += 360 return longitude
[docs] def calc_vincenty_inverse(lat1, lon1, lat2, lon2, a=WGS84_A, f=WGS84_F): """ Vincenty Inverse Solution of Geodesics on the Ellipsoid. Computes the distance between two geographic points on the WGS84 ellipsoid and the forward and backward azimuths between these points. :param lat1: Latitude of point A in degrees (positive for northern, negative for southern hemisphere) :param lon1: Longitude of point A in degrees (positive for eastern, negative for western hemisphere) :param lat2: Latitude of point B in degrees (positive for northern, negative for southern hemisphere) :param lon2: Longitude of point B in degrees (positive for eastern, negative for western hemisphere) :param a: Radius of Earth in m. Uses the value for WGS84 by default. :param f: Flattening of Earth. Uses the value for WGS84 by default. :return: (Great circle distance in m, azimuth A->B in degrees, azimuth B->A in degrees) :raises: This method may have no solution between two nearly antipodal points; an iteration limit traps this case and a ``StopIteration`` exception will be raised. .. note:: This code is based on an implementation incorporated in Matplotlib Basemap Toolkit 0.9.5 http://sourceforge.net/projects/\ matplotlib/files/matplotlib-toolkits/basemap-0.9.5/ (matplotlib/toolkits/basemap/greatcircle.py) Algorithm from Geocentric Datum of Australia Technical Manual. * https://www.icsm.gov.au/publications * https://www.icsm.gov.au/publications/gda2020-technical-manual-v16 It states:: There are a number of formulae that are available to calculate accurate geodetic positions, azimuths and distances on the ellipsoid. Vincenty's formulae (Vincenty, 1975) may be used for lines ranging from a few cm to nearly 20,000 km, with millimetre accuracy. The formulae have been extensively tested for the Australian region, by comparison with results from other formulae (Rainsford, 1955 & Sodano, 1965). * Inverse problem: azimuth and distance from known latitudes and longitudes * Direct problem: Latitude and longitude from known position, azimuth and distance. """ # Check inputs _check_latitude(lat1, 'lat1') lon1 = _normalize_longitude(lon1) _check_latitude(lat2, 'lat2') lon2 = _normalize_longitude(lon2) b = a * (1 - f) # semiminor axis if math.isclose(lat1, lat2) and math.isclose(lon1, lon2): return 0.0, 0.0, 0.0 # convert latitudes and longitudes to radians: lat1 = math.radians(lat1) lon1 = math.radians(lon1) lat2 = math.radians(lat2) lon2 = math.radians(lon2) tan_u1 = (1 - f) * math.tan(lat1) tan_u2 = (1 - f) * math.tan(lat2) u_1 = math.atan(tan_u1) u_2 = math.atan(tan_u2) dlon = lon2 - lon1 last_dlon = -4000000.0 # an impossible value omega = dlon # Iterate until no significant change in dlon or iterlimit has been # reached (http://www.movable-type.co.uk/scripts/latlong-vincenty.html) iterlimit = 100 try: while (last_dlon < -3000000.0 or dlon != 0 and abs((last_dlon - dlon) / dlon) > 1.0e-9): sqr_sin_sigma = pow(math.cos(u_2) * math.sin(dlon), 2) + \ pow((math.cos(u_1) * math.sin(u_2) - math.sin(u_1) * math.cos(u_2) * math.cos(dlon)), 2) sin_sigma = math.sqrt(sqr_sin_sigma) cos_sigma = math.sin(u_1) * math.sin(u_2) + math.cos(u_1) * \ math.cos(u_2) * math.cos(dlon) sigma = math.atan2(sin_sigma, cos_sigma) sin_alpha = math.cos(u_1) * math.cos(u_2) * math.sin(dlon) / \ sin_sigma sqr_cos_alpha = 1 - sin_alpha * sin_alpha if math.isclose(sqr_cos_alpha, 0): # Equatorial line cos2sigma_m = 0 else: cos2sigma_m = cos_sigma - \ (2 * math.sin(u_1) * math.sin(u_2) / sqr_cos_alpha) c = (f / 16) * sqr_cos_alpha * (4 + f * (4 - 3 * sqr_cos_alpha)) last_dlon = dlon dlon = omega + (1 - c) * f * sin_alpha * \ (sigma + c * sin_sigma * (cos2sigma_m + c * cos_sigma * (-1 + 2 * pow(cos2sigma_m, 2)))) iterlimit -= 1 if iterlimit < 0: # iteration limit reached raise StopIteration except ValueError: # usually "math domain error" raise StopIteration u2 = sqr_cos_alpha * (a * a - b * b) / (b * b) _a = 1 + (u2 / 16384) * (4096 + u2 * (-768 + u2 * (320 - 175 * u2))) _b = (u2 / 1024) * (256 + u2 * (-128 + u2 * (74 - 47 * u2))) delta_sigma = _b * sin_sigma * \ (cos2sigma_m + (_b / 4) * (cos_sigma * (-1 + 2 * pow(cos2sigma_m, 2)) - (_b / 6) * cos2sigma_m * (-3 + 4 * sqr_sin_sigma) * (-3 + 4 * pow(cos2sigma_m, 2)))) dist = b * _a * (sigma - delta_sigma) alpha12 = math.atan2( (math.cos(u_2) * math.sin(dlon)), (math.cos(u_1) * math.sin(u_2) - math.sin(u_1) * math.cos(u_2) * math.cos(dlon))) alpha21 = math.atan2( (math.cos(u_1) * math.sin(dlon)), (-math.sin(u_1) * math.cos(u_2) + math.cos(u_1) * math.sin(u_2) * math.cos(dlon))) if alpha12 < 0.0: alpha12 = alpha12 + (2.0 * math.pi) if alpha12 > (2.0 * math.pi): alpha12 = alpha12 - (2.0 * math.pi) alpha21 = alpha21 + math.pi if alpha21 < 0.0: alpha21 = alpha21 + (2.0 * math.pi) if alpha21 > (2.0 * math.pi): alpha21 = alpha21 - (2.0 * math.pi) # convert to degrees: alpha12 = alpha12 * 360 / (2.0 * math.pi) alpha21 = alpha21 * 360 / (2.0 * math.pi) return dist, alpha12, alpha21
[docs] def gps2dist_azimuth(lat1, lon1, lat2, lon2, a=WGS84_A, f=WGS84_F): """ Computes the distance between two geographic points on the WGS84 ellipsoid and the forward and backward azimuths between these points. :param lat1: Latitude of point A in degrees (positive for northern, negative for southern hemisphere) :param lon1: Longitude of point A in degrees (positive for eastern, negative for western hemisphere) :param lat2: Latitude of point B in degrees (positive for northern, negative for southern hemisphere) :param lon2: Longitude of point B in degrees (positive for eastern, negative for western hemisphere) :param a: Radius of Earth in m. Uses the value for WGS84 by default. :param f: Flattening of Earth. Uses the value for WGS84 by default. :return: (Great circle distance in m, azimuth A->B in degrees, azimuth B->A in degrees) .. note:: This function will check if you have installed the Python module `geographiclib <http://geographiclib.sf.net>`_ - a very fast module for converting between geographic, UTM, UPS, MGRS, and geocentric coordinates, for geoid calculations, and for solving geodesic problems. Otherwise the locally implemented Vincenty's Inverse formulae (:func:`obspy.geodetics.base.calc_vincenty_inverse`) is used which has known limitations for two nearly antipodal points and is ca. 4x slower. """ if HAS_GEOGRAPHICLIB: _check_latitude(lat1, 'lat1') _check_latitude(lat2, 'lat2') result = Geodesic(a=a, f=f).Inverse(lat1, lon1, lat2, lon2) azim = result['azi1'] if azim < 0: azim += 360 bazim = result['azi2'] + 180 return (result['s12'], azim, bazim) else: try: values = calc_vincenty_inverse(lat1, lon1, lat2, lon2, a, f) if np.all(np.isnan(values)): raise StopIteration return values except StopIteration: msg = ("Catching unstable calculation on antipodes. " "The currently used Vincenty's Inverse formulae " "has known limitations for two nearly antipodal points. " "Install the Python module 'geographiclib' to solve this " "issue.") warnings.warn(msg) return (20004314.5, 0.0, 0.0) except ValueError as e: raise e
[docs] def kilometers2degrees(kilometer, radius=6371.0): """ Convenience function to convert kilometers to degrees assuming a perfectly spherical Earth. :type kilometer: float :param kilometer: Distance in kilometers :type radius: float, optional :param radius: Radius of the Earth used for the calculation. :rtype: float :return: Distance in degrees as a floating point number. .. rubric:: Example >>> from obspy.geodetics import kilometers2degrees >>> kilometers2degrees(300) 2.6979648177561915 """ return kilometer / (2.0 * radius * math.pi / 360.0)
kilometer2degrees = kilometers2degrees
[docs] def degrees2kilometers(degrees, radius=6371.0): """ Convenience function to convert (great circle) degrees to kilometers assuming a perfectly spherical Earth. :type degrees: float :param degrees: Distance in (great circle) degrees :type radius: float, optional :param radius: Radius of the Earth used for the calculation. :rtype: float :return: Distance in kilometers as a floating point number. .. rubric:: Example >>> from obspy.geodetics import degrees2kilometers >>> degrees2kilometers(1) 111.19492664455873 """ return degrees * (2.0 * radius * math.pi / 360.0)
[docs] def locations2degrees(lat1, long1, lat2, long2): """ Convenience function to calculate the great circle distance between two points on a spherical Earth. This method uses the Vincenty formula in the special case of a spherical Earth. For more accurate values use the geodesic distance calculations of geopy (https://github.com/geopy/geopy). :type lat1: float or :class:`numpy.ndarray` :param lat1: Latitude(s) of point 1 in degrees :type long1: float or :class:`numpy.ndarray` :param long1: Longitude(s) of point 1 in degrees :type lat2: float or :class:`numpy.ndarray` :param lat2: Latitude(s) of point 2 in degrees :type long2: float or :class:`numpy.ndarray` :param long2: Longitude(s) of point 2 in degrees :rtype: float or :class:`numpy.ndarray` :return: Distance in degrees as a floating point number, or numpy array of element-wise distances in degrees .. rubric:: Example >>> from obspy.geodetics import locations2degrees >>> locations2degrees(5, 5, 10, 10) # doctest: +ELLIPSIS 7.03970141917538... """ # broadcast explicitly here so it raises once instead of somewhere in the # middle if things can't be broadcast lat1, lat2, long1, long2 = np.broadcast_arrays(lat1, lat2, long1, long2) # Convert to radians. lat1 = np.radians(np.asarray(lat1)) lat2 = np.radians(np.asarray(lat2)) long1 = np.radians(np.asarray(long1)) long2 = np.radians(np.asarray(long2)) long_diff = long2 - long1 gd = np.degrees( np.arctan2( np.sqrt(( np.cos(lat2) * np.sin(long_diff)) ** 2 + (np.cos(lat1) * np.sin(lat2) - np.sin(lat1) * np.cos(lat2) * np.cos(long_diff)) ** 2), np.sin(lat1) * np.sin(lat2) + np.cos(lat1) * np.cos(lat2) * np.cos(long_diff))) return gd
[docs] def mean_longitude(longitudes): """ Compute sample mean longitude, assuming longitude in degrees from -180 to 180. >>> lons = (-170.5, -178.3, 166) >>> np.mean(lons) # doctest: +SKIP -60.933 >>> mean_longitude(lons) # doctest: +ELLIPSIS 179.08509... :type longitudes: :class:`~numpy.ndarray` (or list, ..) :param longitudes: Geographical longitude values ranging from -180 to 180 in degrees. """ from scipy.stats import circmean mean_longitude = circmean(np.array(longitudes), low=-180, high=180) mean_longitude = _normalize_longitude(mean_longitude) return mean_longitude
[docs] def inside_geobounds(obj, minlatitude=None, maxlatitude=None, minlongitude=None, maxlongitude=None, latitude=None, longitude=None, minradius=None, maxradius=None): """ Check whether an object is within a given latitude and/or longitude range, or within a given distance range from a reference geographic point. The object must have ``latitude`` and ``longitude`` attributes, expressed in degrees. :type obj: object :param obj: An object with `latitude` and `longitude` attributes. :type minlatitude: float :param minlatitude: Minimum latitude in degrees. :type maxlatitude: float :param maxlatitude: Maximum latitude in degrees. If this value is smaller than ``minlatitude``, then 360 degrees are added to this value (i.e., wrapping around latitude of +/- 180 degrees) :type minlongitude: float :param minlongitude: Minimum longitude in degrees. :type maxlongitude: float :param maxlongitude: Minimum longitude in degrees. :type latitude: float :param latitude: Latitude of the reference point, in degrees, for distance range selection. :type longitude: float :param longitude: Longitude of the reference point, in degrees, for distance range selection. :type minradius: float :param minradius: Minimum distance, in degrees, from the reference geographic point defined by the latitude and longitude parameters. :type maxradius: float :param maxradius: Maximum distance, in degrees, from the reference geographic point defined by the latitude and longitude parameters. :return: ``True`` if the object is within the given range, ``False`` otherwise. .. rubric:: Example >>> from obspy.geodetics import inside_geobounds >>> from obspy import read_events >>> ev = read_events()[0] >>> orig = ev.origins[0] >>> inside_geobounds(orig, minlatitude=40, maxlatitude=42) True >>> inside_geobounds(orig, minlatitude=40, maxlatitude=42, ... minlongitude=78, maxlongitude=79) False >>> inside_geobounds(orig, latitude=40, longitude=80, ... minradius=1, maxradius=10) True """ if not hasattr(obj, 'latitude') or not hasattr(obj, 'longitude'): raise AttributeError( 'Object must have "latitude" and "longitude" attributes.') olatitude = obj.latitude _check_latitude(olatitude, 'obj.latitude') _check_latitude(minlatitude, 'minlatitude') _check_latitude(maxlatitude, 'maxlatitude') _check_latitude(latitude, 'latitude') # Make sure longitudes are between -180 to 180 degrees olongitude = _normalize_longitude(obj.longitude) minlongitude = _normalize_longitude(minlongitude) maxlongitude = _normalize_longitude(maxlongitude) longitude = _normalize_longitude(longitude) if minlatitude is not None: if olatitude is None or olatitude < minlatitude: return False if maxlatitude is not None: if olatitude is None or olatitude > maxlatitude: return False # Wrap longitude around +/- 180°, if necessary if None not in [minlongitude, maxlongitude] \ and maxlongitude < minlongitude: maxlongitude += 360 if olongitude is not None and olongitude < minlongitude: olongitude += 360 if minlongitude is not None: if olongitude is None or olongitude < minlongitude: return False if maxlongitude is not None: if olongitude is None or olongitude > maxlongitude: return False if all([coord is not None for coord in (latitude, longitude, olatitude, olongitude)]): distance = locations2degrees(latitude, longitude, olatitude, olongitude) if minradius is not None and distance < minradius: return False if maxradius is not None and distance > maxradius: return False return True
if __name__ == '__main__': import doctest doctest.testmod(exclude_empty=True)