Coverage for services/inference/src/utils/uncertainty.py: 84%

1"""Uncertainty ellipse calculations for Phase 6 Inference Service.""" 

2import numpy as np 

3from typing import Dict 

4import math 

5import logging 

6 

7logger = logging.getLogger(__name__) 

8 

9 

10def compute_uncertainty_ellipse( 

11 sigma_x: float, 

12 sigma_y: float, 

13 covariance_xy: float = 0.0, 

14 confidence_level: float = 0.68, # 1-sigma 

15) -> Dict: 

16 """ 

17 Convert (sigma_x, sigma_y) to uncertainty ellipse parameters. 

18  

19 Computes semi-major and semi-minor axes and rotation angle 

20 from the covariance matrix eigenvalues/eigenvectors. 

21  

22 Args: 

23 sigma_x: Standard deviation in X direction (meters) 

24 sigma_y: Standard deviation in Y direction (meters) 

25 covariance_xy: Covariance between X and Y (default 0 = uncorrelated) 

26 confidence_level: Confidence interval 

27 0.68 = 1-sigma (68%) 

28 0.95 = 2-sigma (95%) 

29 0.997 = 3-sigma (99.7%) 

30  

31 Returns: 

32 Dict with ellipse parameters: 

33 { 

34 "semi_major_axis": float, # Meters 

35 "semi_minor_axis": float, # Meters 

36 "rotation_angle": float, # Degrees (-180 to 180) 

37 "confidence_interval": float, 

38 "area_m2": float, # Area of ellipse 

39 "semi_major_scaled": float, # Scaled for confidence level 

40 "semi_minor_scaled": float, # Scaled for confidence level 

41 } 

42 """ 

43 try: 

44 # Create covariance matrix 

45 cov_matrix = np.array([ 

46 [sigma_x**2, covariance_xy], 

47 [covariance_xy, sigma_y**2], 

48 ]) 

49 

50 # Compute eigenvalues and eigenvectors 

51 eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix) 

52 

53 # Sort eigenvalues in descending order (largest first) 

54 idx = eigenvalues.argsort()[::-1] 

55 eigenvalues = eigenvalues[idx] 

56 eigenvectors = eigenvectors[:, idx] 

57 

58 # Semi-major and semi-minor axes (1-sigma) 

59 semi_major = math.sqrt(max(eigenvalues[0], 0.0)) # Avoid negative due to numerical errors 

60 semi_minor = math.sqrt(max(eigenvalues[1], 0.0)) 

61 

62 # Rotation angle (eigenvector orientation) 

63 rotation_rad = math.atan2(eigenvectors[1, 0], eigenvectors[0, 0]) 

64 rotation_deg = math.degrees(rotation_rad) 

65 

66 # Normalize rotation to [-180, 180] 

67 while rotation_deg > 180: 

68 rotation_deg -= 360 

69 while rotation_deg < -180: 

70 rotation_deg += 360 

71 

72 # Scale for confidence level 

73 # For 1-sigma: scale = 1.0 

74 # For 2-sigma: scale = 2.0 

75 # For 95%: scale ≈ 2.45 (chi-square distribution for 2D) 

76 if confidence_level >= 0.95: 

77 scale_factor = 2.45 # 2-sigma ≈ 95% for Gaussian 

78 elif confidence_level >= 0.90: 

79 scale_factor = 2.15 

80 else: 

81 scale_factor = 1.0 

82 

83 result = { 

84 "semi_major_axis": float(semi_major), 

85 "semi_minor_axis": float(semi_minor), 

86 "rotation_angle": float(rotation_deg), 

87 "confidence_interval": float(confidence_level), 

88 "area_m2": float(math.pi * semi_major * semi_minor), 

89 "semi_major_scaled": float(semi_major * scale_factor), 

90 "semi_minor_scaled": float(semi_minor * scale_factor), 

91 } 

92 

93 logger.debug(f"Uncertainty ellipse computed: {result}") 

94 return result 

95 

96 except Exception as e: 

97 logger.error(f"Failed to compute uncertainty ellipse: {e}", exc_info=True) 

98 raise 

99 

100 

101def ellipse_to_geojson( 

102 center_lat: float, 

103 center_lon: float, 

104 semi_major_m: float, 

105 semi_minor_m: float, 

106 rotation_deg: float, 

107 num_points: int = 64, 

108) -> Dict: 

