bugfix for #234 (#235)

src/inversion_scripts/invert.py

@@ -150,7 +150,10 @@ def do_inversion(
                 for p in obs_GC[:, 4]
             ]
         )
+        # Define observational errors (diagonal entries of S_o matrix)
+        obs_error = np.power(obs_error, 2)
         gP = s_superO_p**2 / s_superO_1**2
+        # scale error variance by gP
         obs_error = gP * obs_err
 
         # check to make sure obs_err isn't negative, set 1 as default value
@@ -179,9 +182,6 @@ def do_inversion(
             scaling_matrix = np.tile(scale_factors, (reps, 1))
             K *= scaling_matrix
 
-        # Define observational errors (diagonal entries of S_o matrix)
-        obs_error = np.power(obs_error, 2)
-
         # Measurement-model mismatch: TROPOMI columns minus GEOS-Chem virtual TROPOMI columns
         # This is (y - F(xA)), i.e., (y - (K*xA + c)) or (y - K*xA) in shorthand
         delta_y = obs_GC[:, 0] - obs_GC[:, 1]  # [ppb]

src/inversion_scripts/lognormal_invert.py

@@ -82,7 +82,7 @@ def lognormal_invert(config, state_vector_filepath, jacobian_sf):
 
     # get the So matrix
     ds = np.load("so_super.npz")
-    so = ds["so"] ** 2
+    so = ds["so"]
 
     # Calculate the difference between tropomi and the background
     # simulation, which has no emissions

src/inversion_scripts/merge_partial_k.py

@@ -86,7 +86,9 @@ def merge_partial_k(satdat_dir, lat_bounds, lon_bounds, obs_err, precomp_K):
                 for p in obs_GC[:, 4]
             ]
         )
+        # scale error variance by gP value following Chen et al. 2023
         gP = s_superO_p**2 / s_superO_1**2
+        obs_error = obs_error ** 2
         obs_error = gP * obs_error
 
         # check to make sure obs_err isn't negative, set 1 as default value