1#!/usr/bin/env python3
2
3"""Use EPSSG ORBITREST containing earth fixed-frames vectors to geolocate timeseries data."""
4
5import logging
6import math
7
8from datetime import timedelta
9from chart.db import ts
10from chart.db.model.table import TableInfo
11from chart.plots.geoloc import CannotGeolocate
12
13logger = logging.getLogger()
14
15# There is an entry in geo_location table every 30 secs.
16# make the window + and - 15 secs to get at least 1 entry
17GEOLOCATION_WINDOW = 15
18DB_RETRIEVE_WINDOW = timedelta(seconds=GEOLOCATION_WINDOW)
19
20# There exist several fields with Lat/Lon information for EPSSG.
21# These fields are in X and Y based on earth fixed-frames.
22GEO_LOC_TABLE = TableInfo("GEO_LOCATION")
23X_FIELD_NAME = "X"
24Y_FIELD_NAME = "Y"
25Z_FIELD_NAME = "Z"
26
27
28def ecef2lla(x, y, z):
29 """Convert ECEF x,y,z coordinates to latitude, longitude and altitude
30
31 - input: (x,y,z)
32 - output: (lat,lon,alt)
33
34 Converted to python from MATLAB code by Michael Kleder, July 2005
35 """
36
37 a = 6378137
38 e = 8.1819190842622e-2
39
40 x /= 100.0
41 y /= 100.0
42 z /= 100.0
43 # calculations:
44 try:
45 b = math.sqrt(a ** 2 * (1 - e ** 2))
46 ep = math.sqrt((a ** 2 - b ** 2) / b ** 2)
47 p = math.sqrt(x ** 2 + y ** 2)
48 th = math.atan2(a * z, b * p)
49 lon = math.atan2(y, x)
50 lat = math.atan2(
51 (z + ep ** 2 * b * math.sin(th) ** 3), (p - e ** 2 * a * math.cos(th) ** 3)
52 )
53 N = a / math.sqrt(1 - e ** 2 * math.sin(lat) ** 2)
54 alt = p / math.cos(lat) - N
55 # return lon in range [0,2*pi)
56 lon = lon % (2 * math.pi)
57 lon = lon * 360.0 / (2 * math.pi) - 180.0
58 lat = lat * 360.0 / (2 * math.pi)
59
60 except (ValueError, ZeroDivisionError):
61 lat = 0.0
62 lon = 0.0
63 alt = 0.0
64
65 return lat, lon, alt
66
67
68class Geoloc:
69 """Create an object which can be used to geolocate one or more timestamps for a satellite.
70
71 Initialised with a spacecraft ID and general time range,
72 then the individual functions can be called (must be called in time order)
73 to retrieve geolocations for specific time values.
74 """
75
76 def __init__(self, sid, sensing_start, sensing_stop, width=None, height=None):
77 """Create a geolocator object.
78
79 `sid` (str) :: Source ID
80 `sensing_start` (datetime) :: Minimum validity time for this geolocator
81 `sensing_stop` (datetime) :: Maximum validity time for this geolocator
82 `width` (int) :: Width of target output canvas if x_y() function is to be used
83 `height` (int) :: Height of target output canvas if x_y() function is to be used
84 """
85 self.sid = sid
86 self.sensing_start = sensing_start
87 self.sensing_stop = sensing_stop
88 self.width = width
89 self.height = height
90
91 def lat_lon(self, sensing_time):
92 """Compute latitude, longitude for `sensing_time`.
93
94 Obtain the data from within a validity window in case concrete
95 timestamp does not exist in the DB due to gaps, etc. If no data exist
96 within this window, raise an exception (CannotGeolocate).
97
98 Returns:
99 Tuple of (lat, lon) in degrees
100 Raises:
101 CannotGeolocate: Failed to perform geo location
102 """
103
104 # get a window of data around the desired sensing_time
105 sensing_start = sensing_time - DB_RETRIEVE_WINDOW
106 sensing_stop = sensing_time + DB_RETRIEVE_WINDOW
107 cur = ts.select(
108 GEO_LOC_TABLE,
109 fields=(X_FIELD_NAME, Y_FIELD_NAME, Z_FIELD_NAME),
110 sensing_start=sensing_start,
111 sensing_stop=sensing_stop,
112 sid=self.sid,
113 )
114 # should switch to db.ts.prepare_select() for speed
115 # but that means updating SID classes to handle timeseries prepared cursors
116 # Or switch to single query for the entire life of the Geoloc object
117
118 rs_dirty = cur.fetchall()
119 rs = []
120 rowcount = 0
121
122 # remove empty rows
123 for rs_i in rs_dirty:
124 if rs_i[0] is not None and rs_i[1] is not None:
125 rowcount += 1
126 # this is where we need to convert rs_i to lon_lat vector
127 lat, lon, _ = ecef2lla(rs_i[0], rs_i[1], rs_i[2])
128 rs.append(tuple((lat, lon)))
129
130 # return data in the middle of the retrieved interval
131 if rowcount == 0:
132 raise CannotGeolocate(self.sid, sensing_time)
133
134 if rowcount == 1:
135 return rs[0]
136
137 else:
138 return rs[(rowcount // 2) - 1]
139
140 def x_y(self, sensing_time):
141 """Convert lat/lon to a cylindrical projection, output in pixels."""
142
143 lat, lon = self.lat_lon(sensing_time)
144 return ((lon / 360) + 0.5) * self.width, ((-lat / 180) + 0.5) * self.height