From 3f909e2310f229de88e14af2c511471d95738e39 Mon Sep 17 00:00:00 2001
From: "Mads M. Pedersen" <mmpe@dtu.dk>
Date: Wed, 29 Mar 2017 10:20:06 +0200
Subject: [PATCH] Power shear nan handling + test

---
 wetb/wind/shear.py            |  5 ++++-
 wetb/wind/tests/test_Shear.py | 18 ++++++++++++++++++
 2 files changed, 22 insertions(+), 1 deletion(-)

diff --git a/wetb/wind/shear.py b/wetb/wind/shear.py
index dd61224..dd3281b 100644
--- a/wetb/wind/shear.py
+++ b/wetb/wind/shear.py
@@ -99,7 +99,7 @@ def fit_power_shear_ref(z_u_lst, z_ref, plt=None):
     """
     def shear_error(x, z_u_lst, z_ref):
         alpha, u_ref = x
-        return np.sum([(u - u_ref * (z / z_ref) ** alpha) ** 2 for z, u in z_u_lst])
+        return np.nansum([(u - u_ref * (z / z_ref) ** alpha) ** 2 for z, u in z_u_lst])
     z_u_lst = [(z, np.mean(u)) for z, u in z_u_lst]
     alpha, u_ref = fmin(shear_error, (.1, 10), (z_u_lst, z_ref), disp=False)
     if plt:
@@ -107,6 +107,9 @@ def fit_power_shear_ref(z_u_lst, z_ref, plt=None):
         plt.plot(u, z, '.')
         z = np.linspace(min(z), max(z), 100)
         plt.plot(power_shear(alpha, z_ref, u_ref)(z), z)
+        plt.margins(.1)
+    if alpha==.1 and u_ref==10: # Initial conditions
+        return np.nan, np.nan
     return alpha, u_ref
 
 
diff --git a/wetb/wind/tests/test_Shear.py b/wetb/wind/tests/test_Shear.py
index 38f1d0d..44a369a 100644
--- a/wetb/wind/tests/test_Shear.py
+++ b/wetb/wind/tests/test_Shear.py
@@ -49,6 +49,24 @@ class TestShear(unittest.TestCase):
         14: WSP gl. coo.,Vz 21
         """
 
+    def test_power_shear_fit(self):
+        z = [30,50,70]
+        a,u_ref = .3,10
+        u = power_shear(a,50,u_ref)(z)
+        np.testing.assert_array_almost_equal(fit_power_shear_ref(zip(z,u), 50), [a,u_ref],3)
+        if 0:
+            import matplotlib.pyplot as plt
+            fit_power_shear_ref(zip(z,u), 50, plt)
+            plt.show()
+        
+    def test_power_shear_fit_nan(self):
+        z = [30,50,70]
+        a,u_ref = .3,10
+        u = power_shear(a,50,u_ref)(z)
+        u[2] = np.nan
+        np.testing.assert_array_almost_equal(fit_power_shear_ref(zip(z,u), 50), [a,u_ref],3)
+        u[:] = np.nan
+        np.testing.assert_array_almost_equal(fit_power_shear_ref(zip(z,u), 50), [np.nan, np.nan],3)
 
     def test_power_shear(self):
         if all:
-- 
GitLab