1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6import numpy as num 

7import math 

8 

9from pyrocko import orthodrome 

10from pyrocko.guts import Object, Float 

11 

12 

13guts_prefix = 'pf' 

14 

15d2r = math.pi / 180. 

16r2d = 1.0 / d2r 

17km = 1000. 

18 

19 

20def latlondepth_to_cartesian(lat, lon, depth): 

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

22 

23 

24class Location(Object): 

25 ''' 

26 Geographical location. 

27 

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 ''' 

35 

36 lat = Float.T( 

37 default=0.0, 

38 optional=True, 

39 help='latitude of reference point [deg]') 

40 

41 lon = Float.T( 

42 default=0.0, 

43 optional=True, 

44 help='longitude of reference point [deg]') 

45 

46 north_shift = Float.T( 

47 default=0., 

48 optional=True, 

49 help='northward cartesian offset from reference point [m]') 

50 

51 east_shift = Float.T( 

52 default=0., 

53 optional=True, 

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

55 

56 elevation = Float.T( 

57 default=0.0, 

58 optional=True, 

59 help='surface elevation, above sea level [m]') 

60 

61 depth = Float.T( 

62 default=0.0, 

63 help='depth, below surface [m]') 

64 

65 def __init__(self, **kwargs): 

66 Object.__init__(self, **kwargs) 

67 self._latlon = None 

68 

69 def __setattr__(self, name, value): 

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

71 self.__dict__['_latlon'] = None 

72 

73 Object.__setattr__(self, name, value) 

74 

75 @property 

76 def effective_latlon(self): 

77 ''' 

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

79 ''' 

80 

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)) 

87 

88 return self._latlon 

89 

90 @property 

91 def effective_lat(self): 

92 ''' 

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

94 ''' 

95 

96 return self.effective_latlon[0] 

97 

98 @property 

99 def effective_lon(self): 

100 ''' 

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

102 ''' 

103 

104 return self.effective_latlon[1] 

105 

106 def same_origin(self, other): 

107 ''' 

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

109 ''' 

110 

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

112 

113 def distance_to(self, other): 

114 ''' 

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

116 

117 

118 ''' 

119 

120 if self.same_origin(other): 

121 other_north_shift, other_east_shift = get_offset(other) 

122 

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

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

125 

126 else: 

127 slat, slon = self.effective_latlon 

128 rlat, rlon = get_effective_latlon(other) 

129 

130 return float(orthodrome.distance_accurate50m_numpy( 

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

132 

133 def distance_3d_to(self, other): 

134 ''' 

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

136 

137 All coordinates are transformed to cartesian coordinates if necessary 

138 then distance is: 

139 

140 .. math:: 

141 

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

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

144 

145 ''' 

146 

147 if self.same_origin(other): 

148 other_north_shift, other_east_shift = get_offset(other) 

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

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

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

152 

153 else: 

154 slat, slon = self.effective_latlon 

155 rlat, rlon = get_effective_latlon(other) 

156 

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

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

159 

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

161 

162 def offset_to(self, other): 

163 if self.same_origin(other): 

164 other_north_shift, other_east_shift = get_offset(other) 

165 return ( 

166 other_north_shift - self.north_shift, 

167 other_east_shift - self.east_shift) 

168 

169 else: 

170 azi, bazi = self.azibazi_to(other) 

171 dist = self.distance_to(other) 

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

173 

174 def azibazi_to(self, other): 

175 ''' 

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

177 ''' 

178 

179 if self.same_origin(other): 

180 other_north_shift, other_east_shift = get_offset(other) 

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

182 other_north_shift - self.north_shift) 

183 

184 bazi = azi + 180. 

185 else: 

186 slat, slon = self.effective_latlon 

187 rlat, rlon = get_effective_latlon(other) 

188 azi, bazi = orthodrome.azibazi_numpy(slat, slon, rlat, rlon) 

189 

190 return float(azi), float(bazi) 

191 

192 def set_origin(self, lat, lon): 

193 lat = float(lat) 

194 lon = float(lon) 

195 elat, elon = self.effective_latlon 

196 n, e = orthodrome.latlon_to_ne_numpy(lat, lon, elat, elon) 

197 self.lat = lat 

198 self.lon = lon 

199 self.north_shift = float(n) 

200 self.east_shift = float(e) 

201 self._latlon = elat, elon # unchanged 

202 

203 @property 

204 def coords5(self): 

205 return num.array([ 

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

207 

208 

209def get_offset(obj): 

210 try: 

211 return obj.north_shift, obj.east_shift 

212 except AttributeError: 

213 return 0.0, 0.0 

214 

215 

216def get_effective_latlon(obj): 

217 try: 

218 return obj.effective_latlon 

219 except AttributeError: 

220 return obj.lat, obj.lon