Coverage for /usr/local/lib/python3.13/dist-packages/pyrocko/model/location.py: 86%
116 statements
« prev ^ index » next coverage.py v7.6.0, created at 2025-12-04 10:41 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7Representation of a geographical location, base class for stations, events,
8etc.
9'''
10import math
12import numpy as num
13from pyrocko import orthodrome
14from pyrocko.guts import Float, Object
16guts_prefix = 'pf'
18d2r = math.pi / 180.
19r2d = 1.0 / d2r
20km = 1000.
23def latlondepth_to_cartesian(lat, lon, depth):
24 return orthodrome.geodetic_to_ecef(lat, lon, -depth)
27class Location(Object):
28 '''
29 Geographical location.
31 The location is given by a reference point at the earth's surface
32 (:py:attr:`lat`, :py:attr:`lon`, :py:attr:`elevation`) and a cartesian
33 offset from this point (:py:attr:`north_shift`, :py:attr:`east_shift`,
34 :py:attr:`depth`). The offset corrected lat/lon coordinates of the location
35 can be accessed though the :py:attr:`effective_latlon`,
36 :py:attr:`effective_lat`, and :py:attr:`effective_lon` properties.
37 '''
39 lat = Float.T(
40 default=0.0,
41 optional=True,
42 help='Latitude of reference point [deg].')
44 lon = Float.T(
45 default=0.0,
46 optional=True,
47 help='Longitude of reference point [deg].')
49 north_shift = Float.T(
50 default=0.,
51 optional=True,
52 help='Northward cartesian offset from reference point [m].')
54 east_shift = Float.T(
55 default=0.,
56 optional=True,
57 help='Eastward cartesian offset from reference point [m].')
59 elevation = Float.T(
60 default=0.0,
61 optional=True,
62 help='Surface elevation, above sea level [m].')
64 depth = Float.T(
65 default=0.0,
66 help='Depth, below surface [m].')
68 def __init__(self, **kwargs):
69 Object.__init__(self, **kwargs)
70 self._latlon = None
72 def __setattr__(self, name, value):
73 if name in ('lat', 'lon', 'north_shift', 'east_shift'):
74 self.__dict__['_latlon'] = None
76 Object.__setattr__(self, name, value)
78 @property
79 def effective_latlon(self):
80 '''
81 Property holding the offset-corrected lat/lon pair of the location.
82 '''
84 if self._latlon is None:
85 if self.north_shift == 0.0 and self.east_shift == 0.0:
86 self._latlon = self.lat, self.lon
87 else:
88 self._latlon = tuple(float(x) for x in orthodrome.ne_to_latlon(
89 self.lat, self.lon, self.north_shift, self.east_shift))
91 return self._latlon
93 @property
94 def effective_lat(self):
95 '''
96 Property holding the offset-corrected latitude of the location.
97 '''
99 return self.effective_latlon[0]
101 @property
102 def effective_lon(self):
103 '''
104 Property holding the offset-corrected longitude of the location.
105 '''
107 return self.effective_latlon[1]
109 def same_origin(self, other):
110 '''
111 Check whether other location object has the same reference location.
112 '''
114 return self.lat == other.lat and self.lon == other.lon
116 def distance_to(self, other):
117 '''
118 Compute surface distance [m] to other location object.
120 :param other: Other location.
121 :type other: :py:class:`~pyrocko.model.location.Location`
123 :return: Distance to another location in [m].
124 :rtype: float
125 '''
127 if self.same_origin(other):
128 other_north_shift, other_east_shift = get_offset(other)
130 return math.sqrt((self.north_shift - other_north_shift)**2 +
131 (self.east_shift - other_east_shift)**2)
133 else:
134 slat, slon = self.effective_latlon
135 rlat, rlon = get_effective_latlon(other)
137 return float(orthodrome.distance_accurate50m_numpy(
138 slat, slon, rlat, rlon)[0])
140 def distance_3d_to(self, other):
141 '''
142 Compute 3D distance [m] to other location object.
144 All coordinates are transformed to cartesian coordinates if necessary
145 then distance is:
147 .. math::
149 \\Delta = \\sqrt{\\Delta {\\bf x}^2 + \\Delta {\\bf y}^2 + \
150 \\Delta {\\bf z}^2}
152 :param other: Other location.
153 :type other: :py:class:`~pyrocko.model.location.Location`
155 :return: 3D distance to another location in [m].
156 :rtype: float
157 '''
159 if self.same_origin(other):
160 other_north_shift, other_east_shift = get_offset(other)
161 return math.sqrt((self.north_shift - other_north_shift)**2 +
162 (self.east_shift - other_east_shift)**2 +
163 (self.depth - other.depth)**2)
165 else:
166 slat, slon = self.effective_latlon
167 rlat, rlon = get_effective_latlon(other)
169 sx, sy, sz = latlondepth_to_cartesian(slat, slon, self.depth)
170 rx, ry, rz = latlondepth_to_cartesian(rlat, rlon, other.depth)
172 return math.sqrt((sx-rx)**2 + (sy-ry)**2 + (sz-rz)**2)
174 def ecef(self):
175 lat, lon = self.effective_latlon
176 return num.array(
177 latlondepth_to_cartesian(lat, lon, -self.elevation + self.depth),
178 dtype=float)
180 def crosstrack_distance_to(self, path_begin, path_end):
181 '''
182 Compute distance to a great-circle arc.
184 :param path_begin: Location of the start of the arc.
185 :param path_end: Location of the end of the arc.
186 :type path_begin: :py:class:`~pyrocko.model.location.Location`
187 :type path_end: :py:class:`~pyrocko.model.location.Location`
189 :return: Distance to a great circle arc in [deg].
190 :rtype: float
191 '''
192 return orthodrome.crosstrack_distance(
193 path_begin.effective_latlon[0], path_begin.effective_latlon[1],
194 path_end.effective_latlon[0], path_end.effective_latlon[1],
195 self.effective_latlon[0], self.effective_latlon[1]
196 )
198 def alongtrack_distance_to(self, path_begin, path_end, meter=False):
199 '''
200 Compute distance along a great-circle arc.
202 :param path_begin: Location of the start of the arc.
203 :param path_end: Location of the end of the arc.
204 :param meter: Return [m] instead of [deg].
205 :type path_begin: :py:class:`~pyrocko.model.location.Location`
206 :type path_end: :py:class:`~pyrocko.model.location.Location`
208 :return: Distance from the start of the great circle arc
209 in [deg] or [m].
210 :rtype: float
211 '''
212 func = orthodrome.alongtrack_distance
213 if meter:
214 func = orthodrome.alongtrack_distance_m
216 return func(
217 path_begin.effective_latlon[0], path_begin.effective_latlon[1],
218 path_end.effective_latlon[0], path_end.effective_latlon[1],
219 self.effective_latlon[0], self.effective_latlon[1]
220 )
222 def offset_to(self, other):
223 if self.same_origin(other):
224 other_north_shift, other_east_shift = get_offset(other)
225 return (
226 other_north_shift - self.north_shift,
227 other_east_shift - self.east_shift)
229 else:
230 azi, bazi = self.azibazi_to(other)
231 dist = self.distance_to(other)
232 return dist*math.cos(azi*d2r), dist*math.sin(azi*d2r)
234 def azibazi_to(self, other):
235 '''
236 Compute azimuth and backazimuth to and from other location object.
238 :param other: Other location.
239 :type other: :py:class:`~pyrocko.model.location.Location`
241 :return: Azimuth and back-azimuth to the location in [deg].
242 :rtype: tuple[float, float]
243 '''
245 if self.same_origin(other):
246 other_north_shift, other_east_shift = get_offset(other)
247 azi = r2d * math.atan2(other_east_shift - self.east_shift,
248 other_north_shift - self.north_shift)
250 bazi = azi + 180.
251 else:
252 slat, slon = self.effective_latlon
253 rlat, rlon = get_effective_latlon(other)
254 azi, bazi = orthodrome.azibazi(slat, slon, rlat, rlon)
256 return float(azi), float(bazi)
258 def set_origin(self, lat, lon):
259 lat = float(lat)
260 lon = float(lon)
261 elat, elon = self.effective_latlon
262 n, e = orthodrome.latlon_to_ne(lat, lon, elat, elon)
263 self.lat = lat
264 self.lon = lon
265 self.north_shift = float(n)
266 self.east_shift = float(e)
267 self._latlon = elat, elon # unchanged
269 @property
270 def coords5(self):
271 return num.array([
272 self.lat, self.lon, self.north_shift, self.east_shift, self.depth])
275def filter_azimuths(locations, center, azimuth, azimuth_width):
276 '''
277 Filter locations by azimuth swath from a center location.
279 :param locations: list[Location]): List of Locations.
280 :param center: Relative center location.
281 :param azimuth: Azimuth in [deg]. -180 to 180 or 0 to 360.
282 :param azimuth_width: Width of the swath.
283 :type locations: list
284 :type center: Location
285 :type azimuth: float
286 :type azimuth_width: float
288 :return: Filtered locations.
289 :rtype: list[Location]
290 '''
291 filt_locations = []
292 for loc in locations:
293 loc_azi, _ = center.azibazi_to(loc)
294 angle = orthodrome.angle_difference(loc_azi, azimuth)
295 if abs(angle) <= azimuth_width / 2.:
296 filt_locations.append(loc)
297 return filt_locations
300def filter_distance(locations, reference, distance_min, distance_max):
301 '''Filter location by distance to a reference point.
303 :param locations: Locations to filter.
304 :param reference: Reference location.
305 :param distance_min: Minimum distance in [m].
306 :param distance_max: Maximum distance in [m].
307 :type locations: list
308 :type reference: Location
309 :type distance_min: float
310 :type distance_max: float
312 :return: Filtered locations.
313 :rtype: list[Location]
314 '''
315 return [
316 loc
317 for loc in locations
318 if distance_min <= loc.distance_to(reference) <= distance_max
319 ]
322def filter_crosstrack_distance(locations, start_path, end_path, distance_max):
323 '''Filter location by distance to a great-circle path.
325 :param locations: Locations to filter.
326 :param start_path: Start of the great circle path.
327 :param end_path: End of the great circle path.
328 :param distance_max: Distance to the great-circle in [deg].
329 :type locations: list
330 :type start_path: Location
331 :type end_path: Location
332 :type distance_max: float
334 :return: Filtered locations.
335 :rtype: list[Location]
336 '''
337 return [
338 loc
339 for loc in locations
340 if abs(loc.crosstrack_distance_to(
341 start_path, end_path)) <= distance_max
342 ]
345def get_offset(obj):
346 try:
347 return obj.north_shift, obj.east_shift
348 except AttributeError:
349 return 0.0, 0.0
352def get_effective_latlon(obj):
353 try:
354 return obj.effective_latlon
355 except AttributeError:
356 return obj.lat, obj.lon