Merge pull request #31 from tammasloughran/30-heatwave-identification…

…-misses-negative-temperature-exceedances

Fixed a potential bug with identifying heatwaves that have negative temperatures

ehfheatwaves/__init__.py

@@ -4,5 +4,5 @@
 
 from ehfheatwaves.ehfheatwaves import *
 
-__version__ = '1.1.5'
+__version__ = '1.2.0'
 

ehfheatwaves/ehfheatwaves.py

@@ -106,8 +106,14 @@ def identify_hw(ehfs:np.ndarray)->tuple:
     """
     # Agregate consecutive days with EHF>0
     # First day contains duration
-    events = (ehfs>0.).astype(int)
-    events[events.mask==True] = 0
+    if np.isnan(ehfs).any(): # then ehfs is a view to tmax or tmin
+        # This is handled differently to EHFs because tmax or tmin could theoretically hold -ve
+        # values. Therefore we identify values of the exceedence index that are not nan values as
+        # heatwaves.
+        events = np.logical_not(np.isnan(ehfs)).astype(int)
+    else: # ehfs is a view to actual EHF values
+        events = (ehfs>0.0).astype(int)
+        events[events.mask==True] = 0
     for i in range(events.shape[0] - 2, -1, -1):
          events[i,events[i,...]>0] = events[i+1,events[i,...]>0]+1
 
@@ -147,8 +153,14 @@ def identify_semi_hw(ehfs:np.ndarray)->tuple:
     """
     # Agregate consecutive days with EHF>0
     # First day contains duration
-    events = (ehfs>0.).astype(int)
-    events[events.mask==True] = 0
+    if np.isnan(ehfs).any(): # then ehfs is a view to tmax or tmin
+        # This is handled differently to EHFs because tmax or tmin could theoretically hold -ve
+        # values. Therefore we identify values of the exceedence index that are not nan values as
+        # heatwaves.
+        events = np.logical_not(np.isnan(ehfs)).astype(int)
+    else: # ehfs is a view to actual EHF values
+        events = (ehfs>0.0).astype(int)
+        events[events.mask==True] = 0
     for i in range(events.shape[0] - 2, -1, -1):
          events[i,events[i,...]>0] = events[i+1,events[i,...]>0]+1
 
@@ -473,12 +485,14 @@ def main():
                 idoy = i-timedata.daysinyear*int((i+1)/timedata.daysinyear)
                 txexceed[i,...] = tmax[i,...]>txpct[idoy,...]
             txexceed[txexceed>0] = tmax[txexceed>0]
+            txexceed[txexceed==0] = np.nan # need nans to be identified correctly
         if options.keeptmin:
             tnexceed = np.ma.ones(tmin.shape)*np.nan
             for i in range(0, tmin.shape[0]):
                 idoy = i-timedata.daysinyear*int((i+1)/timedata.daysinyear)
                 tnexceed[i,...] = tmin[i,...]>tnpct[idoy,...]
             tnexceed[tnexceed>0] = tmin[tnexceed>0]
+            tnexceed[tnexceed==0] = np.nan # need nans to be identified correctly
 
     # Calculate daily output
     if options.dailyout or options.dailyout: event, ends = identify_hw(EHF)

tests/run_tests.py

@@ -182,13 +182,24 @@ class TestIdentifyHW(unittest.TestCase):
     """Tests for the identify_hw function."""
 
     # Some example EHF index data
-    ehfdata = np.ma.array([-1.3, -0.3, 3.4, 8.5, -0.4, # Not a heatwave
-                        5.6, 10.2, 20.4, -1.4, # a three day heatwave starting on day 6
-                        7.8, 15.5, 16.9, 17.9, 30.2, -3.3]) # a 5 day heatwave starting on day 10
+    ehfdata = np.ma.array([
+            -1.3, -0.3, 3.4, 8.5, -0.4, # Not a heatwave
+            5.6, 10.2, 20.4, -1.4, # a three day heatwave starting on day 6
+            7.8, 15.5, 16.9, 17.9, 30.2, -3.3, # a 5 day heatwave starting on day 10
+            ])
     ehfdata[ehfdata<0] = 0
     known_events = np.ma.array([0,0,0,0,0,1,1,1,0,1,1,1,1,1,0])
     known_ends = np.ma.array([0,0,0,0,0,3,0,0,0,5,0,0,0,0,0])
 
+    # Example data for a heatwave with negative temperatures.
+    tmindata = np.ma.array([
+            np.nan,np.nan,np.nan,np.nan, # does not exceed threshold so is masked out.
+            -2.5,-1.9,-0.2,0.5,2.1, # Heatwaves with negative values.
+            np.nan,np.nan, # does not exceed threshold.
+            ])
+    tmin_known_events = np.ma.array([0,0,0,0,1,1,1,1,1,0,0])
+    tmin_known_ends = np.ma.array([0,0,0,0,5,0,0,0,0,0,0])
+
     def testReturnTupple(self):
         """Should return a tupple containging the event indicator and the durations (numpy.ndarray)."""
         result = ehfheatwaves.identify_hw(self.ehfdata)
@@ -209,6 +220,12 @@ def testShape(self):
         self.assertEqual(events.shape, input_shape)
         self.assertEqual(ends.shape, input_shape)
 
+    def testNegativeTemperatures(self):
+        """Tmin or tmax exceedances with negative values should be identified correctly."""
+        events, ends = ehfheatwaves.identify_hw(self.tmindata)
+        self.assertTrue((events==self.tmin_known_events).all())
+        self.assertTrue((ends==self.tmin_known_ends).all())
+
 
 class TestIdentifySemiHW(unittest.TestCase):
     """Tests for the identify_semi_hw function."""