From 244a0c540cd19af89d31473bcc994ea85c8cfb2a Mon Sep 17 00:00:00 2001 From: Radoslaw Guzinski Date: Wed, 31 Mar 2021 13:14:49 +0200 Subject: [PATCH] Use 32 bit floats for most calculations --- pyTSEB/PyTSEB.py | 37 ++++++++++++++++----------------- pyTSEB/TSEB.py | 12 +++++------ pyTSEB/energy_combination_ET.py | 19 ++++++++--------- 3 files changed, 33 insertions(+), 35 deletions(-) diff --git a/pyTSEB/PyTSEB.py b/pyTSEB/PyTSEB.py index aac3b43..c4594e0 100755 --- a/pyTSEB/PyTSEB.py +++ b/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,7 +813,7 @@ 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) @@ -822,11 +821,11 @@ def _set_param_array(self, parameter, dims, band=1): 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: diff --git a/pyTSEB/TSEB.py b/pyTSEB/TSEB.py index e47f886..80a7580 100644 --- a/pyTSEB/TSEB.py +++ b/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: diff --git a/pyTSEB/energy_combination_ET.py b/pyTSEB/energy_combination_ET.py index 40d1aa3..5bc9a6f 100755 --- a/pyTSEB/energy_combination_ET.py +++ b/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))