1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import numpy as num
7import math
9from pyrocko import orthodrome
10from pyrocko.guts import Object, Float
13guts_prefix = 'pf'
15d2r = math.pi / 180.
16r2d = 1.0 / d2r
17km = 1000.
20def latlondepth_to_cartesian(lat, lon, depth):
21 return orthodrome.geodetic_to_ecef(lat, lon, -depth)
24class Location(Object):
25 '''
26 Geographical location.
28 The location is given by a reference point at the earth's surface
29 (:py:attr:`lat`, :py:attr:`lon`, :py:attr:`elevation`) and a cartesian
30 offset from this point (:py:attr:`north_shift`, :py:attr:`east_shift`,
31 :py:attr:`depth`). The offset corrected lat/lon coordinates of the location
32 can be accessed though the :py:attr:`effective_latlon`,
33 :py:attr:`effective_lat`, and :py:attr:`effective_lon` properties.
34 '''
36 lat = Float.T(
37 default=0.0,
38 optional=True,
39 help='latitude of reference point [deg]')
41 lon = Float.T(
42 default=0.0,
43 optional=True,
44 help='longitude of reference point [deg]')
46 north_shift = Float.T(
47 default=0.,
48 optional=True,
49 help='northward cartesian offset from reference point [m]')
51 east_shift = Float.T(
52 default=0.,
53 optional=True,
54 help='eastward cartesian offset from reference point [m]')
56 elevation = Float.T(
57 default=0.0,
58 optional=True,
59 help='surface elevation, above sea level [m]')
61 depth = Float.T(
62 default=0.0,
63 help='depth, below surface [m]')
65 def __init__(self, **kwargs):
66 Object.__init__(self, **kwargs)
67 self._latlon = None
69 def __setattr__(self, name, value):
70 if name in ('lat', 'lon', 'north_shift', 'east_shift'):
71 self.__dict__['_latlon'] = None
73 Object.__setattr__(self, name, value)
75 @property
76 def effective_latlon(self):
77 '''
78 Property holding the offset-corrected lat/lon pair of the location.
79 '''
81 if self._latlon is None:
82 if self.north_shift == 0.0 and self.east_shift == 0.0:
83 self._latlon = self.lat, self.lon
84 else:
85 self._latlon = tuple(float(x) for x in orthodrome.ne_to_latlon(
86 self.lat, self.lon, self.north_shift, self.east_shift))
88 return self._latlon
90 @property
91 def effective_lat(self):
92 '''
93 Property holding the offset-corrected latitude of the location.
94 '''
96 return self.effective_latlon[0]
98 @property
99 def effective_lon(self):
100 '''
101 Property holding the offset-corrected longitude of the location.
102 '''
104 return self.effective_latlon[1]
106 def same_origin(self, other):
107 '''
108 Check whether other location object has the same reference location.
109 '''
111 return self.lat == other.lat and self.lon == other.lon
113 def distance_to(self, other):
114 '''
115 Compute surface distance [m] to other location object.
118 '''
120 if self.same_origin(other):
121 other_north_shift, other_east_shift = get_offset(other)
123 return math.sqrt((self.north_shift - other_north_shift)**2 +
124 (self.east_shift - other_east_shift)**2)
126 else:
127 slat, slon = self.effective_latlon
128 rlat, rlon = get_effective_latlon(other)
130 return float(orthodrome.distance_accurate50m_numpy(
131 slat, slon, rlat, rlon)[0])
133 def distance_3d_to(self, other):
134 '''
135 Compute 3D distance [m] to other location object.
137 All coordinates are transformed to cartesian coordinates if necessary
138 then distance is:
140 .. math::
142 \\Delta = \\sqrt{\\Delta {\\bf x}^2 + \\Delta {\\bf y}^2 + \
143 \\Delta {\\bf z}^2}
145 '''
147 if self.same_origin(other):
148 other_north_shift, other_east_shift = get_offset(other)
149 return math.sqrt((self.north_shift - other_north_shift)**2 +
150 (self.east_shift - other_east_shift)**2 +
151 (self.depth - other.depth)**2)
153 else:
154 slat, slon = self.effective_latlon
155 rlat, rlon = get_effective_latlon(other)
157 sx, sy, sz = latlondepth_to_cartesian(slat, slon, self.depth)
158 rx, ry, rz = latlondepth_to_cartesian(rlat, rlon, other.depth)
160 return math.sqrt((sx-rx)**2 + (sy-ry)**2 + (sz-rz)**2)
162 def offset_to(self, other):
163 if self.same_origin(other):
164 other_north_shift, other_east_shift = get_offset(other)
165 return (
166 other_north_shift - self.north_shift,
167 other_east_shift - self.east_shift)
169 else:
170 azi, bazi = self.azibazi_to(other)
171 dist = self.distance_to(other)
172 return dist*math.cos(azi*d2r), dist*math.sin(azi*d2r)
174 def azibazi_to(self, other):
175 '''
176 Compute azimuth and backazimuth to and from other location object.
177 '''
179 if self.same_origin(other):
180 other_north_shift, other_east_shift = get_offset(other)
181 azi = r2d * math.atan2(other_east_shift - self.east_shift,
182 other_north_shift - self.north_shift)
184 bazi = azi + 180.
185 else:
186 slat, slon = self.effective_latlon
187 rlat, rlon = get_effective_latlon(other)
188 azi, bazi = orthodrome.azibazi_numpy(slat, slon, rlat, rlon)
190 return float(azi), float(bazi)
192 def set_origin(self, lat, lon):
193 lat = float(lat)
194 lon = float(lon)
195 elat, elon = self.effective_latlon
196 n, e = orthodrome.latlon_to_ne_numpy(lat, lon, elat, elon)
197 self.lat = lat
198 self.lon = lon
199 self.north_shift = float(n)
200 self.east_shift = float(e)
201 self._latlon = elat, elon # unchanged
203 @property
204 def coords5(self):
205 return num.array([
206 self.lat, self.lon, self.north_shift, self.east_shift, self.depth])
209def get_offset(obj):
210 try:
211 return obj.north_shift, obj.east_shift
212 except AttributeError:
213 return 0.0, 0.0
216def get_effective_latlon(obj):
217 try:
218 return obj.effective_latlon
219 except AttributeError:
220 return obj.lat, obj.lon