« back to all changes in this revision
Viewing changes to blitzortung/calc.py
- browse files at revision 391.1.1
- compare with another revision
- download diff
- download tarball
- view history from revision 391.1.1
- Committer: Andreas Würl
- Date: 2016-09-14 20:39:51 UTC
- mto: This revision was merged to the branch mainline in revision 392.
- Revision ID: git-v1:0577158d36adaee55ba3f8a63918b601e8a4edd4
remove pandas dependency
- files removed:
- blitzortung/builder/strike_cluster.py
- blitzortung/clustering
- files modified:
- .travis.yml
added
removed
blitzortung/calc.py
1
# -*- coding: utf8 -*-
2
3
"""
4
5
Copyright 2014-2016 Andreas Würl
6
7
Licensed under the Apache License, Version 2.0 (the "License");
8
you may not use this file except in compliance with the License.
9
You may obtain a copy of the License at
10
11
http://www.apache.org/licenses/LICENSE-2.0
12
13
Unless required by applicable law or agreed to in writing, software
14
distributed under the License is distributed on an "AS IS" BASIS,
15
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16
See the License for the specific language governing permissions and
17
limitations under the License.
18
19
"""
20
21
import math
22
import datetime
23
import itertools
24
import collections
25
import numpy as np
26
import pandas as pd
27
from injector import singleton
28
29
from blitzortung import data, builder, types
30
31
32
@singleton
33
class SignalVelocity(object):
34
"""
35
class for time/distance conversion regarding the reduced speed of light
36
"""
37
38
# speed of light in m / ns
39
__c0 = 0.299792458
40
41
__c_reduction_permille = 2.5
42
43
__c = (1 - 0.001 * __c_reduction_permille) * __c0
44
45
def get_distance_time(self, distance):
46
""" return the time in nanoseconds for a given distance in meters """
47
return int(distance / self.__c)
48
49
def get_time_distance(self, time_ns):
50
""" return the distance in meters for a given time interval in nanoseconds """
51
return time_ns * self.__c
52
53
54
class SimulatedData(object):
55
def __init__(self, x_coord_or_point, y_coord=None):
56
self.strike_location = types.Point(x_coord_or_point, y_coord)
57
self.signal_velocity = SignalVelocity()
58
self.event_builder = builder.Event()
59
self.timestamp = pd.Timestamp(datetime.datetime.utcnow())
60
61
def get_timestamp(self):
62
return self.timestamp
63
64
def get_event_at(self, x_coord_or_point, y_coord=None, distance_offset=0.0):
65
event_location = types.Point(x_coord_or_point, y_coord)
66
distance = self.strike_location.distance_to(event_location)
67
68
nanosecond_offset = self.signal_velocity.get_distance_time(distance + distance_offset)
69
70
self.event_builder.set_x(x_coord_or_point)
71
self.event_builder.set_y(y_coord)
72
self.event_builder.set_timestamp(self.timestamp, nanosecond_offset)
73
74
return self.event_builder.build()
75
76
def get_source_event(self):
77
self.event_builder.set_x(self.strike_location.x)
78
self.event_builder.set_y(self.strike_location.y)
79
self.event_builder.set_timestamp(self.timestamp)
80
return self.event_builder.build()
81
82
83
class ThreePointSolution(data.Event):
84
def __init__(self, reference_event, azimuth, distance, signal_velocity):
85
location = reference_event.geodesic_shift(azimuth, distance)
86
distance = reference_event.distance_to(location)
87
88
total_nanoseconds = reference_event.timestamp.value
89
total_nanoseconds -= signal_velocity.get_distance_time(distance)
90
self.signal_velocity = signal_velocity
91
92
super(ThreePointSolution, self).__init__(pd.Timestamp(total_nanoseconds), location)
93
94
def get_residual_time_at(self, event):
95
distance = self.distance_to(event)
96
distance_runtime = self.signal_velocity.get_distance_time(distance)
97
98
measured_runtime = event.timestamp.value - self.timestamp.value
99
100
return (measured_runtime - distance_runtime) / 1000.0
101
102
def get_total_residual_time_of(self, events):
103
residual_time_sum = 0.0
104
for event in events:
105
residual_time = self.get_residual_time_at(event)
106
107
residual_time_sum += residual_time * residual_time
108
109
return math.sqrt(residual_time_sum) / len(events)
110
111
def __str__(self):
112
return super(ThreePointSolution, self).__str__()
113
114
115
class ThreePointSolver(object):
116
"""
117
calculates the exact coordinates of the intersection of two hyperbola defined by three event points/times
118
119
The solution is calculated in a polar coordinate system which has its origin in the location
120
of the first event point
121
"""
122
123
def __init__(self, events):
124
if len(events) != 3:
125
raise ValueError("ThreePointSolution requires three events")
126
127
self.events = events
128
from . import INJECTOR
129
130
self.signal_velocity = INJECTOR.get(SignalVelocity)
131
132
distance_0_1 = events[0].distance_to(events[1])
133
distance_0_2 = events[0].distance_to(events[2])
134
135
azimuth_0_1 = events[0].azimuth_to(events[1])
136
azimuth_0_2 = events[0].azimuth_to(events[2])
137
138
time_difference_0_1 = events[0].ns_difference_to(events[1])
139
time_difference_0_2 = events[0].ns_difference_to(events[2])
140
141
time_distance_0_1 = self.signal_velocity.get_time_distance(time_difference_0_1)
142
time_distance_0_2 = self.signal_velocity.get_time_distance(time_difference_0_2)
143
144
self.solutions = self.solve(time_distance_0_1, distance_0_1, azimuth_0_1, time_distance_0_2,
145
distance_0_2, azimuth_0_2)
146
147
def get_solution(self):
148
solution_count = len(self.solutions)
149
150
if solution_count > 1:
151
solutions = {}
152
for solution in self.solutions:
153
solutions[solution] = solution.get_total_residual_time_of(self.events)
154
print(" 3PointSolution:", solution, solutions[solution])
155
return min(solutions, key=solutions.get)
156
elif solution_count == 1:
157
return self.solutions[0]
158
else:
159
return None
160
161
def solve(self, d1, g1, azimuth1, d2, g2, azimuth2):
162
163
phi1 = self.azimuth_to_angle(azimuth1)
164
phi2 = self.azimuth_to_angle(azimuth2)
165
166
(p1, q1) = self.calculate_p_q(d1, g1)
167
(p2, q2) = self.calculate_p_q(d2, g2)
168
169
(cosine, sine) = self.calculate_angle_projection(phi1, phi2)
170
171
denominator = q1 * q1 + q2 * q2 - 2 * q1 * q2 * cosine
172
173
root_argument = q2 * q2 * (-self.square(p1 - p2) + denominator) * sine * sine
174
175
solutions = []
176
177
if root_argument < 0.0 or denominator == 0.0:
178
print("%.1f %.1f %.1f°, %.1f %.1f %.1f°" % (d1, g1, azimuth1, d2, g2, azimuth2))
179
return solutions
180
181
part_1 = (-p1 * q1 + p2 * q1 + (p1 - p2) * q2 * cosine) / denominator
182
part_2 = math.sqrt(root_argument) / denominator
183
184
solution_angles = []
185
if abs(part_1 + part_2) <= 1.0:
186
solution_angles.append(math.acos(part_1 + part_2))
187
solution_angles.append(-math.acos(part_1 + part_2))
188
if abs(part_1 - part_2) <= 1.0:
189
solution_angles.append(math.acos(part_1 - part_2))
190
solution_angles.append(-math.acos(part_1 - part_2))
191
192
for solution_angle in solution_angles:
193
if self.is_angle_valid_for_hyperbola(solution_angle, d1, g1, phi1) and \
194
self.is_angle_valid_for_hyperbola(solution_angle, d2, g2, phi2):
195
196
solution_distance = self.hyperbola_radius(solution_angle, d1, g1, 0)
197
solution_azimuth = self.angle_to_azimuth(solution_angle + phi1)
198
solution_a = ThreePointSolution(self.events[0], solution_azimuth, solution_distance,
199
self.signal_velocity)
200
201
solution_distance = self.hyperbola_radius(solution_angle, d2, g2, phi2 - phi1)
202
solution_azimuth = self.angle_to_azimuth(solution_angle + phi1)
203
solution_b = ThreePointSolution(self.events[0], solution_azimuth, solution_distance,
204
self.signal_velocity)
205
206
if solution_a.has_same_location(solution_b):
207
solutions.append(solution_a)
208
209
return solutions
210
211
@staticmethod
212
def calculate_angle_projection(phi1, phi2):
213
angle = phi2 - phi1
214
return math.cos(angle), math.sin(angle)
215
216
@staticmethod
217
def hyperbola_radius(theta, d, g, phi):
218
return 0.5 * (g * g - d * d) / (d + g * math.cos(theta - phi))
219
220
def is_angle_valid_for_hyperbola(self, angle, d, g, phi):
221
angle -= phi
222
223
if angle <= math.pi:
224
angle += 2 * math.pi
225
if angle > math.pi:
226
angle -= 2 * math.pi
227
228
asymptotic_angle = self.calculate_asymptotic_angle(d, g)
229
230
return -asymptotic_angle < angle < asymptotic_angle
231
232
@staticmethod
233
def calculate_asymptotic_angle(d, g):
234
return math.acos(-d / g)
235
236
@staticmethod
237
def calculate_p_q(d, g):
238
denominator = g * g - d * d
239
return d / denominator, g / denominator
240
241
@staticmethod
242
def square(x):
243
return x * x
244
245
@staticmethod
246
def azimuth_to_angle(azimuth):
247
return math.pi / 2 - azimuth
248
249
@staticmethod
250
def angle_to_azimuth(angle):
251
return math.pi / 2 - angle
252
253
254
class FitSeed(object):
255
def __init__(self, events, signal_velocity):
256
self.events = events
257
self.signal_velocity = signal_velocity
258
259
self.solutions = {}
260
261
def iterate_combinations(self):
262
263
for combination in itertools.combinations(self.events[1:5], 2):
264
self.find_three_point_solution([self.events[0]] + list(combination))
265
266
def find_three_point_solution(self, selected_events):
267
three_point_solver = ThreePointSolver(selected_events)
268
269
solution = three_point_solver.get_solution()
270
if solution:
271
self.solutions[solution] = solution.get_total_residual_time_of(self.events)
272
273
def get_seed_event(self):
274
275
self.iterate_combinations()
276
277
if self.solutions:
278
for solution, residual_time in self.solutions.iteritems():
279
print("%.1f %s" % (residual_time, str(solution)))
280
return min(self.solutions, key=self.solutions.get)
281
else:
282
return None
283
284
285
class FitParameter:
286
Time, Longitude, Latitude = range(3)
287
288
def __init__(self):
289
pass
290
291
292
class LeastSquareFit(object):
293
TIME_FACTOR = 1000.0
294
295
def __init__(self, three_point_solution, events, signal_velocity):
296
self.n_dim = 3
297
self.m_dim = len(events)
298
self.events = events
299
self.time_reference = events[0].timestamp
300
self.signal_velocity = signal_velocity
301
302
self.parameters = collections.OrderedDict()
303
self.parameters[FitParameter.Time] = self.calculate_time_value(three_point_solution.timestamp)
304
self.parameters[FitParameter.Longitude] = three_point_solution.x
305
self.parameters[FitParameter.Latitude] = three_point_solution.y
306
307
self.a_matrix = np.zeros((self.m_dim, self.n_dim))
308
self.b_vector = np.zeros(self.m_dim)
309
self.residuals = np.zeros(self.m_dim)
310
311
self.least_square_sum = None
312
self.previous_least_square_sum = None
313
314
self.successful = False
315
self.iteration_count = 0
316
317
def get_parameter(self, parameter):
318
return self.parameters[parameter]
319
320
def calculate_time_value(self, timestamp):
321
return (timestamp.value - self.time_reference.value) / self.TIME_FACTOR
322
323
def calculate_partial_derivative(self, event_index, parameter_index):
324
325
delta = 0.001
326
327
self.parameters[parameter_index] += delta
328
329
slope = (self.residuals[event_index] - self.get_residual_time_at(self.events[event_index])) / delta
330
331
self.parameters[parameter_index] -= delta
332
333
return slope
334
335
def perform_fit_step(self):
336
337
self.iteration_count += 1
338
339
self.initialize_data()
340
341
(solution, residues, rank, singular_values) = np.linalg.lstsq(self.a_matrix, self.b_vector)
342
343
self.update_parameters(solution)
344
345
self.calculate_least_square_sum()
346
347
def initialize_data(self):
348
349
for index, event in enumerate(self.events):
350
self.residuals[index] = self.get_residual_time_at(event)
351
352
for event_index in range(self.m_dim):
353
for parameter_index in self.parameters:
354
self.a_matrix[event_index][parameter_index] = self.calculate_partial_derivative(event_index,
355
parameter_index)
356
self.b_vector[event_index] = self.get_residual_time_at(self.events[event_index])
357
358
def update_parameters(self, solution):
359
for parameter_index, delta in enumerate(solution):
360
self.parameters[parameter_index] += delta
361
362
def calculate_least_square_sum(self):
363
residual_time_sum = 0.0
364
365
for event in self.events:
366
residual_time = self.get_residual_time_at(event)
367
368
residual_time_sum += residual_time * residual_time
369
370
overdetermination_factor = max(1, self.m_dim - self.n_dim)
371
372
residual_time_sum /= overdetermination_factor
373
374
self.previous_least_square_sum = self.least_square_sum
375
self.least_square_sum = residual_time_sum
376
377
return residual_time_sum
378
379
def get_residual_time_at(self, event):
380
location = self.get_location()
381
distance = location.distance_to(event)
382
distance_runtime = self.signal_velocity.get_distance_time(distance) / self.TIME_FACTOR
383
384
measured_runtime = self.calculate_time_value(event.timestamp) - self.parameters[FitParameter.Time]
385
386
return measured_runtime - distance_runtime
387
388
def get_location(self):
389
return types.Point(self.parameters[FitParameter.Longitude], self.parameters[FitParameter.Latitude])
390
391
def get_timestamp(self):
392
return pd.Timestamp(self.time_reference.value + int(round(self.parameters[FitParameter.Time] * 1000, 0)))
393
394
def get_least_square_sum(self):
395
return self.least_square_sum
396
397
def get_least_square_change(self):
398
if self.least_square_sum and self.previous_least_square_sum:
399
return self.least_square_sum / self.previous_least_square_sum - 1.0
400
return float("nan")
401
402
def requires_another_iteration(self):
403
if self.least_square_sum and self.previous_least_square_sum:
404
if self.least_square_sum / self.previous_least_square_sum > 1.1:
405
return False
406
if abs(self.get_least_square_change()) < 10e-4:
407
self.successful = True
408
return False
409
410
if self.iteration_count >= 20:
411
return False
412
413
return True
414
415
def is_successful(self):
416
return self.successful
417
418
def get_solution(self):
419
event_builder = builder.Event()
420
421
location = self.get_location()
422
event_builder.set_x(location.x())
423
event_builder.set_y(location.y())
424
event_builder.set_timestamp(self.get_timestamp())
425
426
return event_builder.build()
427
428
def get_number_of_iterations(self):
429
return self.iteration_count
Loggerhead is a web-based interface for Breezy
Version: 2.0.1