Merge pull request #39 from DHI-GRAS/32bit

Use 32 bit floats for most calculations

Author: hectornieto

pyTSEB/PyTSEB.py

@@ -172,7 +172,7 @@ def process_local_image(self):
             elif field == "input_mask":
                 if self.p['input_mask'] == '0':
                     # Create mask from landcover array
-                    mask = np.ones(dims)
+                    mask = np.ones(dims, np.int32)
                     mask[np.logical_or.reduce((in_data['landcover'] == res.WATER,
                                                in_data['landcover'] == res.URBAN,
                                                in_data['landcover'] == res.SNOW))] = 0
@@ -530,12 +530,12 @@ def run(self, in_data, mask=None):
                         "resistance_form": [self.resistance_form, self.res_params]}
 
         if mask is None:
-            mask = np.ones(in_data['LAI'].shape)
+            mask = np.ones(in_data['LAI'].shape, np.int32)
 
         # Create the output dictionary
         out_data = dict()
         for field in self._get_output_structure():
-            out_data[field] = np.zeros(in_data['LAI'].shape) + np.NaN
+            out_data[field] = np.zeros(in_data['LAI'].shape, np.float32) + np.NaN
 
         # Esimate diffuse and direct irradiance
         difvis, difnir, fvis, fnir = rad.calc_difuse_ratio(
@@ -591,10 +591,11 @@ def run(self, in_data, mask=None):
                                f_c=in_data['f_c'][i])
 
         # Net shortwave radiation for vegetation
-        F = np.zeros(in_data['LAI'].shape)
+        F = np.zeros(in_data['LAI'].shape, np.float32)
         F[i] = in_data['LAI'][i] / in_data['f_c'][i]
         # Clumping index
-        omega0, Omega = np.zeros(in_data['LAI'].shape), np.zeros(in_data['LAI'].shape)
+        omega0 = np.zeros(in_data['LAI'].shape, np.float32)
+        Omega = np.zeros(in_data['LAI'].shape, np.float32)
         omega0[i] = CI.calc_omega0_Kustas(
             in_data['LAI'][i],
             in_data['f_c'][i],
@@ -637,7 +638,7 @@ def run(self, in_data, mask=None):
 
         if self.water_stress:
             i = np.array(np.logical_and(~noVegPixels, mask == 1))
-            [_, _, _, _, _, _, out_data['LE_0'][i], _, 
+            [_, _, _, _, _, _, out_data['LE_0'][i], _,
              out_data['LE_C_0'][i], _, _, _, _, _, _, _, _, _, _] = \
                  pet.shuttleworth_wallace(
                               in_data['T_A1'][i],
@@ -670,13 +671,11 @@ def run(self, in_data, mask=None):
 
             out_data['CWSI'][i] = 1.0 - (out_data['LE_C1'][i] / out_data['LE_C_0'][i])
 
-        
         if self.calc_daily_ET:
-            out_data['ET_day'] = met.flux_2_evaporation(in_data['S_dn_24'] * out_data['LE1'] / in_data['S_dn'], 
-                                                        t_k=20+273.15, 
+            out_data['ET_day'] = met.flux_2_evaporation(in_data['S_dn_24'] * out_data['LE1'] / in_data['S_dn'],
+                                                        t_k=20+273.15,
                                                         time_domain=24)
-            
-        
+
         print("Finished processing!")
         return out_data
 
@@ -801,7 +800,7 @@ def _set_param_array(self, parameter, dims, band=1):
 
         # See if the parameter is a number
         try:
-            array = np.zeros(dims) + float(parameter)
+            array = np.zeros(dims, np.float32) + float(parameter)
             return success, array
         except ValueError:
             pass
@@ -814,19 +813,19 @@ def _set_param_array(self, parameter, dims, band=1):
             return success, array
         # If it is then get the value of that parameter
         try:
-            array = np.zeros(dims) + float(inputString)
+            array = np.zeros(dims, np.float32) + float(inputString)
         except ValueError:
             try:
                 fid = gdal.Open(inputString, gdal.GA_ReadOnly)
                 if self.subset:
                     array = fid.GetRasterBand(band).ReadAsArray(self.subset[0],
                                                                 self.subset[1],
                                                                 self.subset[2],
-                                                                self.subset[3])
+                                                                self.subset[3]).astype(np.float32)
                 else:
-                    array = fid.GetRasterBand(band).ReadAsArray()
+                    array = fid.GetRasterBand(band).ReadAsArray().astype(np.float32)
             except AttributeError:
-                print("%s image not present for parameter %s"%(inputString, parameter))
+                print("%s image not present for parameter %s" % (inputString, parameter))
                 success = False
             finally:
                 fid = None
@@ -1065,7 +1064,7 @@ def _get_input_structure(self):
                             # Soil heat flux parameter
                             ("G", "Soil Heat Flux Parameter"),
                             ('S_dn_24', 'Daily shortwave irradiance')])
-            
+
         return input_fields
 
     def _set_special_model_input(self, field, dims):
@@ -1577,9 +1576,9 @@ def _set_special_model_input(self, field, dims):
                 inputs[field] = fid.GetRasterBand(1).ReadAsArray(subset[0],
                                                                  subset[1],
                                                                  subset[2],
-                                                                 subset[3])
+                                                                 subset[3]).astype(np.float32)
             else:
