1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import annotations
7import math
9import numpy as num
10from pyrocko import orthodrome
11from pyrocko.guts import Float, Object
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.
117 :param other: Other location.
118 :type other: :py:class:`~pyrocko.model.location.Location`
120 :return: Distance to another location in [m].
121 :rtype: float
122 '''
124 if self.same_origin(other):
125 other_north_shift, other_east_shift = get_offset(other)
127 return math.sqrt((self.north_shift - other_north_shift)**2 +
128 (self.east_shift - other_east_shift)**2)
130 else:
131 slat, slon = self.effective_latlon
132 rlat, rlon = get_effective_latlon(other)
134 return float(orthodrome.distance_accurate50m_numpy(
135 slat, slon, rlat, rlon)[0])
137 def distance_3d_to(self, other):
138 '''
139 Compute 3D distance [m] to other location object.
141 All coordinates are transformed to cartesian coordinates if necessary
142 then distance is:
144 .. math::
146 \\Delta = \\sqrt{\\Delta {\\bf x}^2 + \\Delta {\\bf y}^2 + \
147 \\Delta {\\bf z}^2}
149 :param other: Other location.
150 :type other: :py:class:`~pyrocko.model.location.Location`
152 :return: 3D distance to another location in [m].
153 :rtype: float
154 '''
156 if self.same_origin(other):
157 other_north_shift, other_east_shift = get_offset(other)
158 return math.sqrt((self.north_shift - other_north_shift)**2 +
159 (self.east_shift - other_east_shift)**2 +
160 (self.depth - other.depth)**2)
162 else:
163 slat, slon = self.effective_latlon
164 rlat, rlon = get_effective_latlon(other)
166 sx, sy, sz = latlondepth_to_cartesian(slat, slon, self.depth)
167 rx, ry, rz = latlondepth_to_cartesian(rlat, rlon, other.depth)
169 return math.sqrt((sx-rx)**2 + (sy-ry)**2 + (sz-rz)**2)
171 def crosstrack_distance_to(self, path_begin, path_end):
172 '''
173 Compute distance to a great-circle arc.
175 :param path_begin: Location of the start of the arc.
176 :param path_end: Location of the end of the arc.
177 :type path_begin: :py:class:`~pyrocko.model.location.Location`
178 :type path_end: :py:class:`~pyrocko.model.location.Location`
180 :return: Distance to a great circle arc in [deg].
181 :rtype: float
182 '''
183 return orthodrome.crosstrack_distance(
184 *path_begin.effective_latlon,
185 *path_end.effective_latlon,
186 *self.effective_latlon
187 )
189 def alongtrack_distance_to(self, path_begin, path_end, meter=False):
190 '''
191 Compute distance along a great-circle arc.
193 :param path_begin: Location of the start of the arc.
194 :param path_end: Location of the end of the arc.
195 :param meter: Return [m] instead of [deg].
196 :type path_begin: :py:class:`~pyrocko.model.location.Location`
197 :type path_end: :py:class:`~pyrocko.model.location.Location`
199 :return: Distance from the start of the great circle arc
200 in [deg] or [m].
201 :rtype: float
202 '''
203 if meter:
204 return orthodrome.alongtrack_distance_m(
205 *path_begin.effective_latlon,
206 *path_end.effective_latlon,
207 *self.effective_latlon
208 )
209 return orthodrome.alongtrack_distance(
210 *path_begin.effective_latlon,
211 *path_end.effective_latlon,
212 *self.effective_latlon
213 )
215 def offset_to(self, other):
216 if self.same_origin(other):
217 other_north_shift, other_east_shift = get_offset(other)
218 return (
219 other_north_shift - self.north_shift,
220 other_east_shift - self.east_shift)
222 else:
223 azi, bazi = self.azibazi_to(other)
224 dist = self.distance_to(other)
225 return dist*math.cos(azi*d2r), dist*math.sin(azi*d2r)
227 def azibazi_to(self, other):
228 '''
229 Compute azimuth and backazimuth to and from other location object.
231 :param other: Other location.
232 :type other: :py:class:`~pyrocko.model.location.Location`
234 :return: Azimuth and back-azimuth to the location in [deg].
235 :rtype: tuple[float, float]
236 '''
238 if self.same_origin(other):
239 other_north_shift, other_east_shift = get_offset(other)
240 azi = r2d * math.atan2(other_east_shift - self.east_shift,
241 other_north_shift - self.north_shift)
243 bazi = azi + 180.
244 else:
245 slat, slon = self.effective_latlon
246 rlat, rlon = get_effective_latlon(other)
247 azi, bazi = orthodrome.azibazi_numpy(slat, slon, rlat, rlon)
249 return float(azi), float(bazi)
251 def set_origin(self, lat, lon):
252 lat = float(lat)
253 lon = float(lon)
254 elat, elon = self.effective_latlon
255 n, e = orthodrome.latlon_to_ne_numpy(lat, lon, elat, elon)
256 self.lat = lat
257 self.lon = lon
258 self.north_shift = float(n)
259 self.east_shift = float(e)
260 self._latlon = elat, elon # unchanged
262 @property
263 def coords5(self):
264 return num.array([
265 self.lat, self.lon, self.north_shift, self.east_shift, self.depth])
268def filter_azimuths(
269 locations: list[Location],
270 center: Location,
271 azimuth: float,
272 azimuth_width: float,
273) -> list[Location]:
274 """Filter locations by azimuth swath.
276 Args:
277 locations (list[Location]): List of Locations.
278 center (Location): Relative center location.
279 azimuth (float): Azimuth in [deg]. -180 to 180 or 0 to 360.
280 azimuth_width (float): Width of the swath.
282 Returns:
283 list[Location]: Filtered locations.
284 """
285 OFFSET = 720.0
286 filt_locations = []
288 azimuth_min = azimuth + OFFSET - azimuth_width
289 azimuth_max = azimuth + OFFSET + azimuth_width
290 for loc in locations:
291 azimuth, _ = center.azibazi_to(loc)
292 if (azimuth_min <= azimuth + OFFSET <= azimuth_max) or (
293 azimuth_min <= azimuth % 360.0 + OFFSET <= azimuth_max
294 ):
295 filt_locations.append(loc)
296 return filt_locations
299def filter_distance(
300 locations: list[Location],
301 reference: Location,
302 distance_min: float,
303 distance_max: float,
304) -> list[Location]:
305 """Filter location by distance to a reference point.
307 Args:
308 locations (list[Location]): Locations to filter.
309 reference (Location): Reference location.
310 distance_min (float): Minimum distance in [m].
311 distance_max (float): Maximum distance in [m].
313 Returns:
314 list[Location]: Filtered locations.
315 """
316 return [
317 loc
318 for loc in locations
319 if distance_min <= loc.distance_to(reference) <= distance_max
320 ]
323def filter_crosstrack_distance(
324 locations: list[Location],
325 start_path: Location,
326 end_path: Location,
327 distance_max: float
328) -> list[Location]:
329 """Filter location by distance to a great-circle path.
331 Args:
332 locations (list[Location]): Locations to filter.
333 start_path (Location): Start of the great circle path.
334 end_path (Location): End of the great circle path.
335 distance_max (float): Distance to the great-circle in [deg].
337 Returns:
338 list[Location]: Filtered locations.
339 """
340 return [
341 loc
342 for loc in locations
343 if abs(loc.crosstrack_distance_to(
344 start_path, end_path)) <= distance_max
345 ]
348def get_offset(obj):
349 try:
350 return obj.north_shift, obj.east_shift
351 except AttributeError:
352 return 0.0, 0.0
355def get_effective_latlon(obj):
356 try:
357 return obj.effective_latlon
358 except AttributeError:
359 return obj.lat, obj.lon