109 """ 

110 Convert uncertainty ellipse to GeoJSON Feature for Mapbox visualization. 

111  

112 Generates a polygon approximating the ellipse for rendering on map. 

113 Handles the conversion from meters to geographic coordinates (lat/lon). 

114  

115 Args: 

116 center_lat: Center latitude in decimal degrees (-90 to 90) 

117 center_lon: Center longitude in decimal degrees (-180 to 180) 

118 semi_major_m: Semi-major axis in meters 

119 semi_minor_m: Semi-minor axis in meters 

120 rotation_deg: Rotation angle in degrees 

121 num_points: Number of points to approximate ellipse (default 64 for smooth curve) 

122  

123 Returns: 

124 GeoJSON Feature dict with Polygon geometry: 

125 { 

126 "type": "Feature", 

127 "geometry": { 

128 "type": "Polygon", 

129 "coordinates": [[[lon, lat], [lon, lat], ...]] 

130 }, 

131 "properties": { 

132 "name": "Uncertainty Ellipse", 

133 "semi_major_m": float, 

134 "semi_minor_m": float, 

135 "rotation_deg": float, 

136 } 

137 } 

138 """ 

139 try: 

140 # Earth radius in meters (WGS84) 

141 EARTH_RADIUS_M = 6371000.0 

142 

143 # Convert meter-based ellipse to degrees 

144 # At equator: 1 degree ≈ 111km 

145 lat_scale = semi_major_m / EARTH_RADIUS_M 

146 lon_scale = semi_minor_m / (EARTH_RADIUS_M * math.cos(math.radians(center_lat))) 

147 

148 # Generate ellipse points 

149 angles = np.linspace(0, 2 * np.pi, num_points, endpoint=False) 

150 rotation_rad = math.radians(rotation_deg) 

151 

152 coordinates = [] 

153 for angle in angles: 

154 # Ellipse parametric equations 

155 x = semi_major_m * math.cos(angle) 

156 y = semi_minor_m * math.sin(angle) 

157 

158 # Apply rotation 

159 x_rot = x * math.cos(rotation_rad) - y * math.sin(rotation_rad) 

160 y_rot = x * math.sin(rotation_rad) + y * math.cos(rotation_rad) 

161 

162 # Convert to lat/lon 

163 lat = center_lat + y_rot / EARTH_RADIUS_M 

164 lon = center_lon + x_rot / (EARTH_RADIUS_M * math.cos(math.radians(center_lat))) 

165 

166 coordinates.append([lon, lat]) 

167 

168 # Close polygon (first point == last point) 

169 coordinates.append(coordinates[0]) 

170 

171 result = { 

172 "type": "Feature", 

173 "geometry": { 

174 "type": "Polygon", 

175 "coordinates": [coordinates], 

176 }, 

177 "properties": { 

178 "name": "Uncertainty Ellipse", 

179 "semi_major_m": float(semi_major_m), 

180 "semi_minor_m": float(semi_minor_m), 

181 "rotation_deg": float(rotation_deg), 

182 "center_lat": float(center_lat), 

183 "center_lon": float(center_lon), 

184 }, 

185 } 

186 

187 logger.debug(f"GeoJSON ellipse generated with {len(coordinates)} points") 

188 return result 

189 

190 except Exception as e: 

191 logger.error(f"Failed to convert ellipse to GeoJSON: {e}", exc_info=True) 

192 raise 

193 

194 

195def create_uncertainty_circle( 

196 center_lat: float, 

197 center_lon: float, 

198 radius_m: float, 

199 num_points: int = 32, 

200) -> Dict: 

201 """ 

202 Create a circle (special case of ellipse) as GeoJSON. 

203  

204 Useful for simpler visualization when sigma_x ≈ sigma_y. 

205  

206 Args: 

207 center_lat: Center latitude 

208 center_lon: Center longitude 

209 radius_m: Radius in meters 

210 num_points: Number of points (default 32) 

211  

212 Returns: 

213 GeoJSON Feature with circular polygon 

214 """ 

215 return ellipse_to_geojson( 

216 center_lat, center_lon, radius_m, radius_m, 0.0, num_points 

217 )