import numpy as num import math
from pyrocko import orthodrome from pyrocko.guts import Object, Float
guts_prefix = 'pf'
d2r = math.pi / 180. r2d = 1.0 / d2r km = 1000.
def latlondepth_to_cartesian(lat, lon, depth): return orthodrome.geodetic_to_ecef(lat, lon, -depth)
class Location(Object): ''' Geographical location.
The location is given by a reference point at the earth's surface (:py:attr:`lat`, :py:attr:`lon`, :py:attr:`elevation`) and a cartesian offset from this point (:py:attr:`north_shift`, :py:attr:`east_shift`, :py:attr:`depth`). The offset corrected lat/lon coordinates of the location can be accessed though the :py:attr:`effective_latlon`, :py:attr:`effective_lat`, and :py:attr:`effective_lon` properties. '''
lat = Float.T( default=0.0, optional=True, help='latitude of reference point [deg]')
lon = Float.T( default=0.0, optional=True, help='longitude of reference point [deg]')
north_shift = Float.T( default=0., optional=True, help='northward cartesian offset from reference point [m]')
east_shift = Float.T( default=0., optional=True, help='eastward cartesian offset from reference point [m]')
elevation = Float.T( default=0.0, optional=True, help='surface elevation, above sea level [m]')
depth = Float.T( default=0.0, help='depth, below surface [m]')
def __init__(self, **kwargs): Object.__init__(self, **kwargs) self._latlon = None
def __setattr__(self, name, value): if name in ('lat', 'lon', 'north_shift', 'east_shift'): self.__dict__['_latlon'] = None
Object.__setattr__(self, name, value)
@property def effective_latlon(self): ''' Property holding the offset-corrected lat/lon pair of the location. '''
if self._latlon is None: if self.north_shift == 0.0 and self.east_shift == 0.0: self._latlon = self.lat, self.lon else: self._latlon = tuple(float(x) for x in orthodrome.ne_to_latlon( self.lat, self.lon, self.north_shift, self.east_shift))
return self._latlon
@property def effective_lat(self): ''' Property holding the offset-corrected latitude of the location. '''
return self.effective_latlon[0]
@property def effective_lon(self): ''' Property holding the offset-corrected longitude of the location. '''
return self.effective_latlon[1]
def same_origin(self, other): ''' Check whether other location object has the same reference location. '''
return self.lat == other.lat and self.lon == other.lon
def distance_to(self, other): ''' Compute surface distance [m] to other location object.
'''
if self.same_origin(other): other_north_shift, other_east_shift = get_offset(other)
return math.sqrt((self.north_shift - other_north_shift)**2 + (self.east_shift - other_east_shift)**2)
else: slat, slon = self.effective_latlon rlat, rlon = get_effective_latlon(other)
return float(orthodrome.distance_accurate50m_numpy( slat, slon, rlat, rlon)[0])
def distance_3d_to(self, other): ''' Compute 3D distance [m] to other location object.
All coordinates are transformed to cartesian coordinates if necessary then distance is:
.. math::
\\Delta = \\sqrt{\\Delta {\\bf x}^2 + \\Delta {\\bf y}^2 + \ \\Delta {\\bf z}^2}
'''
if self.same_origin(other): other_north_shift, other_east_shift = get_offset(other) return math.sqrt((self.north_shift - other_north_shift)**2 + (self.east_shift - other_east_shift)**2 + (self.depth - other.depth)**2)
else: slat, slon = self.effective_latlon rlat, rlon = get_effective_latlon(other)
sx, sy, sz = latlondepth_to_cartesian(slat, slon, self.depth) rx, ry, rz = latlondepth_to_cartesian(rlat, rlon, other.depth)
return math.sqrt((sx-rx)**2 + (sy-ry)**2 + (sz-rz)**2)
def offset_to(self, other): if self.same_origin(other): other_north_shift, other_east_shift = get_offset(other) return ( other_north_shift - self.north_shift, other_east_shift - self.east_shift)
else: azi, bazi = self.azibazi_to(other) dist = self.distance_to(other) return dist*math.cos(azi*d2r), dist*math.sin(azi*d2r)
def azibazi_to(self, other): ''' Compute azimuth and backazimuth to and from other location object. '''
if self.same_origin(other): other_north_shift, other_east_shift = get_offset(other) azi = r2d * math.atan2(other_east_shift - self.east_shift, other_north_shift - self.north_shift)
bazi = azi + 180. else: slat, slon = self.effective_latlon rlat, rlon = get_effective_latlon(other) azi, bazi = orthodrome.azibazi_numpy(slat, slon, rlat, rlon)
return float(azi), float(bazi)
def set_origin(self, lat, lon): lat = float(lat) lon = float(lon) elat, elon = self.effective_latlon n, e = orthodrome.latlon_to_ne_numpy(lat, lon, elat, elon) self.lat = lat self.lon = lon self.north_shift = float(n) self.east_shift = float(e) self._latlon = elat, elon # unchanged
@property def coords5(self): return num.array([ self.lat, self.lon, self.north_shift, self.east_shift, self.depth])
def get_offset(obj): try: return obj.north_shift, obj.east_shift except AttributeError: return 0.0, 0.0
def get_effective_latlon(obj): try: return obj.effective_latlon except AttributeError: return obj.lat, obj.lon |