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 depth = Float.T(
57 default=0.0,
58 help='depth, below surface [m]')
60 elevation = Float.T(
61 default=0.0,
62 optional=True,
63 help='surface elevation, above sea level [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 def copy_location(self):
76 return Location(
77 lat=self.lat,
78 lon=self.lon,
79 north_shift=self.north_shift,
80 east_shift=self.east_shift,
81 depth=self.depth,
82 elevation=self.elevation)
84 @property
85 def location_attribute_tuple(self):
86 return (
87 self.lat,
88 self.lon,
89 self.north_shift,
90 self.east_shift,
91 self.depth,
92 self.elevation)
94 def same_location(self, other):
95 '''
96 Check if locations are equal.
98 Locations are treated as equal when all attributes are equal.
99 '''
101 return self.location_attribute_tuple == other.location_attribute_tuple
103 @property
104 def effective_latlon(self):
105 '''
106 Property holding the offset-corrected lat/lon pair of the location.
107 '''
109 if self._latlon is None:
110 if self.north_shift == 0.0 and self.east_shift == 0.0:
111 self._latlon = self.lat, self.lon
112 else:
113 self._latlon = tuple(float(x) for x in orthodrome.ne_to_latlon(
114 self.lat, self.lon, self.north_shift, self.east_shift))
116 return self._latlon
118 @property
119 def effective_lat(self):
120 '''
121 Property holding the offset-corrected latitude of the location.
122 '''
124 return self.effective_latlon[0]
126 @property
127 def effective_lon(self):
128 '''
129 Property holding the offset-corrected longitude of the location.
130 '''
132 return self.effective_latlon[1]
134 def same_origin(self, other):
135 '''
136 Check whether other location object has the same reference location.
137 '''
139 return self.lat == other.lat and self.lon == other.lon
141 def distance_to(self, other):
142 '''
143 Compute surface distance [m] to other location object.
146 '''
148 if self.same_origin(other):
149 other_north_shift, other_east_shift = get_offset(other)
151 return math.sqrt((self.north_shift - other_north_shift)**2 +
152 (self.east_shift - other_east_shift)**2)
154 else:
155 slat, slon = self.effective_latlon
156 rlat, rlon = get_effective_latlon(other)
158 return float(orthodrome.distance_accurate50m_numpy(
159 slat, slon, rlat, rlon)[0])
161 def distance_3d_to(self, other):
162 '''
163 Compute 3D distance [m] to other location object.
165 All coordinates are transformed to cartesian coordinates if necessary
166 then distance is:
168 .. math::
170 \\Delta = \\sqrt{\\Delta {\\bf x}^2 + \\Delta {\\bf y}^2 + \
171 \\Delta {\\bf z}^2}
173 '''
175 if self.same_origin(other):
176 other_north_shift, other_east_shift = get_offset(other)
177 return math.sqrt((self.north_shift - other_north_shift)**2 +
178 (self.east_shift - other_east_shift)**2 +
179 (self.depth - other.depth)**2)
181 else:
182 slat, slon = self.effective_latlon
183 rlat, rlon = get_effective_latlon(other)
185 sx, sy, sz = latlondepth_to_cartesian(slat, slon, self.depth)
186 rx, ry, rz = latlondepth_to_cartesian(rlat, rlon, other.depth)
188 return math.sqrt((sx-rx)**2 + (sy-ry)**2 + (sz-rz)**2)
190 def offset_to(self, other):
191 if self.same_origin(other):
192 other_north_shift, other_east_shift = get_offset(other)
193 return (
194 other_north_shift - self.north_shift,
195 other_east_shift - self.east_shift)
197 else:
198 azi, bazi = self.azibazi_to(other)
199 dist = self.distance_to(other)
200 return dist*math.cos(azi*d2r), dist*math.sin(azi*d2r)
202 def azibazi_to(self, other):
203 '''
204 Compute azimuth and backazimuth to and from other location object.
205 '''
207 if self.same_origin(other):
208 other_north_shift, other_east_shift = get_offset(other)
209 azi = r2d * math.atan2(other_east_shift - self.east_shift,
210 other_north_shift - self.north_shift)
212 bazi = azi + 180.
213 else:
214 slat, slon = self.effective_latlon
215 rlat, rlon = get_effective_latlon(other)
216 azi, bazi = orthodrome.azibazi_numpy(slat, slon, rlat, rlon)
218 return float(azi), float(bazi)
220 def set_origin(self, lat, lon):
221 lat = float(lat)
222 lon = float(lon)
223 elat, elon = self.effective_latlon
224 n, e = orthodrome.latlon_to_ne_numpy(lat, lon, elat, elon)
225 self.lat = lat
226 self.lon = lon
227 self.north_shift = float(n)
228 self.east_shift = float(e)
229 self._latlon = elat, elon # unchanged
231 @property
232 def coords5(self):
233 return num.array([
234 self.lat, self.lon, self.north_shift, self.east_shift, self.depth])
237def get_offset(obj):
238 try:
239 return obj.north_shift, obj.east_shift
240 except AttributeError:
241 return 0.0, 0.0
244def get_effective_latlon(obj):
245 try:
246 return obj.effective_latlon
247 except AttributeError:
248 return obj.lat, obj.lon