From 9a461f5947ea3258d957598d1d90a4600f72986b Mon Sep 17 00:00:00 2001
From: hgloeckner <helene.gloeckner@posteo.de>
Date: Thu, 23 Jan 2025 12:03:10 +0100
Subject: [PATCH] copy error estimates from JOANNE

---
 pydropsonde/circles.py | 71 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 71 insertions(+)

diff --git a/pydropsonde/circles.py b/pydropsonde/circles.py
index 0b53a1b4..04a41940 100644
--- a/pydropsonde/circles.py
+++ b/pydropsonde/circles.py
@@ -4,6 +4,7 @@
 import tqdm
 import circle_fit as cf
 import pydropsonde.helper.physics as hp
+import pydropsonde.helper.xarray_helper as hx
 
 _no_default = object()
 
@@ -230,6 +231,76 @@ def apply_fit2d(self, variables=None):
         self.circle_ds = ds
         return self
 
+    def add_regression_stderr(self, variables=None):
+        """
+        Calculation of regression standard error, following Lenschow, Donald H and Savic-Jovcic,Verica and Stevens, Bjorn 2007
+
+        """
+        alt_dim = self.alt_dim
+        if variables is None:
+            variables = ["u", "v", "q", "ta", "p", "density"]
+        ds = self.circle_ds
+
+        dx_denominator = ((ds.x - ds.x.mean(dim="sonde_id")) ** 2).sum(dim="sonde_id")
+        dy_denominator = ((ds.y - ds.y.mean(dim="sonde_id")) ** 2).sum(dim="sonde_id")
+
+        for var in variables:
+            var_err = (
+                ds[var]
+                - (ds[f"mean_{var}"] + ds[f"d{var}dx"] * ds.x + ds[f"d{var}dy"] * ds.y)
+            ) ** 2
+            nominator = (var_err) / (var_err.count(dim="sonde_id") - 3)
+
+            se_x = np.sqrt(nominator / dx_denominator)
+            se_y = np.sqrt(nominator / dy_denominator)
+
+            dvardx_std_name = ds[f"d{var}dx"].attrs.get(
+                "standard_name", f"derivative_of_{var}_wrt_x"
+            )
+            unit = ds[f"d{var}dx"].attrs.get("units", "")
+            dvardy_std_name = ds[f"d{var}dy"].attrs.get(
+                "standard_name", f"derivative_of_{var}_wrt_y"
+            )
+
+            ds = ds.assign(
+                {
+                    f"se_d{var}dx": (
+                        ["sonde_id", alt_dim],
+                        se_x.values,
+                        dict(
+                            standard_name=f"{dvardx_std_name} standard_error",
+                            units=unit,
+                        ),
+                    ),
+                    f"se_d{var}dy": (
+                        ["sonde_id", alt_dim],
+                        se_y.values,
+                        dict(
+                            standard_name=f"{dvardy_std_name} standard_error",
+                            units=unit,
+                        ),
+                    ),
+                }
+            )
+            ds = hx.add_ancillary_var(ds, f"d{var}dx", f"se_d{var}dx")
+            ds = hx.add_ancillary_var(ds, f"d{var}dy", f"se_d{var}dy")
+
+        div_std_name = ds["div"].attrs.get("standard_name", "divergence_of_wind")
+        ds = ds.assign(
+            {
+                "se_div": (
+                    ["sonde_id", alt_dim],
+                    np.sqrt(ds.se_dudx**2 + ds.se_dudy**2).values,
+                    dict(standard_name=f"{div_std_name} standard_error", units=unit),
+                )
+            }
+        )
+
+        ds = hx.add_ancillary_var(ds, "div", "se_div")
+        ds = hx.add_ancillary_var(ds, "vor", "se_div")
+
+        return self
+
     def add_density(self, sonde_dim="sonde_id", alt_dim="gpsalt"):
         """
         Calculate and add the density to the circle dataset.