1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5import math
7import numpy as num
8from pyrocko import orthodrome
9from pyrocko.guts import Float, Object
11guts_prefix = 'pf'
13d2r = math.pi / 180.
14r2d = 1.0 / d2r
15km = 1000.
18def latlondepth_to_cartesian(lat, lon, depth):
19 return orthodrome.geodetic_to_ecef(lat, lon, -depth)
22class Location(Object):
23 '''
24 Geographical location.
26 The location is given by a reference point at the earth's surface
27 (:py:attr:`lat`, :py:attr:`lon`, :py:attr:`elevation`) and a cartesian
28 offset from this point (:py:attr:`north_shift`, :py:attr:`east_shift`,
29 :py:attr:`depth`). The offset corrected lat/lon coordinates of the location
30 can be accessed though the :py:attr:`effective_latlon`,
31 :py:attr:`effective_lat`, and :py:attr:`effective_lon` properties.
32 '''
34 lat = Float.T(
35 default=0.0,
36 optional=True,
37 help='latitude of reference point [deg]')
39 lon = Float.T(
40 default=0.0,
41 optional=True,
42 help='longitude of reference point [deg]')
44 north_shift = Float.T(
45 default=0.,
46 optional=True,
47 help='northward cartesian offset from reference point [m]')
49 east_shift = Float.T(
50 default=0.,
51 optional=True,
52 help='eastward cartesian offset from reference point [m]')
54 elevation = Float.T(
55 default=0.0,
56 optional=True,
57 help='surface elevation, above sea level [m]')
59 depth = Float.T(
60 default=0.0,
61 help='depth, below surface [m]')
63 def __init__(self, **kwargs):
64 Object.__init__(self, **kwargs)
65 self._latlon = None
67 def __setattr__(self, name, value):
68 if name in ('lat', 'lon', 'north_shift', 'east_shift'):
69 self.__dict__['_latlon'] = None
71 Object.__setattr__(self, name, value)
73 @property
74 def effective_latlon(self):
75 '''
76 Property holding the offset-corrected lat/lon pair of the location.
77 '''
79 if self._latlon is None:
80 if self.north_shift == 0.0 and self.east_shift == 0.0:
81 self._latlon = self.lat, self.lon
82 else:
83 self._latlon = tuple(float(x) for x in orthodrome.ne_to_latlon(
84 self.lat, self.lon, self.north_shift, self.east_shift))
86 return self._latlon
88 @property
89 def effective_lat(self):
90 '''
91 Property holding the offset-corrected latitude of the location.
92 '''
94 return self.effective_latlon[0]
96 @property
97 def effective_lon(self):
98 '''
99 Property holding the offset-corrected longitude of the location.
100 '''
102 return self.effective_latlon[1]
104 def same_origin(self, other):
105 '''
106 Check whether other location object has the same reference location.
107 '''
109 return self.lat == other.lat and self.lon == other.lon
111 def distance_to(self, other):
112 '''
113 Compute surface distance [m] to other location object.
115 :param other: Other location.
116 :type other: :py:class:`~pyrocko.model.location.Location`
118 :return: Distance to another location in [m].
119 :rtype: float
120 '''
122 if self.same_origin(other):
123 other_north_shift, other_east_shift = get_offset(other)
125 return math.sqrt((self.north_shift - other_north_shift)**2 +
126 (self.east_shift - other_east_shift)**2)
128 else:
129 slat, slon = self.effective_latlon
130 rlat, rlon = get_effective_latlon(other)
132 return float(orthodrome.distance_accurate50m_numpy(
133 slat, slon, rlat, rlon)[0])
135 def distance_3d_to(self, other):
136 '''
137 Compute 3D distance [m] to other location object.
139 All coordinates are transformed to cartesian coordinates if necessary
140 then distance is:
142 .. math::
144 \\Delta = \\sqrt{\\Delta {\\bf x}^2 + \\Delta {\\bf y}^2 + \
145 \\Delta {\\bf z}^2}
147 :param other: Other location.
148 :type other: :py:class:`~pyrocko.model.location.Location`
150 :return: 3D distance to another location in [m].
151 :rtype: float
152 '''
154 if self.same_origin(other):
155 other_north_shift, other_east_shift = get_offset(other)
156 return math.sqrt((self.north_shift - other_north_shift)**2 +
157 (self.east_shift - other_east_shift)**2 +
158 (self.depth - other.depth)**2)
160 else:
161 slat, slon = self.effective_latlon
162 rlat, rlon = get_effective_latlon(other)
164 sx, sy, sz = latlondepth_to_cartesian(slat, slon, self.depth)
165 rx, ry, rz = latlondepth_to_cartesian(rlat, rlon, other.depth)
167 return math.sqrt((sx-rx)**2 + (sy-ry)**2 + (sz-rz)**2)
169 def crosstrack_distance_to(self, path_begin, path_end):
170 '''
171 Compute distance to a great-circle arc.
173 :param path_begin: Location of the start of the arc.
174 :param path_end: Location of the end of the arc.
175 :type path_begin: :py:class:`~pyrocko.model.location.Location`
176 :type path_end: :py:class:`~pyrocko.model.location.Location`
178 :return: Distance to a great circle arc in [deg].
179 :rtype: float
180 '''
181 return orthodrome.crosstrack_distance(
182 path_begin.effective_latlon[0], path_begin.effective_latlon[1],
183 path_end.effective_latlon[0], path_end.effective_latlon[1],
184 self.effective_latlon[0], self.effective_latlon[1]
185 )
187 def alongtrack_distance_to(self, path_begin, path_end, meter=False):
188 '''
189 Compute distance along a great-circle arc.
191 :param path_begin: Location of the start of the arc.
192 :param path_end: Location of the end of the arc.
193 :param meter: Return [m] instead of [deg].
194 :type path_begin: :py:class:`~pyrocko.model.location.Location`
195 :type path_end: :py:class:`~pyrocko.model.location.Location`
197 :return: Distance from the start of the great circle arc
198 in [deg] or [m].
199 :rtype: float
200 '''
201 func = orthodrome.alongtrack_distance
202 if meter:
203 func = orthodrome.alongtrack_distance_m
205 return func(
206 path_begin.effective_latlon[0], path_begin.effective_latlon[1],
207 path_end.effective_latlon[0], path_end.effective_latlon[1],
208 self.effective_latlon[0], self.effective_latlon[1]
209 )
211 def offset_to(self, other):
212 if self.same_origin(other):
213 other_north_shift, other_east_shift = get_offset(other)
214 return (
215 other_north_shift - self.north_shift,
216 other_east_shift - self.east_shift)
218 else:
219 azi, bazi = self.azibazi_to(other)
220 dist = self.distance_to(other)
221 return dist*math.cos(azi*d2r), dist*math.sin(azi*d2r)
223 def azibazi_to(self, other):
224 '''
225 Compute azimuth and backazimuth to and from other location object.
227 :param other: Other location.
228 :type other: :py:class:`~pyrocko.model.location.Location`
230 :return: Azimuth and back-azimuth to the location in [deg].
231 :rtype: tuple[float, float]
232 '''
234 if self.same_origin(other):
235 other_north_shift, other_east_shift = get_offset(other)
236 azi = r2d * math.atan2(other_east_shift - self.east_shift,
237 other_north_shift - self.north_shift)
239 bazi = azi + 180.
240 else:
241 slat, slon = self.effective_latlon
242 rlat, rlon = get_effective_latlon(other)
243 azi, bazi = orthodrome.azibazi_numpy(slat, slon, rlat, rlon)
245 return float(azi), float(bazi)
247 def set_origin(self, lat, lon):
248 lat = float(lat)
249 lon = float(lon)
250 elat, elon = self.effective_latlon
251 n, e = orthodrome.latlon_to_ne_numpy(lat, lon, elat, elon)
252 self.lat = lat
253 self.lon = lon
254 self.north_shift = float(n)
255 self.east_shift = float(e)
256 self._latlon = elat, elon # unchanged
258 @property
259 def coords5(self):
260 return num.array([
261 self.lat, self.lon, self.north_shift, self.east_shift, self.depth])
264def filter_azimuths(locations, center, azimuth, azimuth_width):
265 """Filter locations by azimuth swath from a center location.
267 :param locations:ocation]): List of Locations.
268 :param center: Relative center location.
269 :param azimuth: Azimuth in [deg]. -180 to 180 or 0 to 360.
270 :param azimuth_width: Width of the swath.
271 :type locations: list
272 :type center: Location
273 :type azimuth: float
274 :type azimuth_width: float
276 :return: Filtered locations.
277 :rtype: list[Location]
278 """
279 filt_locations = []
280 for loc in locations:
281 loc_azi, _ = center.azibazi_to(loc)
282 angle = orthodrome.angle_difference(loc_azi, azimuth)
283 if abs(angle) <= azimuth_width / 2.:
284 filt_locations.append(loc)
285 return filt_locations
288def filter_distance(locations, reference, distance_min, distance_max):
289 """Filter location by distance to a reference point.
291 :param locations: Locations to filter.
292 :param reference: Reference location.
293 :param distance_min: Minimum distance in [m].
294 :param distance_max: Maximum distance in [m].
295 :type locations: list
296 :type reference: Location
297 :type distance_min: float
298 :type distance_max: float
300 :return: Filtered locations.
301 :rtype: list[Location]
302 """
303 return [
304 loc
305 for loc in locations
306 if distance_min <= loc.distance_to(reference) <= distance_max
307 ]
310def filter_crosstrack_distance(locations, start_path, end_path, distance_max):
311 """Filter location by distance to a great-circle path.
313 :param locations: Locations to filter.
314 :param start_path: Start of the great circle path.
315 :param end_path: End of the great circle path.
316 :param distance_max: Distance to the great-circle in [deg].
317 :type locations: list
318 :type start_path: Location
319 :type end_path: Location
320 :type distance_max: float
322 :return: Filtered locations.
323 :rtype: list[Location]
324 """
325 return [
326 loc
327 for loc in locations
328 if abs(loc.crosstrack_distance_to(
329 start_path, end_path)) <= distance_max
330 ]
333def get_offset(obj):
334 try:
335 return obj.north_shift, obj.east_shift
336 except AttributeError:
337 return 0.0, 0.0
340def get_effective_latlon(obj):
341 try:
342 return obj.effective_latlon
343 except AttributeError:
344 return obj.lat, obj.lon