import numpy as np from numpy.lib.stride_tricks import sliding_window_view import utils import conf import ancillary_data as c_tools # ############## GROUP 1 TESTS ############## # # 11 micron brightness temperature threshold test def simple_test(rad, threshold, cmin): radshape = rad.shape rad = rad.reshape(np.prod(radshape)) thr = np.array(threshold) confidence = np.ones(rad.shape) bit = np.zeros(rad.shape) if thr[4] == 1: print("simple test running") # the C code has the line below that I don't quite understand the purpose of. # It seems to be setting the bit to 0 if the BT value is greater than the midpoint # # if (m31 >= dobt11[1]) (void) set_bit(13, pxout.testbits); # confidence = utils.conf_test(rad, thr) confidence = conf.conf_test(rad, thr) bit[rad >= thr[1]] = 1 return np.minimum(cmin, confidence.reshape(radshape)), confidence.reshape(radshape), bit.reshape(radshape) def sst_test(rad1, rad2, vza, surf_temp, threshold, cmin): a1 = 1.8860 a2 = 0.9380 a3 = 0.1280 a4 = 1.0940 radshape = rad1.shape b31 = rad1.reshape(np.prod(radshape)) - 273.16 b32 = rad2.reshape(np.prod(radshape)) - 273.16 thr = np.array(threshold) confidence = np.ones(b31.shape) bit = np.zeros(b31.shape) rad_diff = b31 - b32 sstc = surf_temp.reshape(np.prod(radshape)) - 273.16 mu = np.cos(vza.reshape(np.prod(radshape)) * np.pi/180.0) modsst = 273.16 + a1 + a2*b31 + a3*rad_diff*sstc + a4*rad_diff*((1/mu)-1) sfcdif = surf_temp.reshape(np.prod(radshape)) - modsst if thr[4] == 1: print('SST test running') confidence = conf.conf_test(sfcdif, thr) bit[sfcdif < thr[1]] = 1 return np.minimum(cmin, confidence.reshape(radshape)), confidence.reshape(radshape), bit.reshape(radshape) def test_11_12_diff(data, threshold, cmin): radshape = data.M15.shape b31 = data.M15.values.reshape(np.prod(radshape)) b32 = data.M15.values.reshape(np.prod(radshape)) vza = data.sensor_zenith.values.reshape(np.prod(radshape)) thr = np.array(threshold) confidence = np.ones(b31.shape) bit = np.zeros(b31.shape) rad_diff = b31 - b32 # Get secant of viewing zenith angle dtr = np.pi/180 cosvza = np.cos(vza * dtr) schi = np.full(cosvza.shape, 99.0) schi[cosvza > 0.0] = 1.0/cosvza[cosvza > 0.0] # Need to define this in cython btd_thr = c_tools.py_cithr(1, np.array(schi, dtype=np.float32), np.array(b31, dtype=np.float32)) idx = np.nonzero((btd_thr < 0.1) | (np.abs(schi-99.0) < 0.0001)) btd_thr[idx] = thr[0] locut = btd_thr + 0.3*btd_thr hicut = btd_thr - 1.25 corr_thr = np.array([locut, btd_thr, hicut, np.ones(locut.shape)], dtype=np.float) if thr[1] == 1: print('11-12um diff test running') bit[rad_diff < thr[1]] = 1 confidence = conf.conf_test(rad_diff, corr_thr) return np.minimum(cmin, confidence.reshape(radshape)), confidence.reshape(radshape), bit.reshape(radshape) def test_11_4_diff(rad1, rad2, threshold, scene_flags, cmin): radshape = rad1.shape b31 = rad1.reshape(np.prod(radshape)) b20 = rad2.reshape(np.prod(radshape)) thr = np.array(threshold) sunglint = scene_flags['sunglint'].reshape(np.prod(radshape)) confidence = np.ones(b31.shape) if thr[4] == 1: print('11-4um diff test running') confidence[sunglint == 0] = conf.conf_test((b31-b20)[sunglint == 0], thr) return np.minimum(cmin, confidence.reshape(radshape)), confidence.reshape(radshape) def vis_nir_ratio_test(rad1, rad2, threshold, scene, cmin): if threshold['Daytime_Ocean']['vis_nir_ratio'][6] == 1: print("NIR-Visible ratio test running") radshape = rad1.shape rad1 = rad1.reshape(np.prod(radshape)) rad2 = rad2.reshape(np.prod(radshape)) sunglint = scene['sunglint'].reshape(np.prod(radshape)) vrat = rad2/rad1 confidence = np.ones(rad1.shape) tmp = threshold['Daytime_Ocean']['vis_nir_ratio'] thr_no_sunglint = np.array([tmp[0], tmp[1], tmp[2], tmp[3], tmp[4], tmp[5], 1, 1]) tmp = threshold['Sun_Glint']['snglnt'] thr_sunglint = np.array([tmp[0], tmp[1], tmp[2], tmp[3], tmp[4], tmp[5], 1]) # thr_no_sunglint = np.array(threshold['Daytime_Ocean']['vis_nir_ratio']) # thr_sunglint = np.array(threshold['Sun_Glint']['snglnt']) # thr_sunglint = np.append(thr_sunglint, 1) # temp value to avoid linter bitching at me # eventually we would have the test run in two blocks as: # confidence[sunglint == 1] = conf.conf_test_dble(vrat[sunglint == 1], sg_threshold['snglnt']) # confidence[sunglint == 0] = conf.conf_test_dble(vrat[sunglint == 0], threshold['vis_nir_ratio']) # sunglint needs to be defined somewhere # thr = np.full((rad.shape[0], 4), thr[:4]).T # thresh = np.full((rad1.shape[0], thr_no_sunglint.shape[0]), thr_no_sunglint) # thresh[sunglint == 1, :6] = thr_sunglint confidence[sunglint == 0] = conf.conf_test_dble(vrat, thr_no_sunglint)[sunglint == 0] confidence[sunglint == 1] = conf.conf_test_dble(vrat, thr_sunglint)[sunglint == 1] # confidence = conf.conf_test_dble(vrat, thresh.T) return np.minimum(cmin, confidence.reshape(radshape)), confidence.reshape(radshape) def nir_refl_test(rad, threshold, sunglint_thresholds, viirs_data, cmin): print("NIR reflectance test running") sza = viirs_data.solar_zenith.values refang = viirs_data.sunglint_angle.values vza = viirs_data.sensor_zenith.values dtr = np.pi/180 # Keep in mind that band_n uses MODIS band numbers (i.e. 2=0.86um and 7=2.1um) # For VIIRS would be 2=M07 (0.865um) and 7=M11 (2.25um) band_n = 2 vzcpow = 0.75 # THIS NEEDS TO BE READ FROM THE THRESHOLDS FILE radshape = rad.shape rad = rad.reshape(np.prod(radshape)) confidence = np.ones(rad.shape) sza = sza.reshape(rad.shape) vza = vza.reshape(rad.shape) refang = refang.reshape(rad.shape) sunglint_flag = utils.sunglint_scene(refang, sunglint_thresholds).reshape(rad.shape) # ref2 [5] # b2coeffs [4] # b2mid [1] # b2bias_adj [1] # b2lo [1] # vzcpow [3] (in different place) cosvza = np.cos(vza*dtr) coeffs = threshold['b2coeffs'] hicut0 = np.array(coeffs[0] + coeffs[1]*sza + coeffs[2]*np.power(sza, 2) + coeffs[3]*np.power(sza, 3)) hicut0 = (hicut0 * 0.01) + threshold['b2adj'] hicut0 = hicut0 * threshold['b2bias_adj'] midpt0 = hicut0 + (threshold['b2mid'] * threshold['b2bias_adj']) locut0 = midpt0 + (threshold['b2lo'] * threshold['b2bias_adj']) thr = np.array([locut0, midpt0, hicut0, threshold['ref2'][3]*np.ones(rad.shape)]) corr_thr = np.zeros((4, rad.shape[0])) corr_thr[:3, sunglint_flag == 0] = thr[:3, sunglint_flag == 0] * (1./np.power(cosvza[sunglint_flag == 0], vzcpow)) corr_thr[3, sunglint_flag == 0] = thr[3, sunglint_flag == 0] # corr_thr[:3, :] = thr[:3, :] * (1./np.power(cosvza[:], vzcpow)) # corr_thr[3, :] = thr[3, :] for flag in range(1, 4): if len(refang[sunglint_flag == flag]) > 0: sunglint_thr = utils.get_sunglint_thresholds(refang, sunglint_thresholds, band_n, flag, thr) corr_thr[:3, sunglint_flag == flag] = sunglint_thr[:3, sunglint_flag == flag] * (1./np.power(cosvza[sunglint_flag == flag], vzcpow)) corr_thr[3, sunglint_flag == flag] = sunglint_thr[3, sunglint_flag == flag] confidence = conf.conf_test(rad, corr_thr) return np.minimum(cmin, confidence.reshape(radshape)), confidence.reshape(radshape) def nir_high_cloud_test(): pass def test_11um_var(rad, threshold, var_threshold): print("11um variability test running") thr = np.array(threshold['11um_var']) radshape = rad.shape var = np.zeros((radshape[0], radshape[1], 9)) # chk_spatial2() need to figure out what this is # np = rg_var.num_small_diffs * 1.0 test = sliding_window_view(np.pad(rad, [1, 1], mode='constant'), (3, 3)) - np.expand_dims(rad, (2, 3)) var[np.abs(test).reshape(radshape[0], radshape[1], 9) < var_threshold['dovar11']] = 1 var = var.sum(axis=2).reshape(np.prod(radshape)) rad = rad.reshape(np.prod(radshape)) confidence = np.zeros(rad.shape) confidence[var == 9] = conf.conf_test(rad[var == 9], thr) return confidence.reshape(radshape) def test_11_4diff(rad1, rad2, threshold, viirs_data, sg_thresh): print("11um - 4um difference test running") radshape = rad1.shape raddiff = (rad1 - rad2).reshape(np.prod(radshape)) day = np.zeros(radshape) day[viirs_data.solar_zenith <= 85] = 1 day = day.reshape(raddiff.shape) sunglint = np.zeros(rad1.shape) sunglint[viirs_data.sunglint_angle <= sg_thresh] = 1 sunglint = sunglint.reshape(raddiff.shape) thr = np.array(threshold['test11_4lo']) confidence = np.zeros(raddiff.shape) # confidence[(day == 1) & (sunglint == 0)] = utils.conf_test(raddiff[(day == 1) & (sunglint == 0)], thr) confidence[(day == 1) & (sunglint == 0)] = conf.conf_test(raddiff[(day == 1) & (sunglint == 0)], thr) return confidence.reshape(radshape) def vir_refl_test(rad, threshold, viirs_data): print('Visible reflectance test running') thr = threshold['vis_refl_test'] radshape = rad.shape() rad = rad.reshape(np.prod(radshape)) confidence = np.zeros(radshape) vzcpow = 0.75 # THIS NEEDS TO BE READ FROM THE THRESHOLDS FILE vza = viirs_data.sensor_zenith.values dtr = np.pi/180 cosvza = np.cos(vza*dtr) coeffs = utils.get_b1_thresholds() coeffs[:, :3] = coeffs[:, :3] * threshold['b1_bias_adj'] # this quantity is the return of get_b1_thresholds() in the C code # it's defined here to keep a consistent logic with the original source, for now irtn = 0 if irtn != 0: coeffs = thr coeffs[:, :3] = coeffs[:, :3] * 1/np.power(cosvza, vzcpow) confidence = conf.conf_test(rad, coeffs) return confidence.reshape(radshape) class CloudMaskTests(object): def __init__(self, scene, radiance, coefficients): self.scene = scene self.coefficients = coefficients def select_coefficients(self): pass def test_G1(self): pass def test_G2(self): pass def test_G3(self): pass def test_G4(self): pass def overall_confidence(self): pass def test(): rad = np.random.randint(50, size=[4, 8]) # coeffs = [5, 42, 20, 28, 15, 35, 1] # coeffs = [20, 28, 5, 42, 15, 35, 1] coeffs = [35, 15, 20, 1, 1] # confidence = conf_test_dble(rad, coeffs) confidence = test_11um(rad, coeffs) print(rad) print('\n') print(confidence) if __name__ == "__main__": test()