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
11from src.orthodrome import angle_difference
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(locations, center, azimuth, azimuth_width):
269 """Filter locations by azimuth swath from a center location.
271 :param locations:ocation]): List of Locations.
272 :param center: Relative center location.
273 :param azimuth: Azimuth in [deg]. -180 to 180 or 0 to 360.
274 :param azimuth_width: Width of the swath.
275 :type locations: list
276 :type center: Location
277 :type azimuth: float
278 :type azimuth_width: float
280 :return: Filtered locations.
281 :rtype: list[Location]
282 """
283 filt_locations = []
284 for loc in locations:
285 loc_azi, _ = center.azibazi_to(loc)
286 angle = angle_difference(loc_azi, azimuth)
287 if abs(angle) <= azimuth_width / 2.:
288 filt_locations.append(loc)
289 return filt_locations
292def filter_distance(locations, reference, distance_min, distance_max):
293 """Filter location by distance to a reference point.
295 :param locations: Locations to filter.
296 :param reference: Reference location.
297 :param distance_min: Minimum distance in [m].
298 :param distance_max: Maximum distance in [m].
299 :type locations: list
300 :type reference: Location
301 :type distance_min: float
302 :type distance_max: float
304 :return: Filtered locations.
305 :rtype: list[Location]
306 """
307 return [
308 loc
309 for loc in locations
310 if distance_min <= loc.distance_to(reference) <= distance_max
311 ]
314def filter_crosstrack_distance(locations, start_path, end_path, distance_max):
315 """Filter location by distance to a great-circle path.
317 :param locations: Locations to filter.
318 :param start_path: Start of the great circle path.
319 :param end_path: End of the great circle path.
320 :param distance_max: Distance to the great-circle in [deg].
321 :type locations: list
322 :type start_path: Location
323 :type end_path: Location
324 :type distance_max: float
326 :return: Filtered locations.
327 :rtype: list[Location]
328 """
329 return [
330 loc
331 for loc in locations
332 if abs(loc.crosstrack_distance_to(
333 start_path, end_path)) <= distance_max
334 ]
337def get_offset(obj):
338 try:
339 return obj.north_shift, obj.east_shift
340 except AttributeError:
341 return 0.0, 0.0
344def get_effective_latlon(obj):
345 try:
346 return obj.effective_latlon
347 except AttributeError:
348 return obj.lat, obj.lon