From ece8796ee84f41a81b4066c2f81315c2ffba7694 Mon Sep 17 00:00:00 2001 From: Jone Peter Reistad Date: Thu, 18 Jan 2024 20:44:27 +0100 Subject: [PATCH] allow error for each component of data vector to be specified (#33) --- lompe/model/data.py | 13 ++++++++++--- lompe/model/model.py | 7 +++++-- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/lompe/model/data.py b/lompe/model/data.py index e889794..28fb2f5 100644 --- a/lompe/model/data.py +++ b/lompe/model/data.py @@ -124,8 +124,12 @@ def __init__(self, values, coordinates = None, LOS = None, components = 'all', d magnetometer data and iweight=1.0 for ionospheric convection measurements. Keep in mind that this weight is directly applied to the a priori inverse data covariance matrix, so the data error is effectively increased by a factor of 1/sqrt(iweight). - error: array of same length as values, or float, optional + error: either float, or 1D array of same length as values, or 2D array + (shape 2,N or 3,N), optional Measurement error. Used to calculate the data covariance matrix. Use SI units. + If 2D array, it means that the error of each component (east, north, up) is + specified separately. If 1D array, the same value is used for each component. + If float, the same value is used for all observations and components. """ @@ -214,8 +218,11 @@ def subset(self, indices): """ self.values = np.array(self.values , ndmin = 2)[:, indices].squeeze() - self.error = self.error[indices] - + if np.array(self.error, ndmin=2).shape[0] >= 2: # errors for each component + self.error = self.error[:,indices] + else: + self.error = self.error[indices] + for key in self.coords.keys(): self.coords[key] = self.coords[key][indices] diff --git a/lompe/model/model.py b/lompe/model/model.py index 33d38b8..090fccd 100644 --- a/lompe/model/model.py +++ b/lompe/model/model.py @@ -312,8 +312,11 @@ def run_inversion(self, l1 = 0, l2 = 0, spatial_weight = np.ones(ds.values.size) dimensions = np.array(ds.values, ndmin = 2).shape[0] - error = np.tile(ds.error, dimensions) - + if np.array(ds.error, ndmin=2).shape[0]==1: + error = np.tile(ds.error, dimensions) + else: #error is different for different components + error = ds.error.flatten() + w_i = spatial_weight * 1/(error**2) * iweights[ii] if iweights[ii] != 1: print('{}: Measurement uncertainty effectively changed from {} to {}'.format(dtype, np.median(error), np.median(error)/np.sqrt(iweights[ii])))