@read_ecm_v1_5_lut @use_ecm_v1_5_lut @read_standard_hdf @ecm_v1_5_routines ;---------------------------------------------------------- pro run_ecm_v1_5_on_level2 , no_uni=no_uni, no_sub = no_sub, no_obv = no_obv Ecm_Lut_Name = 'n20_viirs_lut.nc' ;--- daytime arctic ;Level2_Name = 'level2_test/clavrx_snpp_viirs.A2015189.2242.001.2017314164212.uwssec_B00019147.level2.hdf' Level2_Name = 'level2_test/clavrx_snpp_viirs.A2019076.0224.001.2019076075638.uwssec_B00038260.level2.hdf' ;--- daytime gulf of mexico ;Level2_Name = 'clavrx_j1_viirs.A2019100.1742.001.2019100213235.uwssec_B00007223.level2.hdf' ;--- nightime central usa ;Level2_Name = 'level2_test/clavrx_snpp_viirs.A2019089.0900.001.2019089135957.uwssec_B00038448.level2.hdf' ;--- daytime scandanavia ;Level2_Name = 'clavrx_snpp_viirs.A2019089.1024.001.2019089154700.uwssec_B00038449.level2.hdf' ;--- nightime california and ocean ;Level2_Name = 'level2_test/clavrx_snpp_viirs.A2019072.0924.001.2019072140850.uwssec_B00038207.level2.hdf' Prior_Name = 'nb_cloud_mask_modis_prior.nc' Missing = -999.0 ;---- read ecm lut Ecm_Lut = read_ecm_v1_5_lut(Ecm_Lut_name) ;---- read level2 l2 = read_level2(Level2_Name) ;--- read prior Prior = read_prior(Prior_Name) ;--------------------------------------------------------------------------------------------- ; make binary masks of snow, land, coast ;--------------------------------------------------------------------------------------------- Land_Mask = make_array(l2.Num_Elem_Idx,l2.Num_Line_Idx,/INT,VALUE=0) idx = where(l2.Land_Class eq 1,cc) if (cc gt 0) then Land_Mask[idx] = 1 Snow_Mask = make_array(l2.Num_Elem_Idx,l2.Num_Line_Idx,/INT,VALUE=0) idx = where(l2.Snow_Class ge 2,cc) if (cc gt 0) then Snow_Mask[idx] = 1 ;--------------------------------------------------------------------------------------------- ; Make some metrics need for obvious cloud mask and rut and tut ;--------------------------------------------------------------------------------------------- t11_sub_flag = 0 ref065_sub_flag = 0 if (max(l2.bt11_max_sub) ne missing and max(l2.bt11_min_sub) ne missing) then t11_sub_flag = 1 if (max(l2.ref065_max_sub) ne missing and max(l2.ref065_min_sub) ne missing) then ref065_sub_flag = 1 if (keyword_set(no_sub)) then begin t11_sub_flag = 0 ref065_sub_flag = 0 endif ;--- make ref 065 metric ref065_metric = l2.ref065 - l2.ref065_clr idx = where(l2.ref065_max_sub gt l2.ref065,cc) if (cc gt 0) then ref065_metric[idx] = l2.ref065_max_sub[idx] - l2.ref065_clr[idx] ;--- make rvct metric rvct_metric = l2.ref065 - l2.ref065_min_3x3 if (ref065_sub_flag eq 1) then begin idx = where(l2.ref065_min_sub ne missing and l2.ref065_min_sub lt l2.ref065_min_3x3,cc) if (cc gt 0) then rvct_metric[idx] = l2.ref065[idx] - l2.ref065_min_sub[idx] endif ;--- make ref_ratio_metric ref_ratio_metric = l2.ref086 / l2.ref065 idx = where(l2.ref065 eq missing or l2.ref086 eq missing,cc) if (cc gt 0) then ref_ratio_metric[idx] = missing ;--- make ndsi_metric ndsi_metric = (l2.ref065 - l2.ref160) / (l2.ref065 + l2.ref160) idx = where(l2.ref065 eq missing or l2.ref160 eq missing,cc) if (cc gt 0) then ndsi_metric[idx] = missing ;--- make tut metric (use sub if possible) tut_metric = l2.bt11_stddev_3x3 if (t11_sub_flag eq 1) then begin tut_metric_sub = (l2.bt11_max_sub - l2.bt11_min_sub)/sqrt(2.0) idx = where(tut_metric_sub gt tut_metric,cc) if (cc gt 0) then tut_metric[idx] = tut_metric_sub[idx] endif ;--- make rut metric (use sub if possible) rut_metric = l2.ref065_stddev_3x3 if (ref065_sub_flag eq 1) then begin rut_metric_sub = (l2.ref065_max_sub - l2.ref065_min_sub)/sqrt(2.0) idx = where(rut_metric_sub gt rut_metric,cc) if (cc gt 0) then rut_metric[idx] = rut_metric_sub[idx] endif ;--- make tmax metric (use sub if possible) t11_rel_3x3_metric = l2.bt11_max_3x3 - l2.bt11 t11_rel_sub_metric = make_array(l2.Num_Elem_Idx,l2.Num_Line_Idx,/FLOAT,VALUE=MISSING) if (t11_sub_flag eq 1) then begin idx = where(l2.bt11_max_sub ne missing and l2.bt11_max_sub gt l2.bt11_max_3x3,cc) if (cc gt 0) then t11_rel_sub_metric[idx] = l2.bt11_max_sub[idx] - l2.bt11_min_sub[idx] endif ;--- make t11 metric t11_metric = l2.bt11 if (t11_sub_flag eq 1) then begin idx = where(l2.bt11_min_sub lt l2.bt11 and l2.bt11_min_sub ne missing,cc) if (cc gt 0) then t11_metric[idx] = l2.bt11_min_sub[idx] endif ;--------------------------------------------------------------------------------------------- ; Step 1: obvious cloud mask ; it will detect obvious cloud and set Post_Prob = 1 for those pixels. ;--------------------------------------------------------------------------------------------- Post_Prob_Obvious = obvious_cloud_detection(l2.Sfc_Idx,l2.coast_mask, l2.solzen, l2.glintzen, $ t11_metric, ref065_metric, $ t11_rel_3x3_metric, t11_rel_sub_metric, t11_sub_flag) ;--------------------------------------------------------------------------------------------- ; Step 2: apply naive bayesian mask ;--------------------------------------------------------------------------------------------- ;-- make arrays needed Classifier_Data = make_array(Ecm_Lut.N_Class,/FLOAT,VALUE=Missing) Post_Prob = make_array(l2.Num_Elem_Idx,l2.Num_Line_Idx,/FLOAT,VALUE=Missing) Prior_Prob = make_array(l2.Num_Elem_Idx,l2.Num_Line_Idx,/FLOAT,VALUE=Missing) Post_Prob_by_Class = make_array(l2.Num_Elem_Idx,l2.Num_Line_Idx,Ecm_Lut.N_Class,/FLOAT,VALUE=Missing) Post_Prob_by_Class_Temp = make_array(Ecm_Lut.N_Class,/FLOAT,VALUE=Missing) ;--- initialize to results of obvious mask Post_Prob = Post_Prob_Obvious if (keyword_set(no_obv)) then Post_Prob[*,*] = Missing for Line_Idx = 0, l2.Num_Line_Idx - 1 do begin for Elem_Idx = 0, l2.Num_Elem_Idx - 1 do begin ;--- check if filled in already by simple mask if (Post_Prob[Elem_Idx,Line_Idx] ne Missing) then goto, cycle ;--- compute prior from our map Prior_Prob[Elem_Idx,Line_Idx] = COMPUTE_PRIOR(Prior, l2.Lon[Elem_Idx,Line_Idx], l2.Lat[Elem_Idx,Line_Idx], l2.Month) ;--- fill classifier data if (Ecm_Lut.T_11_Class_Idx ge 0) then begin Classifier_Data[Ecm_Lut.T_11_Class_Idx] = t11_metric[Elem_Idx,Line_Idx] endif if (Ecm_Lut.Ref_063_Day_Class_Idx ge 0) then begin Classifier_Data[Ecm_Lut.Ref_063_Day_Class_Idx] = ref065_metric[Elem_Idx,Line_Idx] endif if (Ecm_Lut.Ref_063_Min_3x3_Day_Class_Idx ge 0) then begin Classifier_Data[Ecm_Lut.Ref_063_Min_3x3_Day_Class_Idx] = rvct_metric[Elem_Idx,Line_Idx] endif if (Ecm_Lut.Ref_138_Day_Class_Idx ge 0) then begin Classifier_Data[Ecm_Lut.Ref_138_Day_Class_Idx] = l2.ref138[Elem_Idx,Line_Idx] endif if (Ecm_Lut.Ref_Ratio_Day_Class_Idx ge 0) then begin Classifier_Data[Ecm_Lut.Ref_Ratio_Day_Class_Idx] = ref_ratio_metric[Elem_Idx,Line_Idx] endif if (Ecm_Lut.Ndsi_Day_Class_Idx ge 0) then begin Classifier_Data[Ecm_Lut.Ndsi_Day_Class_Idx] = ndsi_metric[Elem_Idx,Line_Idx] endif if (Ecm_Lut.T_Max_Class_Idx ge 0) then begin Classifier_Data[Ecm_Lut.T_Max_Class_Idx] = t11_rel_3x3_metric[Elem_Idx,Line_Idx] endif if (Ecm_Lut.FMFT_Class_Idx ge 0) then begin fmft = (l2.bt11[Elem_Idx,Line_Idx] - l2.bt12[Elem_Idx,Line_Idx]) * $ (l2.bt11[Elem_Idx,Line_Idx]-260.0) / (l2.bt11_clr[Elem_Idx,Line_Idx] - 260.0) if (l2.bt11_clr[ELem_Idx,Line_Idx] le 265.0) then fmft = 0.0 Classifier_Data[Ecm_Lut.FMFT_Class_Idx] = (l2.bt11[Elem_Idx,Line_Idx] - l2.bt12[Elem_Idx,Line_Idx]) - fmft endif if (Ecm_Lut.Emiss_Tropo_Class_Idx ge 0) then begin Classifier_Data[Ecm_Lut.Emiss_Tropo_Class_Idx] = l2.etrop[Elem_Idx,Line_Idx] endif if (Ecm_Lut.Btd_11_85_Class_Idx ge 0) then begin Classifier_Data[Ecm_Lut.Btd_11_85_Class_Idx] = l2.bt11[Elem_Idx,Line_Idx] - l2.bt85[Elem_Idx,Line_Idx] endif if (Ecm_Lut.Btd_375_11_Day_Class_Idx ge 0) then begin Classifier_Data[Ecm_Lut.Btd_375_11_Day_Class_Idx] = $ (l2.bt375[Elem_Idx,Line_Idx] - l2.bt11[Elem_Idx,Line_Idx]) - $ (l2.bt375_clr[Elem_Idx,Line_Idx] - l2.bt11_clr[Elem_Idx,Line_Idx]) endif if (Ecm_Lut.Btd_375_11_Night_Class_Idx ge 0) then begin Classifier_Data[Ecm_Lut.Btd_375_11_Night_Class_Idx] = $ (l2.bt375[Elem_Idx,Line_Idx] - l2.bt11[Elem_Idx,Line_Idx]) - $ (l2.bt375_clr[Elem_Idx,Line_Idx] - l2.bt11_clr[Elem_Idx,Line_Idx]) endif ;--- compute probabilities use_ecm_v1_5_lut, Ecm_Lut, $ l2.Sfc_Idx[Elem_Idx,Line_Idx], $ ;note this 1 based Classifier_Data, $ l2.Solzen[Elem_Idx,Line_Idx], $ l2.Glintzen[Elem_Idx,Line_Idx], $ l2.Tsfc[Elem_Idx,Line_Idx], $ l2.Zsfc[Elem_Idx,Line_Idx], $ l2.Coast_Mask[Elem_Idx,Line_Idx], $ Prior_Prob[Elem_Idx,Line_Idx], $ Post_Prob_By_Class_Temp, $ Post_Prob_Temp Post_Prob_By_Class[Elem_Idx,Line_Idx,*] = Post_Prob_By_Class_Temp Post_Prob[Elem_Idx,Line_Idx] = Post_Prob_Temp cycle: endfor print, 'done with j = ', Line_Idx endfor ;--------------------------------------------------------------------------------------------- ; Step 3: apply tut and rut ;--------------------------------------------------------------------------------------------- Tut = make_array(l2.Num_Elem_Idx,l2.Num_Line_Idx,/INT,VALUE=0) Rut = make_array(l2.Num_Elem_Idx,l2.Num_Line_Idx,/INT,VALUE=0) if (~keyword_set(no_uni)) then begin for Line_Idx = 0, l2.Num_Line_Idx - 1 do begin for Elem_Idx = 0, l2.Num_Elem_Idx - 1 do begin Tut[Elem_Idx,Line_Idx] = TUT_ROUTINE(Land_Mask[Elem_Idx,Line_Idx], $ l2.Coast_Mask[Elem_Idx,Line_Idx], $ Tut_Metric[Elem_Idx,Line_Idx], $ l2.Zsfc_Stddev_3x3[Elem_Idx,Line_Idx]) Rut[Elem_Idx,Line_Idx] = RUT_ROUTINE(Land_Mask[Elem_Idx,Line_Idx], $ Snow_Mask[Elem_Idx,Line_Idx], $ l2.Coast_Mask[Elem_Idx,Line_Idx], $ l2.Ref065_Clr[Elem_Idx,Line_Idx], $ Rut_Metric[Elem_Idx,Line_Idx], $ l2.SolZen[Elem_Idx,Line_Idx]) endfor endfor endif ;------------------------------------------------------------------- ; Make Cloud Mask ;------------------------------------------------------------------- Cloud_Mask = make_array(l2.Num_Elem_Idx,l2.Num_Line_Idx,/INT,VALUE=Missing) Cloud_Mask_Binary = make_array(l2.Num_Elem_Idx,l2.Num_Line_Idx,/INT,VALUE=Missing) for Line_Idx = 0, l2.Num_Line_Idx - 1 do begin for Elem_Idx = 0, l2.Num_Elem_Idx - 1 do begin Make_Cloud_Mask, l2.Sfc_Idx[Elem_Idx,Line_Idx], $ Post_Prob[Elem_Idx,Line_Idx], $ Rut[Elem_Idx,Line_Idx], $ Tut[Elem_Idx,Line_Idx], $ Cloud_Mask_Temp, $ Cloud_Mask_Binary_Temp Cloud_Mask[Elem_Idx,Line_Idx] = Cloud_Mask_Temp Cloud_Mask_Binary[Elem_Idx,Line_Idx] = Cloud_Mask_Binary_Temp endfor endfor save, file = level2_name + ".output.sav", Cloud_Mask, Cloud_Mask_Binary, Post_Prob, Post_Prob_by_Class, Post_Prob_Obvious, Rut, Tut, $ rut_metric, tut_metric, t11_metric, t11_rel_3x3_metric, t11_rel_sub_metric, ref065_metric, l2, snow_mask, land_mask, /compress stop end