"$Id:$" ;--------------------------------------------------------------- ; preamble ;--------------------------------------------------------------- @read_group_attributes @read_1d_group_lut @read_2d_group_lut @read_3d_group_lut @read_prior @compute_prior_from_map @read_level2 @get_prob_mask_phase @grab_prior_from_lut @leapyr @compute_month_day @make_true_phase @make_bce_metric @make_acc_metrics @make_cld_frac_metric @modify_on_flag @compute_clear_cod @ecm_sfc_type @plot_clear_cod_and_hit_rate @plot_latitude_performance @plot_prob_histo @plot_view_angle_performance @read_ecm_lut @apply_ecm_lut ;--------------------------------------------------------------- ; test ecm combined table on a training data file ;--------------------------------------------------------------- pro test_on_training, sensor=sensor, optimize=optimize, global = global, Time_Diff_Max = Time_Diff_Max, $ Cod_Val_Min_Thresh = Cod_Val_Min_Thresh, nomod=nomod, add_thresh=add_thresh,no_plot = no_plot, $ Opt_Metric = Opt_Metric, Report_Metric = Report_Metric, isccp = isccp, force_dnb = force_dnb, $ use_phase = use_phase, phase_opt = phase_opt, Metric_Thresh = Metric_Thresh, Nthin_Frac = Nthin_Frac if (~keyword_set(Cod_Val_Min_Thresh)) then Cod_Val_Min_Thresh = 0.01 if (~keyword_set(Time_Diff_Max)) then Time_Diff_Max = 10.0 ; minutes if (~keyword_set(sensor)) then sensor = 'abhi' if (~keyword_set(Opt_Metric)) then Opt_Metric = 'BWACC' if (~keyword_set(Report_Metric)) then Report_Metric = 'BACC' if (~keyword_set(Metric_Thresh)) then Metric_Thresh = 1.0e-04 if (~keyword_set(Nthin_Frac)) then Nthin_Frac = 1.0 ;-- control training data modification in make_true_phase routine mod_flag = 1 if (keyword_set(nomod)) then mod_flag = 0 if (~keyword_set(use_phase)) then use_phase = 0 if (~keyword_set(phase_opt)) then phase_opt = 0 ;------------------------------------------------------------------------------ ; user input ;------------------------------------------------------------------------------ cc_min = 10 Prior_Source = 1 NClass_Optimal_Thresh = 5 cfrac_min = 0.1 cfrac_max = 0.9 R_Prod_Min = double(0.0) R_Prod_Max = double(1000.0) R_Log_Max = double(40.0) R_Log_Min = double(-40.0) Nsfc = 7 Training_Name = 'training_data/ecm_training_'+sensor+'_test.sav' nc_filename = 'combined_tables/test.nc' ;nc_filename = 'viirs_smoke/optimal_lut.nc' ;Prior_Name = 'prior/nb_cloud_mask_modis_prior.nc' Prior_Name = 'prior/nb_cloud_mask_calipso_prior.nc' Missing = -999.0 Prob_Min_Thresh = 0.001 Sfc_Type_Name = ['global','ocean','shallow water','land','snow','arctic','antarctic','desert'] Output_Name = 'output_on_training.sav' ;------------------------------------------------------------------------------ ; set channel on/off and channel mapping ; Chan_On = channel on / off binary flags ; Chan_Wvl = channel nominal wavelength in nm ; This is set up using CLAVR-x channel numbers ;------------------------------------------------------------------------------ Nchan = 45 Chan_On = make_array(Nchan,/INT,VALUE=1) Chan_On[31-1] = 1 Chan_Wvl = make_array(Nchan,/INT,VALUE=0) Chan_On[1-1] = 1 & Chan_Wvl[1-1] = 650 Chan_On[2-1] = 1 & Chan_Wvl[2-1] = 860 Chan_On[26-1] = 1 & Chan_Wvl[26-1] = 1380 Chan_On[6-1] = 1 & Chan_Wvl[6-1] = 1600 Chan_On[20-1] = 1 & Chan_Wvl[20-1] = 3750 Chan_On[27-1] = 1 & Chan_Wvl[27-1] = 6700 Chan_On[28-1] = 1 & Chan_Wvl[28-1] = 7300 Chan_On[29-1] = 1 & Chan_Wvl[29-1] = 8500 Chan_On[30-1] = 1 & Chan_Wvl[30-1] = 9700 Chan_On[31-1] = 1 & Chan_Wvl[31-1] = 11000 Chan_On[32-1] = 1 & Chan_Wvl[32-1] = 12000 Chan_On[33-1] = 1 & Chan_Wvl[33-1] = 13300 Chan_On[37-1] = 1 & Chan_Wvl[37-1] = 6200 Chan_On[38-1] = 1 & Chan_Wvl[38-1] = 10400 Chan_On[44-1] = 1 & Chan_Wvl[44-1] = 700 ;DNB ;------------------------------------------------------------------- ; read lookup table ;------------------------------------------------------------------- read_ecm_lut, nc_filename, nclass, attrs, luts, channels_used, Class_Use_Flag, Class_Rank, $ Class_Xname, Class_Yname, Class_Zname, Classifier_Names, $ Conf_Clear_Prob_Clear_Thresh, Prob_Clear_Prob_Cloud_Thresh, Prob_Cloud_Conf_Cloud_Thresh, $ cod_clear_ice_thresh, cod_clear_water_thresh, 0 ;------------------------------------------------------------------- ; turn off tables based on chan_on from the ecm table ;------------------------------------------------------------------- for Class_Idx = 0, Nclass -1 do begin if (Class_Use_Flag[Class_Idx] eq 0) then goto, skip_class_0 x = *(Channels_Used[Class_Idx]) for Chan_Idx = 0, Nchan - 1 do begin if (Chan_On[Chan_Idx] eq 0) then begin idx = where(Chan_Wvl[Chan_Idx] eq X.Wvl and Chan_Wvl[Chan_Idx] ne 0,cc) if (cc gt 0) then begin Class_Use_Flag[Class_Idx] = 0 endif endif endfor skip_class_0: ; print, Classifier_Names[Class_Idx], Class_Use_Flag[Class_Idx] endfor Number_Valid_Classes = total(Class_Use_Flag) if (Number_Valid_Classes eq 0) then begin print, 'no active tests' stop endif ;------------------------------------------------------------------- ; read prior table ;------------------------------------------------------------------- read_prior, Prior_Name, Prior ;--------------------------------------------------------------------------------------------- ; read training data ;--------------------------------------------------------------------------------------------- restore, training_name l2 = training_data training_data = !null l2_nx = n_elements(l2.latitude) ;---- adjust nasa dnb to look like noaa dnb [REMOVE WHEN ADJUSTMENT in TRAINING DATA] Dnb_Coef = [-0.373685,0.977945,-0.000261637] idx = where(l2.refl_lunar_dnb_nom ne missing,cc) if (cc gt 0) then begin l2[idx].refl_lunar_dnb_nom = Dnb_Coef[0] + $ Dnb_Coef[1]*l2[idx].refl_lunar_dnb_nom + $ Dnb_Coef[2]*l2[idx].refl_lunar_dnb_nom^2 endif ;------------------------------------------------------------------------------------------- ; ;------------------------------------------------------------------------------------------- if (keyword_set(isccp)) then begin make_ecm_surface_type, L2, sfc_type_mask, /isccp L2.bayes_mask_sfc_type = sfc_type_mask endif ;-------------------------------------------------------------------------------------------- ; randomly pick a certain number for testing ;-------------------------------------------------------------------------------------------- ;N_Thin = Nthin_Frac * L2_Nx ;x = round((l2_nx-1) * randomu(seed,l2_nx)) ;idx_thin = x[0:N_thin-1] ;l2 = l2[idx_thin] ;print, 'number of training points = ', l2_nx ;--- remove bad data idx = where(l2.latitude ne Missing,cc) l2 = l2[idx] l2_nx = cc print, 'number of testing points = ', l2_nx ;-------------------------------------------------------------------------------------- ; filter the calipso data to approximate a better truth ; modify the phase_lidar based on cod thresholds and make Phase_True ;-------------------------------------------------------------------------------------- ;--- calipso filtering cod_cloud_thresh = 1.0 ;--- threshold for passive cod (not used) etrop_cloud_thresh = 0.9 ;--- threshold for passie etrop (not used) zen_min = 0.0 & zen_max = 180.0 solzen_min = 0.0 & solzen_max = 180.0 lunzen_min = 0.0 & lunzen_max = 180.0 solglintzen_min = 0.0 & solglintzen_max = 180.0 lunglintzen_min = 0.0 & lunglintzen_max = 180.0 solscatang_min = 0.0 & solscatang_max = 180.0 snow_class_min = 1 & snow_class_max = 3 solglint_mask_min = 0 & solglint_mask_max = 1 lunglint_mask_min = 0 & lunglint_mask_max = 1 city_mask_min = 0 city_mask_max = 1 moon_illum_frac_min = 0.0 moon_illum_frac_max = 1.0 coast_mask_min = 0 & coast_mask_max = 1 tpw_min = 0.0 & tpw_max = 20.0 tsfc_min = 180.0 & tsfc_max = 340.0 zsfc_min = -500.0 & zsfc_max = 50000.0 zsfc_std_min = 0.0 & zsfc_std_max = 1000.0 lat_min = -90.0 & lat_max = 90.0 phase_lidar = L2.CLOSEST_CALIPSO_CLOUD_PHASE_TOP_LAYER frac_lidar = L2.CLOSEST_CALIPSO_CLOUD_FRACTION_5KM tc_lidar = L2.CLOSEST_CALIPSO_TOP_MID_TEMPERATURE ;--- convert tc_lidar from C to K idx = where(tc_lidar ne missing, cc) if (cc gt 0) then tc_lidar[idx] = tc_lidar[idx] + 273.15 ;phase_lidar = L2.CLOSEST_CALIPSO_PHASE_TOP_LAYER ;frac_lidar = L2.CLOSEST_CALIPSO_CLOUD_FRACTION ;tc_lidar = L2.CLOSEST_CALIPSO_MID_TEMPERATURE Cod_Lidar = L2.CLOSEST_CALIPSO_COD codi_lidar = L2.CLOSEST_CALIPSO_COD_ICE print,'making true phase' make_true_phase, phase_lidar, frac_lidar, tc_lidar, Cod_Lidar, codi_lidar, $ cod_clear_water_thresh, cod_clear_ice_thresh, $ L2.bayes_mask_sfc_type, L2.cld_Temp_opaque, L2.cld_height_opaque, $ L2.solar_zenith_angle, L2.glint_mask, L2.coast_mask, $ L2.surface_elevation, L2.surface_Temperature_nwp, $ L2.refl_0_65um_nom_min_3x3, L2.refl_0_65um_nom, $ L2.temp_11_0um_nom_max_3x3, L2.temp_11_0um_nom, $ L2.cld_opd_mask_0_65um_nom, cod_cloud_thresh, $ L2.emiss_tropo_11_0um_nom, etrop_cloud_thresh, mod_flag, Phase_True, missing ;idx = where(Phase_True eq 0 and L2.solar_zenith_angle lt 80.0 and $ ; L2.refl_0_65um_nom_stddev_sub gt 5 and $ ; L2.refl_0_65um_nom_max_sub gt 10, cc) ;if (cc gt 0) then begin ; Phase_True[idx] = 1 ; Cod_Lidar[idx] = 3.0 ;endif ;print, 'Number of new water = ', cc ;------------------------------------------------------------------------------- ; fill in things that might be missing in training data ;------------------------------------------------------------------------------- idx = where(tag_names(L2) eq "CITY_MASK",cc) if (cc eq 0) then begin city_mask_temp = make_array(l2_nx,/FLOAT, VALUE=MISSING) endif else begin city_mask_temp = L2.city_mask endelse idx = where(tag_names(L2) eq "MOON_ILLUM_FRAC",cc) if (cc eq 0) then begin moon_illum_frac_temp = make_array(l2_nx,/FLOAT, VALUE=MISSING) endif else begin moon_illum_frac_temp = L2.moon_illum_frac endelse slant_tpw = L2.total_precipitable_water_nwp / cos(L2.sensor_zenith_angle * !DTOR) print,'applying filters' ; ; why do this? ; apply_filters, Phase_True, $ L2.CLOSEST_CALIPSO_TIME_DIFFERENCE, Time_Diff_Max, $ L2.sensor_zenith_angle, zen_min, zen_max, $ L2.solar_zenith_angle, solzen_min, solzen_max, $ L2.lunar_zenith_angle, lunzen_min, lunzen_max, $ L2.land_class, $ frac_lidar, cfrac_min, cfrac_max, $ L2.glint_zenith_angle, solglintzen_min, solglintzen_max, $ L2.lunar_glint_zenith_angle, lunglintzen_min, lunglintzen_max, $ L2.scattering_angle, solscatang_min, solscatang_max, $ L2.snow_class, snow_class_min, snow_class_max, $ L2.glint_mask, solglint_mask_min, solglint_mask_max, $ L2.glint_mask_lunar, lunglint_mask_min, lunglint_mask_max, $ city_mask_temp, city_mask_min, city_mask_max, $ moon_illum_frac_temp, moon_illum_frac_min, moon_illum_frac_max, $ L2.coast_mask, coast_mask_min, coast_mask_max, $ slant_tpw, tpw_min, tpw_max, $ L2.surface_Temperature_nwp, tsfc_min, tsfc_max, $ L2.surface_elevation, zsfc_min, zsfc_max, $ L2.surface_elevation_stddev_3x3, zsfc_std_min, zsfc_std_max, $ L2.latitude, lat_min, lat_max, missing ;--------------------------------------------------------------------------------------------- ; compute prior ;--------------------------------------------------------------------------------------------- print,'computing prior' Diurnal = 1 ;--- make month array Month = make_array(L2_Nx,/Int,Value=-1) for i = 0, l2_nx - 1 do begin if (l2[i].doy gt 0 and l2[i].year gt 0) then begin compute_month_day, l2[i].doy, leapyr(fix(l2[i].year)), xx,yy Month[i] = xx endif endfor ;---- call prior routine compute_prior_from_map_1d, l2.Longitude, l2.Latitude, l2.Bad_Pixel_Mask, l2_Nx, Missing, $ prior.lon_min, prior.dLon, prior.Nlon, $ prior.lat_min, prior.dLat, prior.Nlat, $ Prior, Month, Diurnal, $ Pixel_Prior_Prob_Cloud, Pixel_Prior_Prob_Ice_Cloud, Pixel_Prior_Prob_Water_Cloud ;--------------------------------------------------------------------------------------------- ; run data from ecm ; --------------------------------------------------------------------------------------------- apply_ecm_lut, Missing, R_Prod_Max, L2_Nx, Nclass, L2, luts, attrs, $ Class_Use_Flag, Class_Rank, Class_Xname, Class_Yname, Class_Zname, $ Clear_Cond_Ratio, Water_Cond_Ratio, Ice_Cond_Ratio, $ Obs_Prob, Post_Prob_Cloud, Post_Prob_Cloud_By_Class, $ Post_Prob_Clear, Post_Prob_Water, Post_Prob_Ice, $ Cloud_Phase, Cloud_Phase_Uncer, $ Prior_Source, Pixel_Prior_Prob_Cloud, Pixel_Prior_Prob_Ice_Cloud, Pixel_Prior_Prob_Water_Cloud, $ Prior_Prob_Clear, Prior_Prob_Ice, Prior_Prob_Water, $ city_mask_temp, moon_illum_frac_temp ;------------------------------------------------------------------------------------------------- ; make histogram of Post_Prob_Cloud ;------------------------------------------------------------------------------------------------- if (~keyword_set(no_plot)) then begin plot_prob_histo, Nclass, Missing, Post_Prob_Cloud, Post_Prob_Cloud_by_Class, Phase_True endif ;------------------------------------------------------------------------------------------------ ; determine conf/prob boundaries here for each surface type ;------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------- ; correlation of classifiers ;------------------------------------------------------------------------------------------------- Correlation_Matrix = make_array(Nclass,Nclass,/FLOAT,value=missing) npoint = L2_Nx ncorr = 10000 for i_class = 0, Nclass -1 do begin for j_class = 0, Nclass -1 do begin XX = Post_Prob_Cloud_by_Class[*,i_class] YY = Post_Prob_Cloud_by_Class[*,j_class] idx = where(XX ne missing and YY ne missing, npoint) if (npoint gt 0) then begin nstride = max([1,npoint / ncorr]) Correlation_Matrix[i_class, j_class] = correlate(XX[idx[0:npoint-1:nstride]],YY[idx[0:npoint-1:nstride]]) endif endfor endfor openw, lun, 'optimal_output.txt', /get_lun print, '================================================================' print, 'Correlation Matrix' print, '================================================================' printf, lun, '================================================================' printf, lun, 'Correlation Matrix' printf, lun, '================================================================' for j_class = 0, Nclass -1 do begin print, format="(a25, 25F6.2)", Classifier_Names[j_class], Correlation_Matrix[*,j_class] printf, lun, format="(a25, 25F6.2)", Classifier_Names[j_class], Correlation_Matrix[*,j_class] endfor ;------------------------------------------------------------------------------------------------- ; make the binary mask ;------------------------------------------------------------------------------------------------- Cloud_Mask_Binary = make_array(L2_Nx,/INT,VALUE=-1) Cloud_Mask_Binary_By_Class = make_array(L2_Nx,Nclass,/Int,Value=-1) idx = where(Cloud_Phase eq 0,cc) if (cc gt 0) then begin Cloud_Mask_Binary[idx] = 0 endif idx = where(Cloud_Phase gt 0,cc) if (cc gt 0) then begin Cloud_Mask_Binary[idx] = 1 endif ;--- binary mask by class temp_out = make_array(L2_Nx,/INT,VALUE=-1) for Class_Idx = 0, Nclass - 1 do begin temp_out[*] = -1 temp_in = reform(Post_Prob_Cloud_By_Class[*,Class_Idx]) idx = where (temp_in lt 0.5 and temp_in ne Missing,cc) if (cc gt 0) then temp_out[idx] = 0 idx = where (temp_in ge 0.5,cc) if (cc gt 0) then temp_out[idx] = 1 Cloud_Mask_Binary_By_Class[*,Class_Idx] = temp_out endfor ;--- make the 4-level mask Cloud_Mask = make_array(L2_Nx,/INT,VALUE=-1) for Isfc = 1, Nsfc do begin idx_sfc = where(l2.bayes_mask_sfc_type eq Isfc,cc) idx = where (Cloud_Mask_Binary[idx_sfc] eq 0 and Post_Prob_Cloud[idx_sfc] lt Conf_Clear_Prob_Clear_Thresh[Isfc-1],cc) if (cc gt 0) then Cloud_Mask[idx_sfc[idx]] = 0 idx = where (Cloud_Mask_Binary[idx_sfc] eq 0 and Post_Prob_Cloud[idx_sfc] ge Conf_Clear_Prob_Clear_Thresh[Isfc-1], cc) if (cc gt 0) then Cloud_Mask[idx_sfc[idx]] = 1 idx = where (Cloud_Mask_Binary[idx_sfc] eq 1 and Post_Prob_Cloud[idx_sfc] lt Prob_Cloud_Conf_Cloud_Thresh[Isfc-1], cc) if (cc gt 0) then Cloud_Mask[idx_sfc[idx]] = 2 idx = where (Cloud_Mask_Binary[idx_sfc] eq 1 and Post_Prob_Cloud[idx_sfc] ge Prob_Cloud_Conf_Cloud_Thresh[Isfc-1], cc) if (cc gt 0) then Cloud_Mask[idx_sfc[idx]] = 3 endfor ;-------------------------------------------------------------------------------- ; metrics ;-------------------------------------------------------------------------------- metric_mask = make_array(Nsfc+1,/FLOAT,Value=missing) Cld_Frac = make_array(Nsfc+1,/FLOAT,Value=missing) Cld_Frac_True = make_array(Nsfc+1,/FLOAT,Value=missing) for Sfc_Idx = 0, 7 do begin if (Sfc_Idx eq 0) then begin idx = where(Phase_True ge 0 and Cloud_Mask_Binary ge 0 and $ l2.bayes_mask_sfc_type ge 1,cc) cc_all = cc endif else begin idx = where(Phase_True ge 0 and Cloud_Mask_Binary ge 0 and $ l2.bayes_mask_sfc_type eq Sfc_Idx,cc) endelse if (cc gt 0) then begin idx_sub = where(Phase_True[idx] eq 0 or Cod_Lidar[idx] gt Cod_Val_Min_Thresh,cc) idx_sub_day = where(L2[idx_sub].solar_zenith_angle lt 90 and L2[idx_sub].solar_zenith_angle ge 0.0, cc_day) make_acc_metrics, Phase_true[idx[idx_sub]], Post_Prob_Cloud[idx[idx_sub]], $ Acc_Temp,Bacc_Temp, Wacc_Temp,Bwacc_Temp, Bce_Temp, $ Tpr_Temp,Tnr_Temp,Fpr_Temp,Fnr_Temp make_cld_frac_metric, Phase_True[idx[idx_sub]], Post_Prob_Cloud[idx[idx_sub]], $ Cld_Frac_True_Temp,Cld_Frac_Ret_Temp if (Report_Metric eq 'ACC') then Metric_Clear = Acc_Temp if (Report_Metric eq 'BACC') then Metric_Clear = Bacc_Temp if (Report_Metric eq 'WACC') then Metric_Clear = Wacc_Temp if (Report_Metric eq 'BWACC') then Metric_Clear = Bwacc_Temp if (Report_Metric eq 'BCE') then Metric_Clear = Bce_Temp metric_mask[Sfc_Idx] = metric_clear Cld_Frac_True[Sfc_Idx] = Cld_Frac_True_Temp Cld_Frac[Sfc_Idx] = Cld_Frac_Ret_Temp print, Sfc_Idx, Report_Metric, Metric_Clear ; ;--- phase metric ; idx_sub = where(Phase_True[idx] gt 0 and Cod_Lidar[idx] gt Cod_Val_Min_Thresh,cc) ; ii = where(Phase_True ; make_phase_metric, Phase_True[idx[idx_sub]], Post_Prob_Ice[idx[idx_sub]], Post_Prob_Water[idx[idx_sub]], $ ; Phase_Metric, phase_opt ; print, 'Sfc_Idx Representation (all and relative day) = ', float(cc) / float(cc_all), float(cc_day) / float(cc) endif endfor ;----------------------------------------------------- ; now by class ;----------------------------------------------------- metric_mask_by_class = make_array(Nsfc+1,Nclass,/FLOAT,Value=missing) for Class_Idx = 0, Nclass -1 do begin temp = Cloud_Mask_Binary_By_Class[*,Class_Idx] for Sfc_Idx = 0, 7 do begin if (Sfc_Idx eq 0) then begin idx = where(Phase_True ge 0 and temp ge 0 and $ l2.bayes_mask_sfc_type ge 1,cc) endif else begin idx = where(Phase_True ge 0 and temp ge 0 and $ l2.bayes_mask_sfc_type eq Sfc_Idx,cc) endelse if (cc gt 0) then begin idx_sub = where(Phase_True[idx] eq 0 or Cod_Lidar[idx] gt Cod_Val_Min_Thresh,cc) ;-- note, this is not sending the phase, but a binary mask ;-- so mask metrics are fine, not phase prob_temp = Post_Prob_Cloud_By_Class[*,Class_Idx] make_acc_metrics, Phase_true[idx[idx_sub]], Prob_Temp[idx[idx_sub]], Acc_Temp,Bacc_Temp,Wacc_Temp,Bwacc_Temp, Bce_Temp, $ Tpr_Temp,Tnr_Temp,Fpr_Temp,Fnr_Temp make_bce_metric, phase_true[idx[idx_sub]], prob_Temp[idx[idx_sub]], bce_Temp print, Class_Idx, Sfc_Idx, Acc_Temp,Bacc_Temp,Wacc_Temp,Bwacc_Temp, Bce_Temp, Tpr_Temp,Tnr_Temp,Fpr_Temp,Fnr_Temp if (Report_Metric eq 'ACC') then Metric_Clear = Acc_Temp if (Report_Metric eq 'BACC') then Metric_Clear = Bacc_Temp if (Report_Metric eq 'WACC') then Metric_Clear = Wacc_Temp if (Report_Metric eq 'BWACC') then Metric_Clear = Bwacc_Temp if (Report_Metric eq 'BCE') then Metric_Clear = Bce_Temp metric_mask_by_class[Sfc_Idx,Class_Idx] = metric_clear endif endfor endfor ;--- output pc print, '===============================================================================' print, 'Probability of Correct for Cloud by Class and by Sfc, No filter' print, '===============================================================================' printf, lun, '===============================================================================' printf, lun, 'Probability of Correct for Cloud by Class and by Sfc, No filter' printf, lun, '===============================================================================' print, format="(a40,a7,I2,6I5)", ' ', 'global', indgen(7)+1 print, format="(a40,20F5.2)", 'All', Metric_Mask printf, lun, format="(a40,a7,I2,6I5)", ' ', 'global', indgen(7)+1 printf, lun, format="(a40,20F5.2)", 'All', Metric_Mask for Class_Idx = 0,NClass-1 do begin print, format="(a40,20F5.2)", Classifier_Names[Class_Idx], Metric_Mask_By_Class[*,Class_Idx] printf, lun, format="(a40,20F5.2)", Classifier_Names[Class_Idx], Metric_Mask_By_Class[*,Class_Idx] endfor ;-------------------------------------------------------------------------------------------- ;-- By Surface ;-------------------------------------------------------------------------------------------- Metric_Clear_By_Sfc = make_array(Nsfc+1,/FLOAT, VALUE=missing) for Sfc_Idx = 0, Nsfc do begin if (Sfc_Idx eq 0) then begin idx = where (Phase_True ge 0,cc) endif else begin idx = where (Sfc_Idx eq fix(l2.bayes_mask_sfc_type) and Phase_True ge 0,cc) endelse make_acc_metrics, phase_true[idx], Post_Prob_Cloud[idx], $ Acc_Temp,Bacc_Temp,Wacc_Temp,Bwacc_Temp, Bce_Temp, $ Tpr_Temp,Tnr_Temp,Fpr_Temp,Fnr_Temp if (Report_Metric eq 'ACC') then Metric_Clear = Acc_Temp if (Report_Metric eq 'BACC') then Metric_Clear = Bacc_Temp if (Report_Metric eq 'WACC') then Metric_Clear = Wacc_Temp if (Report_Metric eq 'BWACC') then Metric_Clear = Bwacc_Temp if (Report_Metric eq 'BCE') then Metric_Clear = Bce_Temp Metric_Clear_By_Sfc[Sfc_Idx] = Metric_Clear endfor ;====================================================================================== ; plot performance as a function of view angle ;====================================================================================== if (~keyword_set(no_plot)) then begin plot_view_angle_performance, L2.Sensor_Zenith_Angle, Cloud_Phase, Phase_True, Post_Prob_Cloud, Report_Metric endif ;====================================================================================== ; plot performance as a function of latitude ;====================================================================================== if (~keyword_set(no_plot)) then begin plot_latitude_performance, L2.Latitude, Post_Prob_Cloud, Phase_True, Report_Metric endif ;-------------------------------------------------------------- ; Set Prob Thresh for 4-Level Mask ;-------------------------------------------------------------- if (keyword_set(add_thresh)) then begin Sfc_Idx = l2.bayes_mask_sfc_type ;save, file = 'output.sav',Post_Prob_Cloud, Sfc_Idx, Phase_True conf_clear_prob_clear_thresh = make_array(7,/FLOAT,VALUE=MISSING) prob_cloud_conf_cloud_thresh = make_array(7,/FLOAT,VALUE=MISSING) for Isfc = 1,7 do begin idx = where(sfc_idx eq isfc and Post_Prob_Cloud ge 0.0,cc) if (cc gt 0) then begin idx_clear = where(Post_Prob_Cloud[idx] lt 0.5,cc_clear) if (cc_clear gt 1) then begin pick_point, Post_Prob_Cloud[idx[idx_clear]],conf_prob_thresh, /clear ;conf_clear_prob_clear_thresh[Isfc-1] = max([0.10,conf_prob_thresh]) conf_clear_prob_clear_thresh[Isfc-1] = max([0.05,conf_prob_thresh]) endif idx_cloudy = where(Post_Prob_Cloud[idx] ge 0.5,cc_cloudy) if (cc_cloudy gt 1) then begin pick_point, Post_Prob_Cloud[idx[idx_cloudy]], conf_prob_thresh, /cloudy ;prob_cloud_conf_cloud_thresh[Isfc-1] = min([0.90,conf_prob_thresh]) prob_cloud_conf_cloud_thresh[Isfc-1] = min([0.95,conf_prob_thresh]) endif endif print, "conf/prob threshes = ", isfc, conf_clear_prob_clear_thresh[Isfc-1], prob_cloud_conf_cloud_thresh[Isfc-1] endfor ;---------------------------------------------------------------------- ; write these to the table ;---------------------------------------------------------------------- sd_id = ncdf_open(nc_filename,/WRITE) ncdf_attput,sd_id, "prior_name", Prior_Name, /CHAR, /global sds_id = ncdf_varid(sd_id, 'conf_clear_prob_clear_thresh') ncdf_varget, sd_id, sds_id, Conf_Clear_Prob_Clear_Thresh_Old ncdf_varput, sd_id, sds_id, Conf_Clear_Prob_Clear_Thresh sds_id = ncdf_varid(sd_id, 'prob_cloud_conf_cloud_thresh') ncdf_varget, sd_id, sds_id, Prob_Cloud_Conf_Cloud_Thresh_Old ncdf_varput, sd_id, sds_id, Prob_Cloud_Conf_Cloud_Thresh ncdf_close, sd_id print, 'Old Conf Clear/Prob Clear = ', Conf_Clear_Prob_Clear_Thresh_Old print, 'New Conf Clear/Prob Clear = ', Conf_Clear_Prob_Clear_Thresh print, 'Old Prob Cloud / Conf Cloud = ', Prob_Cloud_Conf_Cloud_Thresh_Old print, 'New Prob Cloud / Conf Cloud = ', Prob_Cloud_Conf_Cloud_Thresh endif ;==================================================================== ; make table of performance ;==================================================================== surface_name = ["all","ocean", "other water", "land", "snow", "arctic", "antarctic", "desert"] Sfc_Idx = l2.bayes_mask_sfc_type Metric_Clear_Table = fltarr(3,Nsfc+1) Cld_Phase_Acc_Table = fltarr(3,Nsfc+1) Cld_Frac_True_Table = fltarr(3,Nsfc+1) Cld_Frac_Table = fltarr(3,Nsfc+1) Solzen_Min_Temp = [0,0,90] Solzen_Max_Temp = [180,90,180] print, '===============================================================================' printf, lun, '===============================================================================' print, "PC Table" printf, lun, "PC Table" print, '===============================================================================' printf, lun, '===============================================================================' print, "Cod Thresh = ", Cod_Val_Min_Thresh printf, lun, "Cod Thresh = ", Cod_Val_Min_Thresh print, "surface / clear (a,d,n) / phase (a,d,n) " print, lun, "surface / clear (a,d,n) / phase (a,d,n) " for Isfc = 0,Nsfc do begin for iday = 0, 2 do begin if (isfc eq 0) then begin idx = where (Sfc_Idx ge 1 and $ Phase_True ge 0 and $ (abs(Cod_Lidar) le 0.01 or Cod_Lidar gt Cod_Val_Min_Thresh) and $ L2.solar_zenith_angle ge Solzen_Min_Temp[iday] and $ L2.solar_zenith_angle le Solzen_Max_Temp[iday],cc) endif if (isfc gt 0) then begin idx = where (Sfc_Idx eq isfc and $ Phase_True ge 0 and $ (abs(Cod_Lidar) le 0.01 or Cod_Lidar gt Cod_Val_Min_Thresh) and $ L2.solar_zenith_angle ge Solzen_Min_Temp[iday] and $ L2.solar_zenith_angle le Solzen_Max_Temp[iday],cc) endif make_acc_metrics, phase_true[idx], Post_Prob_Cloud[idx], $ Acc_Temp,Bacc_Temp,Wacc_Temp,Bwacc_Temp, Bce_Temp, $ Tpr_Temp,Tnr_Temp,Fpr_Temp,Fnr_Temp if (Report_Metric eq 'ACC') then Metric_Clear = Acc_Temp if (Report_Metric eq 'BACC') then Metric_Clear = Bacc_Temp if (Report_Metric eq 'WACC') then Metric_Clear = Wacc_Temp if (Report_Metric eq 'BWACC') then Metric_Clear = Bwacc_Temp if (Report_Metric eq 'BCE') then Metric_Clear = Bce_Temp Metric_Clear_Table[iday,isfc] = Metric_Clear make_phase_metrics, phase_True[idx], cloud_phase[idx], $ Acc_Phase,Bacc_Phase,Wacc_Phase,Bwacc_Phase, phase_opt, Acc_Cld_Phase, Bacc_Cld_Phase Cld_Phase_Acc_Table[iday,isfc] = Acc_Cld_Phase make_cld_frac_metric, phase_true[idx], Post_Prob_Cloud[idx], Cld_Frac_True_Temp, Cld_Frac_Ret_Temp Cld_Frac_True_Table[iday,isfc] = Cld_Frac_True_Temp Cld_Frac_Table[iday,isfc] = Cld_Frac_Ret_Temp endfor print, format="(a10,a3, 3F6.2, a3, 3F6.2)", surface_name[isfc], ' / ',Metric_Clear_Table[*,isfc] printf, lun, format="(a10,a3, 3F6.2, a3, 3F6.2)", surface_name[isfc], ' / ',Metric_Clear_Table[*,isfc] endfor ;-------------------------------------------------------------------------------------------------- ; Output Cld Frac Table ;-------------------------------------------------------------------------------------------------- print, '===============================================================================' printf, lun, '===============================================================================' print, "PM and Cld Frac Table" printf, lun, "PM and Cld Frac Table" print, '===============================================================================' printf, lun, '===============================================================================' fmt_string = "(a12,a1,F6.2,a1,F6.2,a1,F6.2,a1,F6.2,a1,F6.2,a1,F6.2,a1,F6.2,a1,F6.2,a1,F6.2)" print, ' ,All ,Day ,Night' print, 'surface , PM, TF, EF, PM, TF, EF, PM, TF, EF' printf,lun, ' ,All ,Day ,Night' printf,lun, 'surface , PM, TF, EF, PM, TF, EF, PM, TF, EF' for Isfc = 0,Nsfc do begin print, format=fmt_string, surface_name[isfc],', ',Metric_Clear_Table[0,isfc], $ ', ',Cld_Frac_True_Table[0,isfc], $ ', ',Cld_Frac_Table[0,isfc], $ ', ',Metric_Clear_Table[1,isfc], $ ', ',Cld_Frac_True_Table[1,isfc], $ ', ',Cld_Frac_True_Table[1,isfc], $ ', ',Metric_Clear_Table[2,isfc], $ ', ',Cld_Frac_True_Table[2,isfc], $ ', ',Cld_Frac_Table[2,isfc] printf,lun, format=fmt_string, surface_name[isfc],', ',Metric_Clear_Table[0,isfc], $ ', ',Cld_Frac_True_Table[0,isfc], $ ', ',Cld_Frac_Table[0,isfc], $ ', ',Metric_Clear_Table[1,isfc], $ ', ',Cld_Frac_True_Table[1,isfc], $ ', ',Cld_Frac_True_Table[1,isfc], $ ', ',Metric_Clear_Table[2,isfc], $ ', ',Cld_Frac_True_Table[2,isfc], $ ', ',Cld_Frac_Table[2,isfc] endfor ;------------------------------------------------------------------------------------------------- ; output cloud phase performance ;------------------------------------------------------------------------------------------------- print, '===============================================================================' printf,lun, '===============================================================================' print, "Phase Accuracy Table" printf, lun, "Phase Accuracy Table" print, '===============================================================================' printf,lun, '===============================================================================' fmt_string = "(a12,a1,F6.2,a1,F6.2,a1,F6.2)" print, ' ,All ,Day ,Night' printf, lun, ' ,All ,Day ,Night' for Isfc = 0,Nsfc do begin print, format=fmt_string, surface_name[isfc],', ',Cld_Phase_Acc_Table[0,isfc], $ ', ',Cld_Phase_Acc_Table[1,isfc], $ ', ',Cld_Phase_Acc_Table[2,isfc] printf,lun, format=fmt_string, surface_name[isfc],', ',Cld_Phase_Acc_Table[0,isfc], $ ', ',Cld_Phase_Acc_Table[1,isfc], $ ', ',Cld_Phase_Acc_Table[2,isfc] endfor ;-------------------------------------------------------------------------------------------------- ; compute cod distribution in clear and hit rate as a function of cod ;-------------------------------------------------------------------------------------------------- if (~keyword_set(no_plot)) then begin plot_clear_cod_and_hit_rate, Missing, Sfc_Idx, Cloud_Mask_Binary, Phase_Lidar, Cod_Lidar, Nsfc endif ;====================================================================================== ; BEGIN OPTIMIZE SECTION ;====================================================================================== if (~keyword_set(optimize)) then begin close, lun free_lun, lun goto, skip_opt endif printf, lun, 'Opt Metric = ', Opt_Metric printf, lun, 'Report Metric = ', Report_Metric ;====================================================================================== ; rank classifiers ;====================================================================================== Ranked_Class_Idx = make_array(Nclass,Nsfc+1,/INT,VALUE=0) Class_Idx_Optimal = make_array(Nclass,Nsfc+1,/INT,VALUE=-1) Metric_Var_Clear_Opt = make_array(Nclass,Nsfc+1,/Float,VALUE=missing) Opt_Var_Clear_Opt = make_array(Nclass,Nsfc+1,/Float,VALUE=missing) Cld_Frac_Opt = make_array(Nclass,Nsfc+1,/Float,VALUE=missing) Nclass_Optimal = make_array(Nsfc+1,/Int,VALUE=0) Nsfc_opt = Nsfc if (keyword_set(global)) then Nsfc_opt = 0 for Sfc_Idx_Optimal = 0, Nsfc_opt do begin if (Sfc_Idx_Optimal eq 0) then begin idx = where(Phase_True ge 0,cc) endif else begin idx = where(Sfc_Idx_Optimal eq fix(l2.bayes_mask_sfc_type) and Phase_True ge 0,cc) endelse for Class_Idx = 0, Nclass-1 do begin compute_best_class, Classifier_Names, Nclass, Nclass_Optimal[Sfc_Idx_Optimal], Class_Idx_Optimal[*,Sfc_Idx_Optimal], $ Clear_Cond_Ratio[idx,*], Ice_Cond_Ratio[idx,*], Water_Cond_Ratio[idx,*], $ Prior_Prob_Clear[idx], Prior_Prob_Ice[idx], Prior_Prob_Water[idx], $ Phase_True[idx], Class_Idx_New, Opt_Var_New, Cld_Frac_New, $ Cod_Lidar, cod_clear_water_thresh, cod_clear_ice_thresh,Cod_Val_Min_Thresh, $ Opt_Metric, Metric_Thresh, Report_Metric, Report_Var_New, Use_Phase, NClass_Optimal_Thresh, Phase_Opt if (Class_Idx_New ge 0) then begin Class_Idx_Optimal[Class_Idx,Sfc_Idx_Optimal] = Class_Idx_New Metric_Var_Clear_Opt[Nclass_Optimal[Sfc_Idx_Optimal],Sfc_Idx_Optimal] = Report_Var_New Opt_Var_Clear_Opt[Nclass_Optimal[Sfc_Idx_Optimal],Sfc_Idx_Optimal] = Opt_Var_New Cld_Frac_Opt[Nclass_Optimal[Sfc_Idx_Optimal],Sfc_Idx_Optimal] = Cld_Frac_New Nclass_Optimal[Sfc_Idx_Optimal] = Nclass_Optimal[Sfc_Idx_Optimal] + 1 endif endfor print, '=============== Optimal Results Prioritized ===========' printf, lun, '=============== Optimal Results Prioritized ===========' print, "Optimal Classifiers for Sfc_Idx = ", Sfc_Idx_Optimal printf, lun, "Optimal Classifiers for Sfc_Idx = ", Sfc_Idx_Optimal print, "Opt_Rank, Class_Idx, Class_Name, Pc_Clear_opt, Cld_Frac_Opt" printf, lun, "Opt_Rank, Class_Idx, Class_Name, Pc_Clear_opt, Cld_Frac_Opt" xname_Temp = make_array(Nclass_Optimal[Sfc_Idx_Optimal],/String,VALUE="") Opt_Var_Temp = make_array(Nclass_Optimal[Sfc_Idx_Optimal],/Float,VALUE=missing) Cld_Frac_Temp = make_array(Nclass_Optimal[Sfc_Idx_Optimal],/Float,VALUE=missing) for Class_Idx = 0,NClass_Optimal[Sfc_Idx_Optimal] - 1 do begin print, format = "(I3,I3,a25,4f9.4)",Class_Idx, Class_Idx_Optimal[Class_Idx,Sfc_Idx_Optimal], $ Classifier_Names[Class_Idx_Optimal[Class_Idx,Sfc_Idx_Optimal]], $ Opt_Var_Clear_Opt[Class_Idx,Sfc_Idx_Optimal], Cld_Frac_Opt[Class_Idx,Sfc_Idx_Optimal] printf, lun, format = "(I3,I3,a25,5f9.4)",Class_Idx, Class_Idx_Optimal[Class_Idx,Sfc_Idx_Optimal], $ Classifier_Names[Class_Idx_Optimal[Class_Idx,Sfc_Idx_Optimal]], $ Opt_Var_Clear_Opt[Class_Idx,Sfc_Idx_Optimal], Cld_Frac_Opt[Class_Idx,Sfc_Idx_Optimal] xname_Temp[Class_Idx] = Classifier_Names[Class_Idx_Optimal[Class_Idx,Sfc_Idx_Optimal]] Opt_Var_Temp[Class_Idx] = Opt_Var_Clear_Opt[Class_Idx,Sfc_Idx_Optimal] Cld_Frac_Temp[Class_Idx] = Cld_Frac_Opt[Class_Idx,Sfc_Idx_Optimal] endfor ;--- plot prioritized classifiers this surface if (~keyword_set(no_plot)) then begin p20 = plot(Opt_Var_Temp, dim=[600,400], color='black', background_color='white', margin=[0.20,0.25,0.1,0.05], $ xtitle = 'Classifier', xtext_orientation = 30, xtickname=xname_Temp, $ ytitle = 'Performance Metric', name = 'PC', symbol="circle") t20 = text(0.2,0.8,/normal,Sfc_Type_Name[Sfc_Idx_Optimal]) p20.save, 'pc_vs_class_'+sfc_type_name[sfc_idx_optimal]+'.png', resolution=200 p20.close p22 = plot(Cld_Frac_Temp,dim=[600,400],color='black',background_color='white',margin=[0.20,0.25,0.1,0.05],$ xtitle = 'Classifier', xtext_orientation = 30, xtickname=xname_Temp, $ ytitle = 'Cloud Fraction', name = 'Cld_Frac', yrange = [0,1], symbol="circle") xxx = replicate(Cld_Frac_True[Sfc_Idx_Optimal],Cld_Frac_Temp.length) p221 = plot(xxx,/overplot,linestyle='-') t22 = text(0.2,0.8,/normal,Sfc_Type_Name[Sfc_Idx_Optimal]) p22.save, 'cldfrac_vs_class_'+sfc_type_name[sfc_idx_optimal]+'.png', resolution=200 p22.close endif endfor print, '=============== Optimal Results UnPrioritized ===============' print, format = "(a10,a33,8I3)", "Class_Idx","Test Name", indgen(Nsfc+1) printf, lun, '=============== Optimal Results UnPrioritized ===============' printf, lun, format = "(a10,a33,8I3)", "Class_Idx","Test Name", indgen(Nsfc+1) remove_class_flag = make_array(Nclass,/INT,VALUE=0) sfc_type_on_flag = make_array(Nsfc+1,Nclass,/INT,VALUE=0) for Class_Idx = 0,NClass - 1 do begin temp = make_array(Nsfc+1,/INT,Value=-1) if (keyword_set(global)) then begin for Sfc_Idx = 0, Nsfc do begin idx = where(Class_Idx eq Class_Idx_Optimal[*,0],cc) if (cc gt 0) then begin temp[0] = idx[0] sfc_type_on_flag[*,Class_Idx] = 1 endif endfor endif if (~keyword_set(global)) then begin for Sfc_Idx = 0, Nsfc do begin idx = where(Class_Idx eq Class_Idx_Optimal[*,Sfc_Idx],cc) if (cc gt 0) then begin temp[Sfc_Idx] = idx[0] sfc_type_on_flag[Sfc_Idx,Class_Idx] = 1 endif endfor endif if (~keyword_set(global)) then begin if (max(temp[1:Nsfc]) lt 0) then remove_class_flag[Class_Idx] = 1 ; skip sfc = 0 which is global endif if (keyword_set(global)) then begin if (max(temp[0]) lt 0) then remove_class_flag[Class_Idx] = 1 ; skip sfc = 0 which is global endif temp_string = string(temp,format="(I3)") idx = where(temp eq -1,cc) if (cc gt 0) then temp_string[idx] = " " print, format = "(I3,a40,8a3)", Class_Idx, Classifier_Names[Class_Idx], temp_string printf, lun, format = "(I3,a40,8a3)", Class_Idx, Classifier_Names[Class_Idx], temp_string endfor ;----- optimized peformance print, '=================================================================' print, 'final stats' print, "Report_Metric = ", Report_Metric print, "Opt_Metric = ", Opt_Metric print, '=================================================================' printf, lun, '=================================================================' printf, lun, 'final stats' printf, lun, "Report_Metric = ", Report_Metric printf, lun, "Opt_Metric = ", Opt_Metric printf, lun, '=================================================================' print, " Sfc_Idx, Report_Metric, Report_Metric After, Opt_Metric After" printf, lun," Sfc_Idx, Report_Metric, Report_Metric After, Opt_Metric After" for Sfc_Idx = 0, Nsfc do begin print, Sfc_Idx, Metric_Clear_By_Sfc[Sfc_Idx],Metric_Var_Clear_Opt[Nclass_Optimal[Sfc_Idx]-1,Sfc_Idx], Opt_Var_Clear_Opt[Nclass_Optimal[Sfc_Idx]-1,Sfc_Idx] printf, lun, Sfc_Idx, Metric_Clear_By_Sfc[Sfc_Idx], Metric_Var_Clear_Opt[Nclass_Optimal[Sfc_Idx]-1,Sfc_Idx], Opt_Var_Clear_Opt[Nclass_Optimal[Sfc_Idx]-1,Sfc_Idx] endfor close, lun free_lun, lun stop ;-------------------------------------------------------------------------------------------- ; trim table list and make script to make optimal table ;-------------------------------------------------------------------------------------------- openw, lun, 'make_optimal_tables.bat', /get_lun printf, lun, '.rnew prob_mask_phase_1d' printf, lun, '.rnew prob_mask_phase_2d' printf, lun, '.rnew prob_mask_phase_3d' Sfc_Idx_Min = 1 Sfc_Idx_Max = Nsfc - 1 if (keyword_set(global)) then begin Sfc_Idx_Min = 0 Sfc_Idx_Max = 0 endif for Class_Idx = 0, NClass - 1 do begin if (remove_class_flag[Class_Idx] eq 0) then begin sfc_flag_on_string = "[" for Sfc_Idx = 1, Nsfc-1 do begin sfc_flag_on_string = sfc_flag_on_string + string(sfc_type_on_flag[Sfc_Idx,Class_Idx],format="(I1)")+", " endfor sfc_flag_on_string = sfc_flag_on_string + string(sfc_type_on_flag[Nsfc,Class_Idx],format="(I1)")+"]" single_table_filename = "./single_tables/"+sensor+"_all_"+Classifier_Names[Class_Idx]+".nc" temp = reform(sfc_type_on_flag[1:NSfc,Class_Idx]) if (Class_Rank[Class_Idx] eq 1) then begin modify_on_flag, 1, single_table_filename, temp print, 'prob_mask_phase_1d, sensor="abhi", /new, /noplot, test_name="', + Classifier_Names[Class_Idx], + '", sfc_type_on_flag = '+sfc_flag_on_string printf, lun, format = "(200a)", 'prob_mask_phase_1d, sensor="abhi", /new, /noplot, test_name="' + Classifier_Names[Class_Idx] + '", sfc_type_on_flag = '+sfc_flag_on_string endif if (Class_Rank[Class_Idx] eq 2) then begin modify_on_flag, 2, single_table_filename, temp print, 'prob_mask_phase_2d, sensor="abhi", /new, /noplot, test_name="', + Classifier_Names[Class_Idx], + '", sfc_type_on_flag = '+sfc_flag_on_string printf, lun, format = "(200a)", 'prob_mask_phase_2d, sensor="abhi", /new, /noplot, test_name="' + Classifier_Names[Class_Idx] + '", sfc_type_on_flag = '+sfc_flag_on_string endif if (Class_Rank[Class_Idx] eq 3) then begin modify_on_flag, 3, single_table_filename, temp print, 'prob_mask_phase_3d, sensor="abhi", /new, /noplot, test_name="', + Classifier_Names[Class_Idx], + '", sfc_type_on_flag = '+sfc_flag_on_string printf, lun, format = "(200a)", 'prob_mask_phase_3d, sensor="abhi", /new, /noplot, test_name="' + Classifier_Names[Class_Idx] + '", sfc_type_on_flag = '+sfc_flag_on_string endif endif endfor close, lun free_lun, lun for Class_Idx = 0,NClass - 1 do begin if (remove_class_flag[Class_Idx] eq 1) then begin single_table_filename = "./single_tables/"+sensor+"_all_"+Classifier_Names[Class_Idx]+".nc" print, "removing ", single_table_filename spawn, "rm "+single_table_filename endif endfor skip_opt: end ;============== ; END OF MAIN ;============== ;================================================================================================== ; Class_Idx_Optimal = Classes already chosen to be optimal (in order) ; Class_Idx_Optimal_New = the new class to be added to the optimal list ; Sfc_Type = 1-7 ;================================================================================================== pro compute_best_class, Classifier_Names, Nclass, Nclass_Optimal, Class_Idx_Optimal, $ Clear_Cond_Ratio, Ice_Cond_Ratio, Water_Cond_Ratio, $ Prior_Prob_Clear, Prior_Prob_Ice, Prior_Prob_Water, $ Phase_true, Class_Idx_New, Opt_Var_Opt, Cld_Frac_Opt, $ Cod_Lidar, cod_clear_water_thresh, cod_clear_ice_thresh,Cod_Val_Min_Thresh, $ Opt_Metric, Metric_Thresh, Report_Metric, Report_Var_Opt, Use_Phase, NClass_Optimal_Thresh, Phase_Opt Nx = n_elements(Prior_Prob_Clear) Class_Idx_New = -1 if (Nclass_Optimal eq 0) then begin ACC_Clear_Base = 0.0 BACC_Clear_Base = 0.0 WACC_Clear_Base = 0.0 BWACC_Clear_Base = 0.0 BCE_Clear_Base = -100.0 ; this is really -1.0* bce. so bigger is better. number should approach 0 ACC_Phase_Base = 0.0 BACC_Phase_Base = 0.0 WACC_Phase_Base = 0.0 BWACC_Phase_Base = 0.0 endif R_Clear_Base = 1.0 R_Ice_Base = 1.0 R_Water_Base = 1.0 if (Nclass_Optimal gt 0) then begin ;------------------------------------------ ; make value for known optimal classes ;------------------------------------------ For Class_Idx = 0,Nclass-1 do begin idx = where(Class_Idx eq Class_Idx_Optimal,cc) if (cc eq 0) then goto, skip R_Clear_Base = R_Clear_Base * Clear_Cond_Ratio[*,Class_Idx] R_Ice_Base = R_Ice_Base * Ice_Cond_Ratio[*,Class_Idx] R_Water_Base = R_Water_Base * Water_Cond_Ratio[*,Class_Idx] skip: endfor ;--------------------------------------------------------------------------------- ; compute base performance ;--------------------------------------------------------------------------------- Post_Prob_Clear_Base = 1.0 / (1.0 + R_Clear_Base/Prior_Prob_Clear - R_Clear_Base) Post_Prob_Ice_Base = 1.0 / (1.0 + R_Ice_Base/Prior_Prob_Ice - R_Ice_Base) Post_Prob_Water_Base = 1.0 / (1.0 + R_Water_Base/Prior_Prob_Water - R_Water_Base) Post_Sum = Post_Prob_Clear_Base + Post_Prob_Ice_Base + Post_Prob_Water_Base Post_Prob_Clear_Base = Post_Prob_Clear_Base / Post_Sum Post_Prob_Ice_Base = Post_Prob_Ice_Base / Post_Sum Post_Prob_Water_Base = Post_Prob_Water_Base / Post_Sum Post_Prob_Cloud_Base = (1.0-Post_Prob_Clear_Base) Cloud_Phase_Base = make_array(Nx,/Int,Value=0) Cloud_Phase_Base_Uncer = make_array(Nx,/Float,Value=missing) compute_cloud_phase_vector, Post_Prob_Clear_Base, Post_Prob_Water_Base, Post_Prob_Ice_Base, Cloud_Phase_Base, Cloud_Phase_Base_Uncer idx = where (Phase_True eq 0 or Cod_Lidar gt Cod_Val_Min_Thresh,cc) make_acc_metrics, Phase_true[idx], Post_Prob_Cloud_Base[idx], Acc_Clear_Base, $ Bacc_Clear_Base,Wacc_Clear_Base,Bwacc_Clear_Base,Bce_Clear_Base, $ Tpr_Base, Tnr_Base, Fpr_Base, Fnr_Base make_phase_metrics, Phase_True[idx], Cloud_Phase_Base[idx], Acc_Phase_Base,Bacc_Phase_Base,Wacc_Phase_Base,Bwacc_Phase_Base, phase_opt, Acc_Cld_Phase, Bacc_Cld_Phase make_cld_frac_metric, Phase_True[idx], Post_Prob_Cloud_Base[idx], Cld_Frac_True_Base, Cld_Frac_True_Temp endif ;------------------------------------------------------------------------------------ ; incrementally add in other classisifiers and look for improvement ;------------------------------------------------------------------------------------ Cloud_Phase_Temp = make_array(Nx,/Int,Value=0) Cld_Frac_Opt = 0.0 Class_New_Idx = -1 if (use_phase eq 0) then begin if (Opt_Metric eq 'ACC') then Opt_Var_Opt = Acc_Clear_Base if (Opt_Metric eq 'BACC') then Opt_Var_Opt = Bacc_Clear_Base if (Opt_Metric eq 'WACC') then Opt_Var_Opt = Wacc_Clear_Base if (Opt_Metric eq 'BWACC') then Opt_Var_Opt = Bwacc_Clear_Base if (Opt_Metric eq 'BCE') then Opt_Var_Opt = Bce_Clear_Base endif if (use_phase eq 1) then begin if (Opt_Metric eq 'ACC') then Opt_Var_Opt = Acc_Phase_Base if (Opt_Metric eq 'BACC') then Opt_Var_Opt = Bacc_Phase_Base if (Opt_Metric eq 'WACC') then Opt_Var_Opt = Wacc_Phase_Base if (Opt_Metric eq 'BWACC') then Opt_Var_Opt = Bwacc_Phase_Base endif ;print, 'Starting Base = ', Opt_Var_Opt For Class_Idx = 0,Nclass-1 do begin idx = where(Class_Idx eq Class_Idx_Optimal,cc) if (cc ne 0) then goto, skip2 R_Clear_Temp = R_Clear_Base * Clear_Cond_Ratio[*,Class_Idx] R_Ice_Temp = R_Ice_Base * Ice_Cond_Ratio[*,Class_Idx] R_Water_Temp = R_Water_Base * Water_Cond_Ratio[*,Class_Idx] Post_Prob_Clear_Temp = 1.0 / (1.0 + R_Clear_Temp/Prior_Prob_Clear - R_Clear_Temp) Post_Prob_Ice_Temp = 1.0 / (1.0 + R_Ice_Temp/Prior_Prob_Ice - R_Ice_Temp) Post_Prob_Water_Temp = 1.0 / (1.0 + R_Water_Temp/Prior_Prob_Water - R_Water_Temp) Post_Sum = Post_Prob_Clear_Temp + Post_Prob_Ice_Temp + Post_Prob_Water_Temp Post_Prob_Clear_Temp = Post_Prob_Clear_Temp / Post_Sum Post_Prob_Ice_Temp = Post_Prob_Ice_Temp / Post_Sum Post_Prob_Water_Temp = Post_Prob_Water_Temp / Post_Sum Post_Prob_Cloud_Temp = 1.0 - Post_Prob_Clear_Temp Cloud_Phase_Temp = make_array(Nx,/Int,Value=0) Cloud_Phase_Temp_Uncer = make_array(Nx,/Float,Value=missing) compute_cloud_phase_vector, Post_Prob_Clear_Temp, Post_Prob_Water_Temp, Post_Prob_Ice_Temp, $ Cloud_Phase_Temp, Cloud_Phase_Temp_Uncer idx = where (Phase_True eq 0 or Cod_Lidar gt Cod_Val_Min_Thresh,cc) make_acc_metrics, Phase_true[idx], Post_Prob_Cloud_Temp[idx], Acc_Clear_Temp, $ Bacc_Clear_Temp,Wacc_Clear_Temp,Bwacc_Clear_Temp,Bce_Clear_Temp, $ Tpr_Temp, Tnr_Temp, Fpr_Temp, Fnr_Temp make_phase_metrics, Phase_True[idx], Cloud_Phase_Temp[idx], Acc_Phase_Temp,Bacc_Phase_Temp,Wacc_Phase_Temp,Bwacc_Phase_Temp, phase_opt, $ Acc_Cld_Phase, Bacc_Cld_Phase make_cld_frac_metric, Phase_True[idx], Post_Prob_Cloud_Temp[idx], Cld_Frac_True_Temp, Cld_Frac_Ret_Temp ; print, 'acc metrics ', Acc_Clear_Temp, Acc_Phase_Temp, total(Post_Prob_Cloud_Temp[idx]), $ ; total(Post_Prob_Water_Temp[idx]), $ ; total(Post_Prob_Ice_Temp[idx]), $ ; total(Cloud_Phase_Temp[idx]) ;--- test for a new optimal state optimal = 0 if (use_phase eq 0) then begin if (Opt_Metric eq 'ACC') then Opt_Var_Temp = Acc_Clear_Temp if (Opt_Metric eq 'BACC') then Opt_Var_Temp = Bacc_Clear_Temp if (Opt_Metric eq 'WACC') then Opt_Var_Temp = Wacc_Clear_Temp if (Opt_Metric eq 'BWACC') then Opt_Var_Temp = Bwacc_Clear_Temp if (Opt_Metric eq 'BCE') then Opt_Var_Temp = Bce_Clear_Temp if (Report_Metric eq 'ACC') then Report_Var_Temp = Acc_Clear_Temp if (Report_Metric eq 'BACC') then Report_Var_Temp = Bacc_Clear_Temp if (Report_Metric eq 'WACC') then Report_Var_Temp = Wacc_Clear_Temp if (Report_Metric eq 'BWACC') then Report_Var_Temp = Bwacc_Clear_Temp if (Report_Metric eq 'BCE') then Report_Var_Temp = Bce_Clear_Temp endif if (use_phase eq 1) then begin if (Opt_Metric eq 'ACC') then Opt_Var_Temp = Acc_Phase_Temp if (Opt_Metric eq 'BACC') then Opt_Var_Temp = Bacc_Phase_Temp if (Opt_Metric eq 'WACC') then Opt_Var_Temp = Wacc_Phase_Temp if (Opt_Metric eq 'BWACC') then Opt_Var_Temp = Bwacc_Phase_Temp ;if (Opt_Metric eq 'BCE') then Opt_Var_Temp = Bce_Phase_Temp if (Report_Metric eq 'ACC') then Report_Var_Temp = Acc_Phase_Temp if (Report_Metric eq 'BACC') then Report_Var_Temp = Bacc_Phase_Temp if (Report_Metric eq 'WACC') then Report_Var_Temp = Wacc_Phase_Temp if (Report_Metric eq 'BWACC') then Report_Var_Temp = Bwacc_Phase_Temp ;if (Report_Metric eq 'BCE') then Report_Var_Temp = Bce_Phase_Temp endif ;--- take any improvement until you reach the minimum of classes if ((Opt_Var_Temp - Opt_Var_Opt ge 0.0) and (NClass_Optimal le NClass_Optimal_Thresh)) then optimal = 1 ;--- require a minimum improvement if (Opt_Var_Temp - Opt_Var_Opt ge Metric_Thresh) then optimal = 1 if (optimal eq 1) then begin Class_Idx_New = Class_Idx Report_Var_Opt = Report_Var_Temp Opt_Var_Opt = Opt_Var_Temp Cld_Frac_Opt = Cld_Frac_Ret_Temp endif skip2: endfor end ;----------------------------------------------------------------------------- ; pp = post prob ; conf_prob_thresh = output threshold of prob boundary between conf and prob ;----------------------------------------------------------------------------- pro pick_point,pp,conf_prob_thresh, clear=clear,cloudy=cloudy nbins = 500 if (keyword_set(clear)) then begin binmin = 0.0 binmax = 0.5 cdf_thresh = 0.75 endif if (keyword_set(cloudy)) then begin binmin = 0.5 binmax = 1.0 cdf_thresh = 0.25 endif binsize = (binmax-binmin) / nbins x = binmin + findgen(nbins)*binsize + binsize/2.0 histo = histogram(pp,min=binmin,binsize=binsize,nbins=nbins) z = total(histo, /cumulative) / total(histo) idx_sub = where(z gt cdf_thresh,cc) if (cc le 1) then print, 'pick point error' isub = min([nbins-2,max([1,idx_sub[0]])]) conf_prob_thresh = x[isub] end ;------------------------------------------------------------------------------------------------------- ; compute the cloud phase from Post_Probs ; ; 0 = clear ; 1 = water ; 2 = ice ;------------------------------------------------------------------------------------------------------- pro compute_cloud_phase, Post_Prob_Clear, Post_Prob_Water, Post_Prob_Ice, Cloud_Phase, Cloud_Phase_Uncer Post_Prob_Cloud = 1.0 - Post_Prob_Clear Post_Prob_Water_Norm = Post_Prob_Water / Post_Prob_Cloud Post_Prob_Ice_Norm = Post_Prob_Ice / Post_Prob_Cloud ;--- make cloud phase and its uncertainty ;--- initialize as clear Cloud_Phase = 0 Cloud_Phase_Uncer = Post_Prob_Clear if (Post_Prob_Clear lt 0.50) then begin Cloud_Phase = 1 Cloud_Phase_Uncer = 1.0 - Post_Prob_Water_Norm if (Post_Prob_Ice_Norm gt Post_Prob_Water_Norm) then begin Cloud_Phase = 2 Cloud_Phase_Uncer = 1.0 - Post_Prob_Ice_Norm endif endif end pro compute_cloud_phase_vector, Post_Prob_Clear, Post_Prob_Water, Post_Prob_Ice, Cloud_Phase, Cloud_Phase_Uncer Post_Prob_Cloud = 1.0 - Post_Prob_Clear Post_Prob_Water_Norm = Post_Prob_Water / Post_Prob_Cloud Post_Prob_Ice_Norm = Post_Prob_Ice / Post_Prob_Cloud ;--- make cloud phase and its uncertainty ;--- initialize as clear Cloud_Phase[*] = 0 Cloud_Phase_Uncer = Post_Prob_Clear idx = where(Post_Prob_Clear lt 0.50,cc) if (cc gt 0) then begin Cloud_Phase[idx] = 1 Cloud_Phase_Uncer[idx]= 1.0 - Post_Prob_Water_Norm[idx] endif idx = where(Post_Prob_Clear lt 0.50 and Post_Prob_Ice_Norm gt Post_Prob_Water_Norm, cc) if (cc gt 0) then begin Cloud_Phase[idx] = 2 Cloud_Phase_Uncer[idx] = 1.0 - Post_Prob_Ice_Norm[idx] endif end