Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/model/location.py: 100%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6''' 

7Representation of a geographical location, base class for stations, events, 

8etc. 

9''' 

10 

11import numpy as num 

12import math 

13 

14from pyrocko import orthodrome 

15from pyrocko.guts import Object, Float 

16 

17 

18guts_prefix = 'pf' 

19 

20d2r = math.pi / 180. 

21r2d = 1.0 / d2r 

22km = 1000. 

23 

24 

25def latlondepth_to_cartesian(lat, lon, depth): 

26 return orthodrome.geodetic_to_ecef(lat, lon, -depth) 

27 

28 

29class Location(Object): 

30 ''' 

31 Geographical location. 

32 

33 The location is given by a reference point at the earth's surface 

34 (:py:attr:`lat`, :py:attr:`lon`, :py:attr:`elevation`) and a cartesian 

35 offset from this point (:py:attr:`north_shift`, :py:attr:`east_shift`, 

36 :py:attr:`depth`). The offset corrected lat/lon coordinates of the location 

37 can be accessed though the :py:attr:`effective_latlon`, 

38 :py:attr:`effective_lat`, and :py:attr:`effective_lon` properties. 

39 ''' 

40 

41 lat = Float.T( 

42 default=0.0, 

43 optional=True, 

44 help='Latitude of reference point [deg].') 

45 

46 lon = Float.T( 

47 default=0.0, 

48 optional=True, 

49 help='Longitude of reference point [deg].') 

50 

51 north_shift = Float.T( 

52 default=0., 

53 optional=True, 

54 help='Northward cartesian offset from reference point [m].') 

55 

56 east_shift = Float.T( 

57 default=0., 

58 optional=True, 

59 help='Eastward cartesian offset from reference point [m].') 

60 

61 elevation = Float.T( 

62 default=0.0, 

63 optional=True, 

64 help='Surface elevation, above sea level [m].') 

65 

66 depth = Float.T( 

67 default=0.0, 

68 help='Depth, below surface [m].') 

69 

70 def __init__(self, **kwargs): 

71 Object.__init__(self, **kwargs) 

72 self._latlon = None 

73 

74 def __setattr__(self, name, value): 

75 if name in ('lat', 'lon', 'north_shift', 'east_shift'): 

76 self.__dict__['_latlon'] = None 

77 

78 Object.__setattr__(self, name, value) 

79 

80 @property 

81 def effective_latlon(self): 

82 ''' 

83 Property holding the offset-corrected lat/lon pair of the location. 

84 ''' 

85 

86 if self._latlon is None: 

87 if self.north_shift == 0.0 and self.east_shift == 0.0: 

88 self._latlon = self.lat, self.lon 

89 else: 

90 self._latlon = tuple(float(x) for x in orthodrome.ne_to_latlon( 

91 self.lat, self.lon, self.north_shift, self.east_shift)) 

92 

93 return self._latlon 

94 

95 @property 

96 def effective_lat(self): 

97 ''' 

98 Property holding the offset-corrected latitude of the location. 

99 ''' 

100 

101 return self.effective_latlon[0] 

102 

103 @property 

104 def effective_lon(self): 

105 ''' 

106 Property holding the offset-corrected longitude of the location. 

107 ''' 

108 

109 return self.effective_latlon[1] 

110 

111 def same_origin(self, other): 

112 ''' 

113 Check whether other location object has the same reference location. 

114 ''' 

115 

116 return self.lat == other.lat and self.lon == other.lon 

117 

118 def distance_to(self, other): 

119 ''' 

120 Compute surface distance [m] to other location object. 

121 

122 

123 ''' 

124 

125 if self.same_origin(other): 

126 other_north_shift, other_east_shift = get_offset(other) 

127 

128 return math.sqrt((self.north_shift - other_north_shift)**2 + 

129 (self.east_shift - other_east_shift)**2) 

130 

131 else: 

132 slat, slon = self.effective_latlon 

133 rlat, rlon = get_effective_latlon(other) 

134 

135 return float(orthodrome.distance_accurate50m_numpy( 

136 slat, slon, rlat, rlon)[0]) 

137 

138 def distance_3d_to(self, other): 

139 ''' 

140 Compute 3D distance [m] to other location object. 

141 

142 All coordinates are transformed to cartesian coordinates if necessary 

143 then distance is: 

144 

145 .. math:: 

146 

147 \\Delta = \\sqrt{\\Delta {\\bf x}^2 + \\Delta {\\bf y}^2 + \ 

148 \\Delta {\\bf z}^2} 

149 

150 ''' 

151 

152 if self.same_origin(other): 

153 other_north_shift, other_east_shift = get_offset(other) 

154 return math.sqrt((self.north_shift - other_north_shift)**2 + 

155 (self.east_shift - other_east_shift)**2 + 

156 (self.depth - other.depth)**2) 

157 

158 else: 

159 slat, slon = self.effective_latlon 

160 rlat, rlon = get_effective_latlon(other) 

161 

162 sx, sy, sz = latlondepth_to_cartesian(slat, slon, self.depth) 

163 rx, ry, rz = latlondepth_to_cartesian(rlat, rlon, other.depth) 

164 

165 return math.sqrt((sx-rx)**2 + (sy-ry)**2 + (sz-rz)**2) 

166 

167 def offset_to(self, other): 

168 if self.same_origin(other): 

169 other_north_shift, other_east_shift = get_offset(other) 

170 return ( 

171 other_north_shift - self.north_shift, 

172 other_east_shift - self.east_shift) 

173 

174 else: 

175 azi, bazi = self.azibazi_to(other) 

176 dist = self.distance_to(other) 

177 return dist*math.cos(azi*d2r), dist*math.sin(azi*d2r) 

178 

179 def azibazi_to(self, other): 

180 ''' 

181 Compute azimuth and backazimuth to and from other location object. 

182 ''' 

183 

184 if self.same_origin(other): 

185 other_north_shift, other_east_shift = get_offset(other) 

186 azi = r2d * math.atan2(other_east_shift - self.east_shift, 

187 other_north_shift - self.north_shift) 

188 

189 bazi = azi + 180. 

190 else: 

191 slat, slon = self.effective_latlon 

192 rlat, rlon = get_effective_latlon(other) 

193 azi, bazi = orthodrome.azibazi(slat, slon, rlat, rlon) 

194 

195 return float(azi), float(bazi) 

196 

197 def set_origin(self, lat, lon): 

198 lat = float(lat) 

199 lon = float(lon) 

200 elat, elon = self.effective_latlon 

201 n, e = orthodrome.latlon_to_ne(lat, lon, elat, elon) 

202 self.lat = lat 

203 self.lon = lon 

204 self.north_shift = float(n) 

205 self.east_shift = float(e) 

206 self._latlon = elat, elon # unchanged 

207 

208 @property 

209 def coords5(self): 

210 return num.array([ 

211 self.lat, self.lon, self.north_shift, self.east_shift, self.depth]) 

212 

213 

214def get_offset(obj): 

215 try: 

216 return obj.north_shift, obj.east_shift 

217 except AttributeError: 

218 return 0.0, 0.0 

219 

220 

221def get_effective_latlon(obj): 

222 try: 

223 return obj.effective_latlon 

224 except AttributeError: 

225 return obj.lat, obj.lon