1#!/usr/bin/env python3
2
3"""Use the ORB_HISTO table to geolocate timeseries data."""
4
5import importlib
6import logging
7import collections
8from datetime import timedelta
9
10import numpy as np
11
12from chart.db.connection import db_connect
13from chart.common.util import timedelta_to_us
14from chart.project import settings
15
16def Geoloc(*args, **kwargs):
17 """Return the correct Geoloc object for this project."""
18 if settings.GEOLOC_CONSTRUCTOR is None:
19 return None
20
21 elif settings.GEOLOC_CONSTRUCTOR is True:
22 return EPSGeoloc(*args, **kwargs)
23
24 elif isinstance(settings.GEOLOC_CONSTRUCTOR, str):
25 module_name, _, class_name = settings.GEOLOC_CONSTRUCTOR.rpartition('.')
26 Geoloc_class = getattr(importlib.import_module(module_name), class_name)
27 return Geoloc_class(*args, **kwargs)
28
29 else:
30 raise ValueError('Unknown type of geoloc constructor')
31
32class CannotGeolocate(Exception):
33 """Requested sid+time combination cannot be found in ORB_HISTO table."""
34
35 def __init__(self, sid, sensing_time, message=None):
36 super(CannotGeolocate, self).__init__()
37 self.sid = sid
38 self.sensing_time = sensing_time
39 self.message = message
40
41 def __str__(self):
42 if self.message is None:
43 return 'Cannot geolocate {sid} at {time}'.format(
44 sid=self.sid, time=self.sensing_time)
45
46 else:
47 return 'Cannot geolocate {sid} at {time}: {message}'.format(
48 sid=self.sid, time=self.sensing_time, message=self.message)
49
50
51class EPSGeoloc:
52 """Create an object which can be used to geolocate one or more timestamps for a satellite.
53
54 Initialised with a spacecraft ID and general time range,
55 then the individual functions can be called (must be called in time order)
56 to retrieve geolocations for specific time values.
57 """
58
59 db_conn = None
60
61 def __init__(self, sid, sensing_start, sensing_stop, width=None, height=None, pos_vel=False):
62 """Create a geolocator object.
63
64 `sid` (str) :: Source ID
65 `sensing_start` (datetime) :: Minimum validity time for this geolocator
66 `sensing_stop` (datetime) :: Maximum validity time for this geolocator
67 `width` (int) ::
68 `height` (int) ::
69 `pos_vel` (boolean) :: Allow velocity information to be retrieved in addition to position
70 """
71
72 # logging.debug('Construct geoloc for {sid} from {start} to {stop}'.format(
73 # sid=sid, start=sensing_start, stop=sensing_stop))
74
75 self.sid = sid
76
77 if EPSGeoloc.db_conn is None:
78 EPSGeoloc.db_conn = db_connect('ORB_HISTO')
79
80 if pos_vel:
81 self.cur = EPSGeoloc.db_conn.query(
82 'SELECT utc_time, x, y, z, x_velocity, y_velocity, z_velocity '
83 'FROM orb_histo '
84 'WHERE scid=:scid AND utc_time>=:start_time AND utc_time<=:stop_time '
85 'ORDER BY utc_time',
86 scid=sid.scid,
87 start_time=sensing_start - timedelta(minutes=1),
88 stop_time=sensing_stop + timedelta(minutes=1))
89 self.this_utc = None
90 self.this_posvel = None
91 self.prev_utc = None
92 self.prev_posvel = None
93
94 else:
95 self.cur = EPSGeoloc.db_conn.query(
96 'SELECT utc_time, longitude, geod_latitude, geoc_radius '
97 'FROM orb_histo '
98 'WHERE scid=:scid AND utc_time>=:start_time AND utc_time<=:stop_time '
99 'ORDER BY utc_time',
100 scid=sid.scid,
101 start_time=sensing_start - timedelta(minutes=1),
102 stop_time=sensing_stop + timedelta(minutes=1))
103 # (invalid attribute name)
104 self.Geo = collections.namedtuple('Geo', 'utc lon lat height') # pylint: disable=C0103
105 self.delta_lon = None
106 self.prev_lon = None
107 self.delta_lat = None
108 self.delta_height = None
109 self.prev_lat = None
110 self.prev_geo = None
111 self.this_geo = None
112
113 self.time_diff = None
114
115 if width is not None:
116 self.width_mul = width / 360.0
117
118 if height is not None:
119 self.height_mul = height / 180.0
120
121 self.height = height
122 self.width = width
123
124 def lat_lon_height(self, sensing_time):
125 """Compute latitude, longitude and height for `sensing_time`."""
126
127 while self.this_geo is None or self.prev_geo is None or sensing_time >= self.this_geo.utc:
128 self.prev_geo = self.this_geo
129 row = self.cur.fetchone()
130 if row is None:
131 raise CannotGeolocate(self.sid, sensing_time)
132
133 self.this_geo = self.Geo(*row)
134 if self.prev_geo is not None:
135 self.delta_lon = ((self.this_geo.lon + 180.0) - (self.prev_geo.lon + 180.0)) % 360.0
136 if self.delta_lon > 180.0:
137 self.delta_lon -= 360.0
138
139 self.prev_lon = self.prev_geo.lon + 180.0
140
141 self.delta_lat = self.this_geo.lat - self.prev_geo.lat
142
143 self.delta_height = self.this_geo.height - self.prev_geo.height
144
145 self.time_diff = timedelta_to_us(self.this_geo.utc - self.prev_geo.utc)
146
147 fract = timedelta_to_us(sensing_time - self.prev_geo.utc) / self.time_diff
148
149 lon = ((self.prev_lon + self.delta_lon * fract) % 360.0) - 180.0
150 lat = self.prev_geo.lat + self.delta_lat * fract
151 height = self.prev_geo.height + self.delta_height * fract
152
153 return lat, lon, height
154
155 def lat_lon(self, sensing_time):
156 """Compute latitude, longitude for `sensing_time`.
157
158 Note that for any Geolocator this function must be called in chronological order.
159 """
160
161 # logging.debug('Geolocating {time}'.format(time=sensing_time))
162
163 while self.this_geo is None or self.prev_geo is None or sensing_time >= self.this_geo.utc:
164 # logging.debug('this geo {current} prev_geo {prev}'.format(
165 # current=self.this_geo, prev=self.prev_geo))
166
167 self.prev_geo = self.this_geo
168 row = self.cur.fetchone()
169 if row is None:
170 raise CannotGeolocate(self.sid, sensing_time)
171
172 self.this_geo = self.Geo(*row)
173 # logging.debug('Retr {geo}'.format(geo=self.this_geo))
174 if self.prev_geo is not None:
175 self.delta_lon = ((self.this_geo.lon + 180.0) - (self.prev_geo.lon + 180.0)) % 360.0
176 if self.delta_lon > 180.0:
177 self.delta_lon -= 360.0
178
179 self.prev_lon = self.prev_geo.lon + 180.0
180
181 self.delta_lat = self.this_geo.lat - self.prev_geo.lat
182
183 self.time_diff = float(timedelta_to_us(self.this_geo.utc - self.prev_geo.utc))
184
185 if self.time_diff is None:
186 raise CannotGeolocate(
187 message='No time diff set', sid=self.sid, sensing_time=sensing_time)
188
189 fract = timedelta_to_us(sensing_time - self.prev_geo.utc) / self.time_diff
190
191 lon = ((self.prev_lon + self.delta_lon * fract) % 360.0) - 180.0
192 lat = self.prev_geo.lat + self.delta_lat * fract
193
194 # if lat > 180 or lat < -180:
195 # logging.warn('lat of {lat} at {time}'.format(lat=lat, time=sensing_time))
196
197 return lat, lon
198
199 def x_y(self, sensing_time):
200 """Compute x, y point for `sensing_time`.
201
202 Class must have been initialised with the optional `width` and `height` parameters.
203 """
204
205 while self.this_geo is None or self.prev_geo is None or sensing_time >= self.this_geo.utc:
206 self.prev_geo = self.this_geo
207 row = self.cur.fetchone()
208 if row is None:
209 raise CannotGeolocate(self.sid, sensing_time)
210
211 self.this_geo = self.Geo(*row)
212 if self.prev_geo is not None:
213 self.delta_lon = ((self.this_geo.lon + 180.0) - (self.prev_geo.lon + 180.0)) % 360.0
214 if self.delta_lon > 180.0:
215 self.delta_lon -= 360.0
216
217 self.prev_lon = self.prev_geo.lon + 180.0
218
219 self.delta_lat = self.this_geo.lat - self.prev_geo.lat
220
221 self.prev_lat = self.prev_geo.lat
222
223 self.time_diff = float(timedelta_to_us(self.this_geo.utc - self.prev_geo.utc))
224
225 fract = timedelta_to_us(sensing_time - self.prev_geo.utc) / self.time_diff
226
227 return (((self.prev_lon + self.delta_lon * fract) % 360.0) * self.width_mul,
228 self.height - ((self.prev_lat + self.delta_lat * fract) + 90.0) * self.height_mul)
229
230 def pos_vel(self, sensing_time):
231 """Return 3d position and vector for `sensing_time`.
232
233 Only valid if `pos_vel` was set to True in constructor.
234 """
235 while self.this_posvel is None or self.prev_posvel is None or sensing_time >= self.this_utc:
236 self.prev_utc = self.this_utc
237 self.prev_posvel = self.this_posvel
238 row = self.cur.fetchone()
239 if row is None:
240 raise CannotGeolocate(self.sid, sensing_time)
241
242 self.this_utc = row[0]
243 self.this_posvel = np.array(row[1:7], dtype=float)
244
245 if self.prev_utc is not None:
246 self.time_diff = float(timedelta_to_us(self.this_utc - self.prev_utc))
247
248 fract = timedelta_to_us(sensing_time - self.prev_utc) / self.time_diff
249
250 res = self.prev_posvel + (self.this_posvel - self.prev_posvel) * fract
251
252 return res[0:3], res[3:6]