-                inputs[field] = fid.GetRasterBand(1).ReadAsArray()
+                inputs[field] = fid.GetRasterBand(1).ReadAsArray().astype(np.float32)
             inputs['scale'] = [geo_LR, prj_LR, self.geo, self.prj]
             success = True
         else:

pyTSEB/TSEB.py

@@ -318,7 +318,7 @@ def TSEB_2T(T_C,
     # calcG_params[1] = None
     # Create the output variables
     [flag, Ln_S, Ln_C, LE_C, H_C, LE_S, H_S, G, R_S, R_x,
-        R_A, iterations] = [np.zeros(T_S.shape)+np.NaN for i in range(12)]
+        R_A, iterations] = [np.zeros(T_S.shape, np.float32)+np.NaN for i in range(12)]
     T_AC = T_A_K.copy()
 
     # iteration of the Monin-Obukhov length
@@ -675,10 +675,10 @@ def TSEB_PT(Tr_K,
     # iteration of the Monin-Obukhov length
     if const_L is None:
         # Initially assume stable atmospheric conditions and set variables for
-        L = np.asarray(np.zeros(Tr_K.shape, np.float32) + np.inf, dtype=np.float32)
+        L = np.zeros(Tr_K.shape) + np.inf
         max_iterations = ITERATIONS
     else:  # We force Monin-Obukhov lenght to the provided array/value
-        L = np.asarray(np.ones(Tr_K.shape, np.float32) * const_L, dtype=np.float32)
+        L = np.ones(Tr_K.shape) * const_L
         max_iterations = 1  # No iteration
     # Calculate the general parameters
     rho = met.calc_rho(p, ea, T_A_K)  # Air density
@@ -696,7 +696,7 @@ def TSEB_PT(Tr_K,
     u_friction = MO.calc_u_star(u, z_u, L, d_0, z_0M)
     u_friction = np.asarray(np.maximum(U_FRICTION_MIN, u_friction), dtype=np.float32)
     L_queue = deque([np.array(L, np.float32)], 6)
-    L_converged = np.asarray(np.zeros(Tr_K.shape), dtype=bool)
+    L_converged = np.zeros(Tr_K.shape, bool)
     L_diff_max = np.inf
 
     # First assume that canopy temperature equals the minumum of Air or
@@ -1133,7 +1133,7 @@ def DTD(Tr_K_0,
     resistance_form = resistance_form[0]
     # Create the output variables
     [flag, T_S, T_C, T_AC, Ln_S, Ln_C, LE_C, H_C, LE_S, H_S, G, R_S, R_x,
-        R_A, H, iterations] = [np.zeros(Tr_K_1.shape) + np.NaN for i in range(16)]
+        R_A, H, iterations] = [np.zeros(Tr_K_1.shape, np.float32) + np.NaN for i in range(16)]
 
     # Calculate the general parameters
     rho = met.calc_rho(p, ea, T_A_K_1)  # Air density
@@ -1480,7 +1480,7 @@ def OSEB(Tr_K,
                          calcG_params[1]],
                         [Tr_K] * 12)
     # Create the output variables
-    [flag, Ln, LE, H, G, R_A] = [np.zeros(Tr_K.shape) + np.NaN for i in range(6)]
+    [flag, Ln, LE, H, G, R_A] = [np.zeros(Tr_K.shape, np.float32) + np.NaN for i in range(6)]
 
     # iteration of the Monin-Obukhov length
     if const_L is None:

pyTSEB/energy_combination_ET.py

@@ -154,7 +154,7 @@ def penman_monteith(T_A_K,
                       [T_A_K] * 14)
 
     # Create the output variables
-    [flag, Ln, LE, H, G, R_A, iterations] = [np.zeros(T_A_K.shape) + np.NaN for i in
+    [flag, Ln, LE, H, G, R_A, iterations] = [np.zeros(T_A_K.shape, np.float32) + np.NaN for i in
                                              range(7)]
 
     # Calculate the general parameters
@@ -175,10 +175,10 @@ def penman_monteith(T_A_K,
     # iteration of the Monin-Obukhov length
     if const_L is None:
         # Initially assume stable atmospheric conditions and set variables for
-        L = np.asarray(np.zeros(T_A_K.shape) + np.inf)
+        L = np.zeros(T_A_K.shape) + np.inf
         max_iterations = ITERATIONS
     else:  # We force Monin-Obukhov lenght to the provided array/value
-        L = np.asarray(np.ones(T_A_K.shape) * const_L)
+        L = np.ones(T_A_K.shape) * const_L
         max_iterations = 1  # No iteration
     u_friction = TSEB.MO.calc_u_star(u, z_u, L, d_0, z_0M)
     u_friction = np.asarray(np.maximum(TSEB.U_FRICTION_MIN, u_friction))
@@ -484,8 +484,7 @@ def shuttleworth_wallace(T_A_K,
 
     # Create the output variables
     [flag, vpd_0, LE, H, LE_C, H_C, LE_S, H_S, G, R_S, R_x, R_A,
-     Rn, Rn_C, Rn_S, C_s, C_c, PM_C, PM_S, iterations] = [np.full(T_A_K.shape,
-                                                                  np.NaN)
+     Rn, Rn_C, Rn_S, C_s, C_c, PM_C, PM_S, iterations] = [np.full(T_A_K.shape, np.NaN, np.float32)
                                                           for i in range(20)]
 
     # Calculate the general parameters
@@ -511,10 +510,10 @@ def shuttleworth_wallace(T_A_K,
     # iteration of the Monin-Obukhov length
     if const_L is None:
         # Initially assume stable atmospheric conditions and set variables for
-        L = np.asarray(np.zeros(T_A_K.shape) + np.inf)
+        L = np.zeros(T_A_K.shape) + np.inf
         max_iterations = ITERATIONS
     else:  # We force Monin-Obukhov lenght to the provided array/value
-        L = np.asarray(np.ones(T_A_K.shape) * const_L)
+        L = np.ones(T_A_K.shape) * const_L
         max_iterations = 1  # No iteration
     u_friction = TSEB.MO.calc_u_star(u, z_u, L, d_0, z_0M)
     u_friction = np.asarray(np.maximum(TSEB.U_FRICTION_MIN, u_friction))
@@ -871,7 +870,7 @@ def penman(T_A_K,
                       [T_A_K] * 11)
 
     # Create the output variables
-    [flag, Ln, LE, H, G, R_A, iterations] = [np.zeros(T_A_K.shape) + np.NaN for i in
+    [flag, Ln, LE, H, G, R_A, iterations] = [np.zeros(T_A_K.shape, np.float32) + np.NaN for i in
                                              range(7)]
 
     # Calculate the general parameters
@@ -889,10 +888,10 @@ def penman(T_A_K,
     # iteration of the Monin-Obukhov length
     if const_L is None:
         # Initially assume stable atmospheric conditions and set variables for
-        L = np.asarray(np.zeros(T_A_K.shape) + np.inf)
+        L = np.zeros(T_A_K.shape) + np.inf
         max_iterations = ITERATIONS
     else:  # We force Monin-Obukhov lenght to the provided array/value
-        L = np.asarray(np.ones(T_A_K.shape) * const_L)
+        L = np.ones(T_A_K.shape) * const_L
         max_iterations = 1  # No iteration
     u_friction = TSEB.MO.calc_u_star(u, z_u, L, d_0, z_0M)
     u_friction = np.asarray(np.maximum(TSEB.U_FRICTION_MIN, u_friction))