;$Id$ @apply_filters @write_1d_table @fmft @ndsi @compute_class_cond @make_true_phase @ecm_sfc_type ;-------------------------------------------------------------------------------------------- ; probabilistic phase for 1d naive bayesian ; ;-------------------------------------------------------------------------------------------- pro prob_mask_phase_1d, $ day = day, night=night, sensor_name = sensor_name, noplot=noplot, new=new, old=old, $ test_name = test_name, denis=denis, sfc_type_on_flag_input = sfc_type_on_flag_input, $ path_dir = path_dir, nomod=nomod, Time_Diff_Max = Time_Diff_Max, $ cod_clear_ice_thresh = cod_clear_ice_thresh, cod_clear_water_thresh = cod_clear_water_thresh, yesplot=yesplot, isccp=isccp if (~keyword_set(sensor_name)) then sensor_name='myd02ssh' if (~keyword_set(path_dir)) then path_dir='./' if (~keyword_set(Time_Diff_Max)) then Time_Diff_Max = 10.0 ; minutes if (~keyword_set(cod_clear_ice_thresh)) then cod_clear_ice_thresh = 0.1 if (~keyword_set(cod_clear_water_thresh)) then cod_clear_water_thresh = 0.1 if (sensor_name eq 'abi') then training_source = 'ecm_training_abi_train.sav' if (sensor_name eq 'abhi') then training_source = 'ecm_training_abhi_train.sav' if (sensor_name eq 'avhrr') then training_source = 'ecm_training_avhrr_train.sav' if (sensor_name eq 'myd02ssh') then training_source = 'ecm_training_myd02ssh_train.sav' if (sensor_name eq 'viirs') then training_source = 'ecm_training_viirs_train_snpp_5+333_2013.sav' if (sensor_name eq 'mtsat2') then training_source = 'ecm_training_mtsat2_train.sav' if (sensor_name eq 'vgac') then training_source = 'ecm_training_vgac_train.sav' if (sensor_name eq 'seviri') then training_source = 'ecm_training_seviri_train.sav' ;--- check for sensor and test inconsistencies if (strpos(sensor_name,'viirs') lt 0 and strpos(test_name,'dnb') gt 0) then begin print, 'DNB classifiers impossible' return endif ;-- control training data modification in make_true_phase routine mod_flag = 1 if (keyword_set(nomod)) then mod_flag = 0 ;------------------------------- ; read in data ;------------------------------- training_path = './training_data/' print, 'restoring ', training_path + training_source restore, training_path + training_source f = training_data training_data = !null ntrain = f.length nb_rank = 1 nsfc = 7 if (~keyword_set(sfc_type_on_flag_input)) then sfc_type_on_flag_input = make_array(nsfc,/INT,VALUE=1) ;---- some parameters missing = -999.0 laplace_flag = 0 max_r = 1000.0 count_min = 10.0 rut_solzen_thresh = 80.0 ;--- partial cloud filter cfrac_min = 0.1 cfrac_max = 0.9 ;--- set any geographical limits lat_min = -90.0 lat_max = 90.0 ;---- adjust nasa dnb to look like noaa dnb [REMOVE WHEN ADJUSTMENT in TRAINING DATA] ;Dnb_Coef = [-0.373685,0.977945,-0.000261637] ;Andy Dnb_Coef = [-0.118767,0.962452,-0.000144502] ; Denis idx = where(f.refl_lunar_dnb_nom ne missing,cc) if (cc gt 0) then begin f[idx].refl_lunar_dnb_nom = Dnb_Coef[0] + $ Dnb_Coef[1]*f[idx].refl_lunar_dnb_nom + $ Dnb_Coef[2]*f[idx].refl_lunar_dnb_nom^2 endif ;--- product filtering cod_cloud_thresh = 1.0 etrop_cloud_thresh = 0.9 ;--- default classifier limits - only modify in classifier definition if needed zen_min = 0.0 zen_max = 90.0 solzen_min = 0.0 solzen_max = 180.0 solglintzen_min = 0.0 solglintzen_max = 180.0 solscatang_min = 0.0 solscatang_max = 180.0 solglint_mask_min = 0 solglint_mask_max = 1 lunzen_min = 0.0 lunzen_max = 180.0 lunglintzen_min = 0.0 lunglintzen_max = 180.0 lunscatang_min = 0.0 lunscatang_max = 180.0 lunglint_mask_min = 0 lunglint_mask_max = 1 coast_mask_min = 0 coast_mask_max = 1.0 city_mask_min = 0 city_mask_max = 1.0 moon_illum_frac_min = 0.0 moon_illum_frac_max = 100.0 slant_tpw_min = 0.0 slant_tpw_max = 100.0 tsfc_min = 180.0 tsfc_max = 340.0 zsfc_min = -100.0 zsfc_max = 20000.0 zsfc_std_min = 00.0 zsfc_std_max = 1000.0 snow_class_min = 1 snow_class_max = 3 sfc_type_on_flag = sfc_type_on_flag_input ;----- set solar limits day_string = 'all' solzen_min = 0.0 solzen_max = 180.0 solzen_day_max = 85.0 zsfc_std_max_iruni = 300.0 solscatang_day_min = 75.0 if (keyword_set(day)) then begin day_string = 'day' solzen_min = 0.0 solzen_max = solzen_day_max endif if (keyword_set(night)) then begin day_string = 'night' solzen_min = 90.0 solzen_max = 180.0 endif ;------------------------------------------------------------------------------ ; small positive (near zero) dnb reflectances should be missing ;------------------------------------------------------------------------------ idx = where(f.refl_lunar_dnb_nom lt 0.1,cc) & if (cc gt 0) then f[idx].refl_lunar_dnb_nom = missing ;------------------------------------------------------------------------------ ; loop through classifiers ;------------------------------------------------------------------------------ nc_header = '1D Mask and Phase LUT' print , "1D Classifier = ", test_name case test_name of "bt11ratio": begin classifier_name = 'bt11ratio' xtitle = 'bt11ratio' x = btratio(f.temp_11_0um_nom,f.surface_temperature_nwp,f.tropopause_temperature_nwp,missing) x_min = -0.1 & x_bin = 0.01 & nbins_x = 120 nchan_used = 1 & channel_wvl_used = [11000] end "bt11": begin classifier_name = 'bt11' xtitle = 'bt11' x = f.temp_11_0um_nom x_min = 160.0 & x_bin = 2.0 & nbins_x = 100 nchan_used = 1 & channel_wvl_used = [1100] end "bt133": begin classifier_name = 'bt133' xtitle = 'bt133' x = f.temp_13_3um_nom x_min = 180.0 & x_bin = 2.0 & nbins_x = 60 nchan_used = 1 & channel_wvl_used = [1330] end "refl065cv": begin classifier_name = 'refl065cv' xtitle = 'refl065cv' x = f.refl_0_65um_nom_stddev_sub / f.refl_0_65um_nom idx = where(f.refl_0_65um_nom_stddev_sub eq missing or f.refl_0_65um_nom eq missing, cc) if (cc gt 0) then x[idx] = missing x_min = 0.0 & x_bin = 0.02 & nbins_x = 50.0 nchan_used = 1 & channel_wvl_used = [650] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min ;zsfc_max = 2000.0 coast_mask_max = 0 zsfc_std_max = zsfc_std_max_iruni end "logrefl065cv": begin classifier_name = 'logrefl065cv' xtitle = 'logrefl065cv' x = f.refl_0_65um_nom_stddev_sub / f.refl_0_65um_nom idx = where(f.refl_0_65um_nom_stddev_sub eq missing or f.refl_0_65um_nom eq missing, cc) if (cc gt 0) then x[idx] = missing x_min = -2.0 & x_bin = 0.05 & nbins_x = 50.0 idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx]) idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min idx = where(x ne missing and x lt x_min) & if (cc gt 0) then x[idx] = x_min nchan_used = 1 & channel_wvl_used = [650] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min ;zsfc_max = 2000.0 coast_mask_max = 0 zsfc_std_max = zsfc_std_max_iruni end "refl065std": begin classifier_name = 'refl065std' xtitle = 'refl065std' x = f.refl_0_65um_nom_stddev_3x3 x_min = 0.0 & x_bin = 0.1 & nbins_x = 100.0 nchan_used = 1 & channel_wvl_used = [650] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min ;zsfc_max = 2000.0 coast_mask_max = 0 zsfc_std_max = zsfc_std_max_iruni end "btdclr10": begin classifier_name = 'btdclr10' xtitle = 'btdclr10' x = f.temp_10_4um_nom_clear_sky - f.temp_10_4um_nom x_min = -15 & x_bin = 0.5 & nbins_x = 100.0 nchan_used = 1 & channel_wvl_used = [10400] solzen_min = 0.0 solzen_max = 180.0 solglintzen_min = 0.0 solglintzen_max = 180.0 end "btdclr11": begin classifier_name = 'btdclr11' xtitle = 'btdclr11' x = f.temp_11_0um_nom_clear_sky - f.temp_11_0um_nom x_min = -15 & x_bin = 0.5 & nbins_x = 100.0 nchan_used = 1 & channel_wvl_used = [11000] end "btdclr11_day": begin classifier_name = 'btdclr11_day' xtitle = 'btdclr11' x = f.temp_11_0um_nom_clear_sky - f.temp_11_0um_nom x_min = -15 & x_bin = 0.5 & nbins_x = 100.0 nchan_used = 1 & channel_wvl_used = [11000] solzen_min = 0.0 solzen_max = 90.0 end "btdclr11_night": begin classifier_name = 'btdclr11_night' xtitle = 'btdclr11' x = f.temp_11_0um_nom_clear_sky - f.temp_11_0um_nom x_min = -15 & x_bin = 0.5 & nbins_x = 100.0 nchan_used = 1 & channel_wvl_used = [11000] solzen_min = 90.0 solzen_max = 180.0 end "btdclr12": begin classifier_name = 'btdclr12' xtitle = 'btdclr12' x = f.temp_12_0um_nom_clear_sky - f.temp_12_0um_nom x_min = -15 & x_bin = 0.5 & nbins_x = 100.0 nchan_used = 1 & channel_wvl_used = [12000] solzen_min = 0.0 solzen_max = 180.0 solglintzen_min = 0.0 solglintzen_max = 180.0 end "etropo11": begin classifier_name = 'etropo11' xtitle = 'etropo11' x = f.emiss_tropo_11_0um_nom x_min = -0.4 & x_bin = 0.01 & nbins_x = 150 nchan_used = 1 & channel_wvl_used = [11000] tsfc_min = 230.0 end "etropo11_day": begin classifier_name = 'etropo11_day' xtitle = 'etropo11' x = f.emiss_tropo_11_0um_nom x_min = -0.4 & x_bin = 0.01 & nbins_x = 150 nchan_used = 1 & channel_wvl_used = [11000] tsfc_min = 230.0 solzen_max = 95.0 end "etropo11_night": begin classifier_name = 'etropo11_night' xtitle = 'etropo11' x = f.emiss_tropo_11_0um_nom x_min = -0.4 & x_bin = 0.01 & nbins_x = 150 nchan_used = 1 & channel_wvl_used = [11000] tsfc_min = 230.0 solzen_min = 85.0 end "etropo10": begin classifier_name = 'etropo10' xtitle = 'etropo10' x = f.emiss_tropo_10_4um_nom x_min = -0.4 & x_bin = 0.01 & nbins_x = 150 nchan_used = 1 & channel_wvl_used = [10400] tsfc_min = 230.0 end "zopa": begin classifier_name = 'zopa' xtitle = 'zopa' x = f.cld_height_opaque x_min = -100.0 & x_bin = 200.0 & nbins_x = 80 nchan_used = 1 & channel_wvl_used = [0] end "dzopasfc": begin classifier_name = 'zopa' xtitle = 'dzopasfc' x = f.cld_height_opaque - f.surface_elevation x_min = -2000.0 & x_bin = 200.0 & nbins_x = 90 nchan_used = 1 & channel_wvl_used = [0] end "logzopa": begin classifier_name = 'logzopa' xtitle = 'logzopa' x = f.cld_height_opaque x_min = 1.5 & x_bin = 0.1 & nbins_x = 30.0 idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx]) idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min nchan_used = 1 & channel_wvl_used = [0] end "topa": begin classifier_name = 'topa' xtitle = 'topa' x = f.cld_temp_opaque x_min = 200.0 & x_bin = 2.0 & nbins_x = 60 nchan_used = 1 & channel_wvl_used = [0] end "dbt11maxsub": begin classifier_name = 'dbt11maxsub' xtitle = 'dbt11maxsub' x = f.temp_11_0um_nom_max_sub - f.temp_11_0um_nom x_min = 0.0 & x_bin = 0.2 & nbins_x = 100.0 nchan_used = 1 & channel_wvl_used = [11000] ;zsfc_max = 2000.0 coast_mask_max = 0 end "logdbt11maxminsub": begin classifier_name = 'logdbt11maxminsub' xtitle = 'logdbt11maxminsub' x = f.temp_11_0um_nom_max_sub - f.temp_11_0um_nom_min_sub x_min = -2.0 & x_bin = 0.04 & nbins_x = 100.0 idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx]) idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min nchan_used = 1 & channel_wvl_used = [11000] coast_mask_max = 0 end "logdbt11maxsub": begin classifier_name = 'logdbt11maxsub' xtitle = 'logdbt11maxsub' x = f.temp_11_0um_nom_max_sub - f.temp_11_0um_nom x_min = -2.0 & x_bin = 0.04 & nbins_x = 100.0 idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx]) idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min nchan_used = 1 & channel_wvl_used = [11000] ;zsfc_max = 2000.0 coast_mask_max = 0 end "dbt11max3x3": begin classifier_name = 'dbt11max3x3' xtitle = 'dbt11max3x3' x = f.temp_11_0um_nom_max_3x3 - f.temp_11_0um_nom x_min = -1.0 & x_bin = 0.1 & nbins_x = 100.0 nchan_used = 1 & channel_wvl_used = [11000] coast_mask_max = 0 zsfc_std_max = zsfc_std_max_iruni end "logdbt11max3x3": begin classifier_name = 'logdbt11max3x3' xtitle = 'logdbt11max3x3' x = f.temp_11_0um_nom_max_3x3 - f.temp_11_0um_nom x_min = -0.12 & x_bin = 0.02 & nbins_x = 100.0 idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx]) idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min nchan_used = 1 & channel_wvl_used = [11000] coast_mask_max = 0 zsfc_std_max = zsfc_std_max_iruni end "bt11std": begin classifier_name = 'bt11std' xtitle = 'bt11std' x = f.temp_11_0um_nom_stddev_3x3 x_min = 0.0 & x_bin = 0.1 & nbins_x = 100.0 nchan_used = 1 & channel_wvl_used = [11000] coast_mask_max = 0 ;zsfc_max = 2000.0 zsfc_std_max = zsfc_std_max_iruni end "logbt11std": begin classifier_name = 'logbt11std' xtitle = 'logbt11std' x = f.temp_11_0um_nom_stddev_3x3 x_min = -1.5 & x_bin = 0.1 & nbins_x = 36.0 idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx]) idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min nchan_used = 1 & channel_wvl_used = [11000] coast_mask_max = 0 ;zsfc_max = 2000.0 zsfc_std_max = zsfc_std_max_iruni end "bt10std": begin classifier_name = 'bt10std' xtitle = 'bt10std' x = f.temp_10_4um_nom_stddev_3x3 x_min = 0.0 & x_bin = 0.1 & nbins_x = 100.0 nchan_used = 1 & channel_wvl_used = [10400] coast_mask_max = 0 ;zsfc_max = 2000.0 zsfc_std_max = zsfc_std_max_iruni end "logbt10std": begin classifier_name = 'logbt10std' xtitle = 'logbt10std' x = f.temp_10_4um_nom_stddev_3x3 x_min = -1.5 & x_bin = 0.1 & nbins_x = 36.0 idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx]) idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min nchan_used = 1 & channel_wvl_used = [10400] coast_mask_max = 0 ;zsfc_max = 2000.0 zsfc_std_max = zsfc_std_max_iruni end "dbt10max3x3": begin classifier_name = 'dbt10max3x3' xtitle = 'dbt10max3x3' x = f.temp_10_4um_nom_max_3x3 - f.temp_10_4um_nom x_min = 0.0 & x_bin = 0.2 & nbins_x = 100.0 nchan_used = 1 & channel_wvl_used = [10400] coast_mask_max = 0 ;zsfc_max = 2000.0 zsfc_std_max = zsfc_std_max_iruni end "btd3811_all": begin classifier_name = 'btd3811_all' xtitle = 'btd3811' x = f.temp_3_75um_nom - f.temp_11_0um_nom x_min = -5 & x_bin = 0.5 & nbins_x = 60.0 nchan_used = 2 & channel_wvl_used = [3750,11000] solzen_min = 0.0 solzen_max = 180.0 solglintzen_min = 0.0 solglintzen_max = 180.0 tsfc_min = 240.0 end "btd3811_day": begin classifier_name = 'btd3811_day' xtitle = 'btd3811' x = f.temp_3_75um_nom - f.temp_11_0um_nom x_min = -5 & x_bin = 0.5 & nbins_x = 60.0 nchan_used = 2 & channel_wvl_used = [3750,11000] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = 95.0 solglint_mask_min = 0 solglint_mask_max = 0 tsfc_min = 240.0 end "btd3811_night": begin classifier_name = 'btd3811_night' xtitle = 'btd3811' x = f.temp_3_75um_nom - f.temp_11_0um_nom x_min = -10 & x_bin = 0.5 & nbins_x = 60.0 nchan_used = 2 & channel_wvl_used = [3750,11000] solzen_min = 90.0 solzen_max = 180.0 tsfc_min = 240.0 end "dbtd3811clr_day": begin classifier_name = 'dbtd3811clr_day' xtitle = 'dbtd3811clr' x = (f.temp_3_75um_nom - f.temp_11_0um_nom) - $ (f.temp_3_75um_nom_clear_sky - f.temp_11_0um_nom_clear_sky) x_min = -10 & x_bin = 0.5 & nbins_x = 80.0 nchan_used = 2 & channel_wvl_used = [3750,11000] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = 95.0 solglint_mask_min = 0 solglint_mask_max = 0 tsfc_min = 240.0 end "dbtd3811clr_night": begin classifier_name = 'dbtd3811clr_night' xtitle = 'dbtd3811clr' x = (f.temp_3_75um_nom - f.temp_11_0um_nom) - $ (f.temp_3_75um_nom_clear_sky - f.temp_11_0um_nom_clear_sky) x_min = -10 & x_bin = 0.5 & nbins_x = 80.0 nchan_used = 2 & channel_wvl_used = [3750,11000] solzen_min = 90.0 solzen_max = 180.0 tsfc_min = 240.0 end "dbtd3811clr_all": begin classifier_name = 'dbtd3811clr_all' xtitle = 'dbtd3811clr' x = (f.temp_3_75um_nom - f.temp_11_0um_nom) - $ (f.temp_3_75um_nom_clear_sky - f.temp_11_0um_nom_clear_sky) x_min = -10 & x_bin = 0.5 & nbins_x = 80.0 nchan_used = 2 & channel_wvl_used = [3750,11000] tsfc_min = 240.0 end "btd3810_all": begin classifier_name = 'btd3810_all' xtitle = 'btd3810' x = f.temp_3_75um_nom - f.temp_10_4um_nom x_min = -5 & x_bin = 0.5 & nbins_x = 60.0 nchan_used = 2 & channel_wvl_used = [3750,11000] solzen_min = 0.0 solzen_max = 180.0 solglintzen_min = 0.0 solglintzen_max = 180.0 tsfc_min = 240.0 end "btd3810_day": begin classifier_name = 'btd3810_day' xtitle = 'btd3810' x = f.temp_3_75um_nom - f.temp_10_4um_nom x_min = -5 & x_bin = 0.5 & nbins_x = 60.0 nchan_used = 2 & channel_wvl_used = [3750,10400] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = 95.0 solglint_mask_min = 0 solglint_mask_max = 0 tsfc_min = 240.0 end "btd3810_night": begin classifier_name = 'btd3810_night' xtitle = 'btd3810' x = f.temp_3_75um_nom - f.temp_10_4um_nom x_min = -5 & x_bin = 0.5 & nbins_x = 60.0 nchan_used = 2 & channel_wvl_used = [3750,10400] solzen_min = 90.0 solzen_max = 180.0 tsfc_min = 240.0 end "btd8511": begin classifier_name = 'btd8511' xtitle = 'btd8511' x = f.temp_8_5um_nom - f.temp_11_0um_nom x_min = -12.0 & x_bin = 0.2 & nbins_x = 100 nchan_used = 2 & channel_wvl_used = [8500,11000] tsfc_min = 230.0 end "btd1112": begin classifier_name = 'btd1112' xtitle = 'btd1112' x = f.temp_11_0um_nom - f.temp_12_0um_nom x_min = -2.0 & x_bin = 0.1 & nbins_x = 100 nchan_used = 2 & channel_wvl_used = [11000,12000] tsfc_min = 240.0 end "fmft": begin classifier_name = 'fmft' xtitle = 'fmft' x = fmft(f.temp_11_0um_nom,f.temp_12_0um_nom, $ f.temp_11_0um_nom_clear_sky,f.temp_12_0um_nom_clear_sky,missing) x_min = -2.0 & x_bin = 0.1 & nbins_x = 100 nchan_used = 2 & channel_wvl_used = [11000,12000] tsfc_min = 230.0 end "fmft_day": begin classifier_name = 'fmft_day' xtitle = 'fmft' x = fmft(f.temp_11_0um_nom,f.temp_12_0um_nom, $ f.temp_11_0um_nom_clear_sky,f.temp_12_0um_nom_clear_sky,missing) x_min = -2.0 & x_bin = 0.1 & nbins_x = 100 nchan_used = 2 & channel_wvl_used = [11000,12000] tsfc_min = 230.0 solzen_max = 95.0 end "fmft_night": begin classifier_name = 'fmft_night' xtitle = 'fmft' x = fmft(f.temp_11_0um_nom,f.temp_12_0um_nom, $ f.temp_11_0um_nom_clear_sky,f.temp_12_0um_nom_clear_sky,missing) x_min = -2.0 & x_bin = 0.1 & nbins_x = 100 nchan_used = 2 & channel_wvl_used = [11000,12000] tsfc_min = 230.0 solzen_min = 85.0 end "drefl065clr": begin classifier_name = 'drefl065clr' xtitle = 'drefl065clr' x = f.refl_0_65um_nom - f.refl_0_65um_nom_clear_sky x_min = -40.0 & x_bin = 1.0 & nbins_x = 120 nchan_used = 1 & channel_wvl_used = [650] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min solglint_mask_min = 0 solglint_mask_max = 0 snow_class_max = 1 end "drefl160clr": begin classifier_name = 'drefl160clr' xtitle = 'drefl160clr' x = f.refl_1_60um_nom - f.refl_1_60um_nom_clear_sky x_min = -40.0 & x_bin = 1.0 & nbins_x = 120 nchan_used = 1 & channel_wvl_used = [650] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min solglint_mask_min = 0 solglint_mask_max = 0 end "drefl065maxsubclr": begin classifier_name = 'drefl065maxsubclr' xtitle = 'drefl065maxsubclr' x = f.refl_0_65um_nom_max_sub - f.refl_0_65um_nom_clear_sky x_min = -40.0 & x_bin = 1.0 & nbins_x = 120 nchan_used = 1 & channel_wvl_used = [650] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min solglint_mask_min = 0 solglint_mask_max = 0 coast_mask_max = 0 end "drefl065min3x3": begin classifier_name = 'drefl065min3x3' xtitle = 'drefl065min3x3' x = f.refl_0_65um_nom - f.refl_0_65um_nom_min_3x3 x_min = 0.0 & x_bin = 0.5 & nbins_x = 100 nchan_used = 1 & channel_wvl_used = [650] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min coast_mask_max = 0 end "logdrefl065minsub": begin classifier_name = 'logdrefl065minsub' xtitle = 'logdrefl065minsub' x = f.refl_0_65um_nom - f.refl_0_65um_nom_min_sub x_min = -2 & x_bin = 0.04 & nbins_x = 100 idx = where(x gt 0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx]) idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min nchan_used = 1 & channel_wvl_used = [650] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min coast_mask_max = 0 ;zsfc_max = 2000.0 end "drefl065minsub": begin classifier_name = 'drefl065minsub' xtitle = 'drefl065minsub' x = f.refl_0_65um_nom - f.refl_0_65um_nom_min_sub x_min = 0.0 & x_bin = 0.5 & nbins_x = 100 nchan_used = 1 & channel_wvl_used = [650] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min coast_mask_max = 0 ;zsfc_max = 2000.0 end ;"refl065maxsub": begin ; classifier_name = 'refl065maxsub' ; xtitle = 'refl065maxsub' ; x = f.refl_0_65um_nom_max_sub ; x_min = 0.0 & x_bin = 0.5 & nbins_x = 100 ; nchan_used = 1 & channel_wvl_used = [650] ; solzen_min = 0.0 ; solzen_max = solzen_day_max ; solglint_mask_min = 0 ; solglint_mask_max = 0 ; coast_mask_max = 0 ; ;zsfc_max = 2000.0 ;end "refl065maxsub": begin classifier_name = 'refl065maxsub' xtitle = 'refl065maxsub' x = f.refl_0_65um_nom x_min = 0.0 & x_bin = 0.5 & nbins_x = 100 nchan_used = 1 & channel_wvl_used = [650] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min solglint_mask_min = 0 solglint_mask_max = 0 coast_mask_max = 0 ;zsfc_max = 2000.0 end "drefl065maxminsub": begin classifier_name = 'drefl065maxminsub' xtitle = 'drefl065maxminsub' x = f.refl_0_65um_nom_max_sub - f.refl_0_65um_nom_min_sub x_min = 0.0 & x_bin = 0.5 & nbins_x = 100 nchan_used = 1 & channel_wvl_used = [650] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min ;zsfc_max = 2000.0 coast_mask_max = 0 end "logdrefl065maxminsub": begin classifier_name = 'logdrefl065maxminsub' xtitle = 'logdrefl065maxminsub' x = f.refl_0_65um_nom_max_sub - f.refl_0_65um_nom_min_sub x_min = -2 & x_bin = 0.04 & nbins_x = 100 idx = where(x gt 0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx]) idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min nchan_used = 1 & channel_wvl_used = [650] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min ;zsfc_max = 2000.0 coast_mask_max = 0 end "refrat086065": begin classifier_name = 'refrat086065' xtitle = 'refrat086065' x = f.refl_0_86um_nom / f.refl_0_65um_nom x_min = 0.0 & x_bin = 0.02 & nbins_x = 100 nchan_used = 2 & channel_wvl_used = [650,860] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min solglint_mask_min = 0 solglint_mask_max = 0 snow_class_max = 1 end "refrat138065": begin classifier_name = 'refrat138065' xtitle = 'refrat138065' x = f.refl_1_38um_nom / f.refl_0_65um_nom x_min = 0.0 & x_bin = 0.02 & nbins_x = 50 nchan_used = 2 & channel_wvl_used = [065,1380] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min slant_tpw_min = 0.5 ;zsfc_max = 2000.0 snow_class_max = 1 end "logcoddnb": begin classifier_name = 'logcoddnb' xtitle = 'logcoddnb' x = f.cld_opd_mask_dnb_nom x_min = -2 & x_bin = 0.1 & nbins_x = 45.0 idx = where(x gt 0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx]) idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min nchan_used = 1 & channel_wvl_used = [700] solzen_min = 90.0 solzen_max = 180.0 lunzen_min = 0.0 lunzen_max = 80.0 lunglint_mask_min = 0 lunglint_mask_max = 0 snow_class_max = 1 end "logcod065": begin classifier_name = 'logcod065' xtitle = 'logcod065' x = f.cld_opd_mask_0_65um_nom x_min = -2 & x_bin = 0.1 & nbins_x = 45.0 idx = where(x gt 0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx]) idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min nchan_used = 1 & channel_wvl_used = [650] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min solglint_mask_min = 0 solglint_mask_max = 0 snow_class_max = 1 end "logcod138": begin classifier_name = 'logcod138' xtitle = 'logcod138' x = f.cld_opd_mask_1_38um_nom x_min = -2 & x_bin = 0.1 & nbins_x = 45.0 idx = where(x gt 0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx]) idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min nchan_used = 1 & channel_wvl_used = [1380] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min slant_tpw_min = 0.5 ;zsfc_max = 2000.0 end "emiss3811_day": begin classifier_name = 'emiss3811_day' xtitle = 'emiss3811' x = f.emiss_3_75um_nom x_min = 0.25 & x_bin = 0.05 & nbins_x = 70.0 nchan_used = 2 & channel_wvl_used = [3750,11000] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min solglint_mask_min = 0 solglint_mask_max = 0 tsfc_min = 240.0 tsfc_min = 230.0 end "emiss3811_night": begin classifier_name = 'emiss3811_night' xtitle = 'emiss3811' x = f.emiss_3_75um_nom x_min = 0.25 & x_bin = 0.05 & nbins_x = 70.0 ; x = f.emiss_3_75um_nom - f.emiss_3_75um_nom_clear ; x_min = -0.5 & x_bin = 0.05 & nbins_x = 70.0 nchan_used = 2 & channel_wvl_used = [3750,11000] solzen_min = 90.0 solzen_max = 180.0 tsfc_min = 240.0 tsfc_min = 230.0 end "emiss3810_day": begin classifier_name = 'emiss3810_day' xtitle = 'emiss3810' x = f.emiss_3_75um_nom_rel_10_4um_nom x_min = 0.25 & x_bin = 0.05 & nbins_x = 70.0 nchan_used = 2 & channel_wvl_used = [3750,10400] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min solglint_mask_min = 0 solglint_mask_max = 0 tsfc_min = 240.0 tsfc_min = 230.0 end "emiss3810_night": begin classifier_name = 'emiss3810_night' xtitle = 'emiss3810' x = f.emiss_3_75um_nom_rel_10_4um_nom x_min = 0.25 & x_bin = 0.05 & nbins_x = 70.0 nchan_used = 2 & channel_wvl_used = [3750,10400] solzen_min = 90.0 solzen_max = 180.0 tsfc_min = 240.0 tsfc_min = 230.0 end "refldnb": begin classifier_name = 'refldnb' xtitle = 'refldnb' x = f.refl_lunar_dnb_nom x_min = -5 & x_bin = 1.0 & nbins_x = 80 nchan_used = 1 & channel_wvl_used = [650] solzen_min = 90.0 solzen_max = 180.0 lunzen_min = 0.0 lunzen_max = 80.0 lunglint_mask_min = 0 lunglint_mask_max = 0 moon_illum_frac_min = 25.0 moon_illum_frac_max = 100.0 snow_class_min = 1 snow_class_max = 1 end "ndsi_dnb_38": begin classifier_name = 'ndsi_dnb_38' xtitle = 'ndsi_dnb_38' x = ndsi(abs(f.refl_lunar_dnb_nom), abs(f.refl_3_75um_nom), missing) x_min = -0.5 & x_bin = 0.025 & nbins_x = 60 nchan_used = 2 & channel_wvl_used = [650,3750] solzen_min = 90.0 solzen_max = 180.0 lunzen_min = 0.0 lunzen_max = 80.0 lunglint_mask_min = 0 lunglint_mask_max = 0 moon_illum_frac_min = 25.0 moon_illum_frac_max = 100.0 snow_class_min = 1 snow_class_max = 3 end "ndsi": begin classifier_name = 'ndsi' xtitle = 'ndsi' x = ndsi(f.refl_0_65um_nom, f.refl_1_60um_nom, missing) x_min = -0.5 & x_bin = 0.025 & nbins_x = 60 nchan_used = 2 & channel_wvl_used = [650,860] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min solglint_mask_min = 0 solglint_mask_max = 0 coast_mask_max = 0.0 end "refl047": begin classifier_name = 'refl047' xtitle = 'refl047' x = f.refl_0_47um_nom x_min = -5 & x_bin = 2.0 & nbins_x = 60 nchan_used = 1 & channel_wvl_used = [470] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min solglint_mask_min = 0 solglint_mask_max = 0 snow_class_max = 1 end "refl065": begin classifier_name = 'refl065' xtitle = 'refl065' x = f.refl_0_65um_nom x_min = -5 & x_bin = 1.0 & nbins_x = 80 nchan_used = 1 & channel_wvl_used = [650] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min solglint_mask_min = 0 solglint_mask_max = 0 snow_class_min = 1 snow_class_max = 1 end "refl160": begin classifier_name = 'refl160' xtitle = 'refl160' x = f.refl_1_60um_nom x_min = -5 & x_bin = 1.0 & nbins_x = 80 nchan_used = 1 & channel_wvl_used = [1600] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min solglint_mask_min = 0 solglint_mask_max = 0 end "refl138": begin classifier_name = 'refl138' xtitle = 'refl138' x = f.refl_1_38um_nom x_min = -5 & x_bin = 0.25 & nbins_x = 80 nchan_used = 1 & channel_wvl_used = [1380] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min slant_tpw_min = 0.5 ;zsfc_max = 2000.0 end "refl38": begin classifier_name = 'refl38' xtitle = 'refl38' x = f.refl_3_75um_nom x_min = -5 & x_bin = 0.4 & nbins_x = 80 nchan_used = 1 & channel_wvl_used = [3750] solzen_min = 0.0 solzen_max = solzen_day_max solscatang_min = solscatang_day_min tsfc_min = 240.0 end "bt6711covar": begin classifier_name = 'bt6711covar' xtitle = 'bt6711covar' x = f.temp_11um_vs_67um_covar_5x5 x_min = -4 & x_bin = 0.15 & nbins_x = 80 nchan_used = 2 & channel_wvl_used = [6700,11000] slant_tpw_min = 0.5 ;zsfc_max = 2000.0 end "bt6211covar": begin classifier_name = 'bt6211covar' xtitle = 'bt6211covar' x = f.temp_11um_vs_62um_covar_5x5 x_min = -4 & x_bin = 0.15 & nbins_x = 80 nchan_used = 2 & channel_wvl_used = [6200,11000] slant_tpw_min = 0.5 ;zsfc_max = 2000.0 end "btd7362": begin classifier_name = 'btd7362' xtitle = 'btd7362' x = f.temp_7_3um_nom - f.temp_6_2um_nom x_min = -10.0 & x_bin = 0.5 & nbins_x = 100 nchan_used = 2 & channel_wvl_used = [6200,7300] ;zsfc_max = 2000.0 end "btd1167": begin classifier_name = 'btd1167' xtitle = 'btd1167' x = f.temp_11_0um_nom - f.temp_6_7um_nom x_min = -5.0 & x_bin = 2.0 & nbins_x = 40 nchan_used = 2 & channel_wvl_used = [6700,11000] ;zsfc_max = 2000.0 end "btd1162": begin classifier_name = 'btd1162' xtitle = 'btd1162' x = f.temp_11_0um_nom - f.temp_6_2um_nom x_min = -5.0 & x_bin = 2.0 & nbins_x = 40 nchan_used = 2 & channel_wvl_used = [6200,11000] ;zsfc_max = 2000.0 end "btd1167_day": begin classifier_name = 'btd1167_day' xtitle = 'btd1167' x = f.temp_11_0um_nom - f.temp_6_7um_nom x_min = -5.0 & x_bin = 2.0 & nbins_x = 40 nchan_used = 2 & channel_wvl_used = [6700,11000] ;zsfc_max = 2000.0 solzen_max = 95.0 end "btd1167_night": begin classifier_name = 'btd1167_night' xtitle = 'btd1167' x = f.temp_11_0um_nom - f.temp_6_7um_nom x_min = -10.0 & x_bin = 2.0 & nbins_x = 45 nchan_used = 2 & channel_wvl_used = [6700,11000] ;zsfc_max = 2000.0 solzen_min = 85.0 end "btd8573": begin classifier_name = 'btd8573' xtitle = 'btd8573' x = f.temp_8_5um_nom - f.temp_7_3um_nom x_min = -5.0 & x_bin = 1.0 & nbins_x = 50 nchan_used = 2 & channel_wvl_used = [7300,11000] ;zsfc_max = 2000.0 end "btd13373": begin classifier_name = 'btd13373' xtitle = 'btd13373' x = f.temp_13_3um_nom - f.temp_7_3um_nom x_min = -10.0 & x_bin = 1.0 & nbins_x = 30 nchan_used = 2 & channel_wvl_used = [7300,13300] ;zsfc_max = 2000.0 end "btd11133": begin classifier_name = 'btd11133' xtitle = 'btd11133' x = f.temp_11_0um_nom - f.temp_13_3um_nom x_min = -10.0 & x_bin = 2.0 & nbins_x = 40 nchan_used = 2 & channel_wvl_used = [1100,13300] ; ;zsfc_max = 2000.0 end "btd1173": begin classifier_name = 'btd1173' xtitle = 'btd1173' x = f.temp_11_0um_nom - f.temp_7_3um_nom x_min = -10.0 & x_bin = 2.0 & nbins_x = 45 nchan_used = 2 & channel_wvl_used = [7300,11000] ;zsfc_max = 2000.0 end "btd1173_day": begin classifier_name = 'btd1173_day' xtitle = 'btd1173' x = f.temp_11_0um_nom - f.temp_7_3um_nom x_min = -10.0 & x_bin = 2.0 & nbins_x = 45 nchan_used = 2 & channel_wvl_used = [7300,11000] ;zsfc_max = 2000.0 solzen_max = 95.0 end "btd1173_night": begin classifier_name = 'btd1173_night' xtitle = 'btd1173' x = f.temp_11_0um_nom - f.temp_7_3um_nom x_min = -5.0 & x_bin = 2.0 & nbins_x = 40 nchan_used = 2 & channel_wvl_used = [7300,11000] ;zsfc_max = 2000.0 solzen_min = 85.0 end else: begin print, "error: unknown test" goto, skip_all end endcase ;----- ;end of class ;----- ;--- but valid values outside border at lowest bin ;idx = where(x lt x_min and x ne missing,cc) & if (cc gt 0) then x[idx] = x_min print, "+++" print, "GENERATING LUT FOR ", classifier_name print, "+++" ;------------------------------------------------ ; fake other dims ;------------------------------------------------ ytitle = "blank" nbins_y = 1 y_min = 0.0 y_bin = 0.0 ztitle = "blank" nbins_z = 1 z_min = 0.0 z_bin = 0.0 ;------------------------------------------------ ; make bin bounds ;------------------------------------------------ x_x = x_min + findgen(nbins_x)*x_bin + x_bin/2 ;----------------------------------------------- ; make surface type ;---------------------------------------------- sfc_type_mask = f.bayes_mask_sfc_type if (keyword_set(isccp)) then begin make_ecm_surface_type, f, sfc_type_mask_2, /isccp endif ;----------------------------------------------- ; apply filters ;---------------------------------------------- cod_lidar = f.CLOSEST_CALIPSO_COD slant_tpw = f.total_precipitable_water_nwp / cos(f.sensor_zenith_angle * !DTOR) if (keyword_set(new)) then begin phase_lidar = f.CLOSEST_CALIPSO_CLOUD_PHASE_TOP_LAYER frac_lidar = f.CLOSEST_CALIPSO_CLOUD_FRACTION_5KM tc_lidar = f.CLOSEST_CALIPSO_TOP_MID_TEMPERATURE frac_lidar = f.CLOSEST_CALIPSO_CLOUD_FRACTION_5km codi_lidar = f.CLOSEST_CALIPSO_COD_ICE endif if (keyword_set(old) or keyword_set(denis)) then begin frac_lidar =f.CLOSEST_CALIPSO_CLOUD_FRACTION_x5 codi_lidar = f.CLOSEST_CALIPSO_COD codi_lidar[*] = missing if (~keyword_set(denis)) then phase_lidar = f.CLOSEST_CALIPSO_CLOUD_PHASE_TOP_LAYER if (keyword_set(denis)) then begin phase_lidar = make_array(frac_lidar.length,/INT,VALUE=-1) idx = where(frac_lidar lt 0.5 and frac_lidar gt 0.0,cc) if (cc gt 0) then phase_lidar[idx] = 0 idx = where(frac_lidar ge 0.5,cc) & if (cc gt 0) then phase_lidar[idx] = 2 idx = where(phase_lidar eq 2 and $ f.cld_temp_opaque ge 0.0 and f.cld_temp_opaque lt 260.0,cc) if (cc gt 0) then phase_lidar[idx] = 1 endif endif ;--- 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 make_true_phase, phase_lidar, frac_lidar, tc_lidar, cod_lidar, codi_lidar, $ cod_clear_water_thresh, cod_clear_ice_thresh, $ f.bayes_mask_sfc_type, f.cld_temp_opaque, f.cld_height_opaque, $ f.solar_zenith_angle, f.glint_mask, f.coast_mask, $ f.surface_elevation, f.surface_temperature_nwp, $ f.refl_0_65um_nom_min_3x3, f.refl_0_65um_nom, $ f.temp_11_0um_nom_max_3x3, f.temp_11_0um_nom, $ f.cld_opd_mask_0_65um_nom, cod_cloud_thresh, $ f.emiss_tropo_11_0um_nom, etrop_cloud_thresh, mod_flag, phase_true, missing ;--- begin do no filtering ;phase_true = phase_lidar ;idx = where(phase_lidar eq 2, cc) ;if (cc gt 0) then phase_true[idx] = 1 ;idx = where(phase_lidar eq 1 or phase_lidar eq 3, cc) ;if (cc gt 0) then phase_true[idx] = 2 ;idx = where(phase_lidar lt 0 or phase_lidar gt 3, cc) ;if (cc gt 0) then phase_true[idx] = -1 ; ;idx = where(phase_true eq -1, cc_bad) ;idx = where(phase_true eq 0, cc0) ;idx = where(phase_true eq 1, cc1) ;idx = where(phase_true eq 2, cc2) ;n = float(cc0+cc1+cc2) ; ;print, "initial phase fractions(bad/c/i/w) = ", float(cc_bad)/(n+cc_bad),float(cc0)/n, float(cc1)/n, float(cc2)/n ;--- end do no filtering ;------------------------------------------------------------------------------- ; fill in things that might be missing in training data ;------------------------------------------------------------------------------- idx = where(tag_names(f) eq "CITY_MASK",cc) if (cc eq 0) then begin city_mask_temp = make_array(ntrain,/FLOAT, VALUE=MISSING) endif else begin city_mask_temp = f.city_mask endelse idx = where(tag_names(f) eq "MOON_ILLUM_FRAC",cc) if (cc eq 0) then begin moon_illum_frac_temp = make_array(ntrain,/FLOAT, VALUE=MISSING) endif else begin moon_illum_frac_temp = f.moon_illum_frac endelse ;------------------------------------------------------------------------------- ; apply filters ;------------------------------------------------------------------------------- ;print, 'before apply filters ', n_elements(city_mask_temp), city_mask_min, city_mask_max apply_filters, phase_true, $ f.CLOSEST_CALIPSO_TIME_DIFFERENCE, Time_Diff_Max, $ f.sensor_zenith_angle, zen_min, zen_max, $ f.solar_zenith_angle, solzen_min, solzen_max, $ f.lunar_zenith_angle, lunzen_min, lunzen_max, $ f.land_class, $ frac_lidar, cfrac_min, cfrac_max, $ f.glint_zenith_angle, solglintzen_min, solglintzen_max, $ f.lunar_glint_zenith_angle, lunglintzen_min, lunglintzen_max, $ f.scattering_angle, solscatang_min, solscatang_max, $ f.snow_class, snow_class_min, snow_class_max, $ f.glint_mask, solglint_mask_min, solglint_mask_max, $ f.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, $ f.coast_mask, coast_mask_min, coast_mask_max, $ slant_tpw, slant_tpw_min, slant_tpw_max, $ f.surface_temperature_nwp, tsfc_min, tsfc_max, $ f.surface_elevation, zsfc_min, zsfc_max, $ f.surface_elevation_stddev_3x3, zsfc_std_min, zsfc_std_max, $ ;f.latitude, lat_min, lat_max, missing,/phase_uni_flag f.latitude, lat_min, lat_max, missing ;------------------------------- ; select valid data within these limits ;------------------------------- idx = where(phase_true ge 0 and x ne missing, cc) ;--- through away filtered data x = x[idx] phase_true = phase_true[idx] sfc_type_mask = sfc_type_mask[idx] ;------------------------------------------------------------- ; make arrays ;------------------------------------------------------------- ice_table = make_array(nbins_x, nsfc, /FLOAT, VALUE=MISSING) water_table = make_array(nbins_x, nsfc, /FLOAT, VALUE=MISSING) clear_table = make_array(nbins_x, nsfc, /FLOAT, VALUE=MISSING) obs_table = make_array(nbins_x, nsfc, /FLOAT, VALUE=0.0) cloud_fraction = make_array(nsfc,/FLOAT,VALUE=missing) clear_fraction = make_array(nsfc,/FLOAT,VALUE=missing) ice_fraction = make_array(nsfc,/FLOAT,VALUE=missing) water_fraction = make_array(nsfc,/FLOAT,VALUE=missing) observation_count = make_array(nsfc,/FLOAT,VALUE=missing) ;-------------------------------------------------------------------- ; loop over surface type ;-------------------------------------------------------------------- count_min = 100 for isfc = 0, nsfc - 1 do begin ;sfc_loop idx = where(sfc_type_mask eq isfc + 1 and phase_true ge 0,cc) ;print, 'isfc, count = ', isfc, cc if (cc gt count_min) then begin x_sfc = x[idx] phase_true_sfc = phase_true[idx] observation_count[isfc] = cc ; print, 'isfc, count = ', isfc, cc endif else begin goto, cycle_sfc endelse ;--- make priors idx = where(phase_true_sfc ge 0, count_all) idx = where(phase_true_sfc eq 2, count_yes_ice_all) idx = where(phase_true_sfc eq 1, count_yes_water_all) idx = where(phase_true_sfc eq 0, count_yes_clear_all) cloud_fraction[isfc] = float(count_yes_water_all + count_yes_ice_all) / float(count_all) ice_fraction[isfc] = float(count_yes_ice_all) / float(count_all) water_fraction[isfc] = float(count_yes_water_all) / float(count_all) clear_fraction[isfc] = 1.0 - cloud_fraction[isfc] cond_ratio_ice = make_array(nbins_x, /FLOAT, VALUE=MISSING) cond_ratio_water = make_array(nbins_x, /FLOAT, VALUE=MISSING) cond_ratio_clear = make_array(nbins_x, /FLOAT, VALUE=MISSING) post_prob_ice = make_array(nbins_x, /FLOAT, VALUE=MISSING) post_prob_water = make_array(nbins_x, /FLOAT, VALUE=MISSING) post_prob_clear = make_array(nbins_x, /FLOAT, VALUE=MISSING) obs_prob_lidar = make_array(nbins_x, /FLOAT, VALUE=0.0) for i = 0, nbins_x -1 do begin ;point_loop x_min_bin = x_min + i*x_bin x_max_bin = x_min + (i+1)*x_bin idx_all_bin = where(x_sfc ge x_min_bin and x_sfc lt x_max_bin and $ phase_true_sfc ge 0, count_all_bin) if (count_all_bin eq 0) then goto, skip_bin if (count_all_bin lt count_min) then goto, skip_bin idx_ice = where(phase_true_sfc[idx_all_bin] eq 2, count_yes_ice_bin) idx_water = where(phase_true_sfc[idx_all_bin] eq 1, count_yes_water_bin) idx_clear = where(phase_true_sfc[idx_all_bin] eq 0, count_yes_clear_bin) ;--- clear count_yes_bin = count_yes_clear_bin count_yes_all = count_yes_clear_all count_no_bin = count_all_bin - count_yes_bin count_no_all = count_all - count_yes_all prior_yes = clear_fraction[isfc] compute_class_cond, nb_rank, laplace_flag, $ count_yes_all, count_no_all, count_yes_bin, count_no_bin, $ max_r, prior_yes, cond_yes, cond_no, r, posterior_prob cond_ratio_clear[i] = r post_prob_clear[i] = posterior_prob ;--- water count_yes_bin = count_yes_water_bin count_yes_all = count_yes_water_all count_no_bin = count_all_bin - count_yes_bin count_no_all = count_all - count_yes_all prior_yes = water_fraction[isfc] compute_class_cond, nb_rank, laplace_flag, $ count_yes_all, count_no_all, count_yes_bin, count_no_bin, $ max_r, prior_yes, cond_yes, cond_no, r, posterior_prob cond_ratio_water[i] = r post_prob_water[i] = posterior_prob ;--- ice count_yes_bin = count_yes_ice_bin count_yes_all = count_yes_ice_all count_no_bin = count_all_bin - count_yes_bin count_no_all = count_all - count_yes_all prior_yes = ice_fraction[isfc] compute_class_cond, nb_rank, laplace_flag, $ count_yes_all, count_no_all, count_yes_bin, count_no_bin, $ max_r, prior_yes, cond_yes, cond_no, r, posterior_prob cond_ratio_ice[i] = r post_prob_ice[i] = posterior_prob ;---- test ; print, "counts ", count_all_bin, count_yes_clear_bin, count_yes_water_bin, count_yes_ice_bin, $ ; count_yes_clear_bin + count_yes_water_bin + count_yes_ice_bin ; print, "probs ", post_prob_clear[i], post_prob_water[i], post_prob_ice[i] ; print, "prob sum ", post_prob_clear[i] + post_prob_water[i] + post_prob_ice[i] ;---- end test if (count_all gt count_min) then begin obs_prob_lidar[i] = float(count_all_bin) endif skip_bin: endfor ;end of pixel loop ;--- normalize obs_prob idx = where(obs_prob_lidar gt 0,cc) if (cc gt 0) then obs_prob_lidar[idx] = obs_prob_lidar[idx] / total(obs_prob_lidar) ;------------------- ; fill in holes ;------------------- ;goto, skip_fill ;idx = where(obs_prob_lidar eq 0,cc) & print, 'number with zero = ', cc obs_prob_interp_value = 1.0e-06 for ix = 0, nbins_x - 1 do begin if (obs_prob_lidar[ix] eq 0) then begin ix_start = max([0,ix-1]) ix_end = min([nbins_x-1,ix+1]) v = obs_prob_lidar[ix_start:ix_end] idx = where(v gt 0,cc) if (cc gt 0) then begin obs_prob_lidar[ix] = mean(v[idx]) v = cond_ratio_clear[ix_start:ix_end] cond_ratio_clear[ix] = mean(v[idx]) v = cond_ratio_water[ix_start:ix_end] cond_ratio_water[ix] = mean(v[idx]) v = cond_ratio_ice[ix_start:ix_end] cond_ratio_ice[ix] = mean(v[idx]) v = post_prob_clear[ix_start:ix_end] post_prob_clear[ix] = mean(v[idx]) v = post_prob_water[ix_start:ix_end] post_prob_water[ix] = mean(v[idx]) v = post_prob_ice[ix_start:ix_end] post_prob_ice[ix] = mean(v[idx]) endif endif endfor ;--- renormalize obs_prob idx = where(obs_prob_lidar gt 0,cc) if (cc gt 0) then obs_prob_lidar[idx] = obs_prob_lidar[idx] / total(obs_prob_lidar) skip_fill: ;----------------------------------------------------------------------------------------- ; smooth ;----------------------------------------------------------------------------------------- ;goto, skip_smooth smooth_width = 3 cond_ratio_clear = smooth_nan(cond_ratio_clear,smooth_width,missing) cond_ratio_water = smooth_nan(cond_ratio_water,smooth_width,missing) cond_ratio_ice = smooth_nan(cond_ratio_ice,smooth_width,missing) post_prob_clear = smooth_nan(post_prob_clear,smooth_width,missing) post_prob_water = smooth_nan(post_prob_water,smooth_width,missing) post_prob_ice = smooth_nan(post_prob_ice,smooth_width,missing) skip_smooth: ;--- normalize obs_prob idx = where(obs_prob_lidar gt 0,cc) if (cc gt 0) then obs_prob_lidar[idx] = obs_prob_lidar[idx] / total(obs_prob_lidar) ;------------------------------------------------------------------------------------ ; apply table ;------------------------------------------------------------------------------------ n = n_elements(phase_true_sfc) phase_retrieved = make_array(n, /FLOAT, VALUE = missing) for k = 0, n - 1 do begin i = min([nbins_x -1,max([0,fix((x_sfc[k] - x_min) / x_bin)])]) ; max_prob = max([post_prob_clear[i], post_prob_water[i] , post_prob_ice[i]]) ; if (post_prob_clear[i] eq max_prob) then phase_retrieved[k] = 0 ; if (post_prob_water[i] eq max_prob) then phase_retrieved[k] = 1 ; if (post_prob_ice[i] eq max_prob) then phase_retrieved[k] = 2 phase_retrieved[k] = 0 if (post_prob_clear[i] lt 0.50) then begin phase_retrieved[k] = 1 if (post_prob_ice[i] gt post_prob_water[i]) then begin phase_retrieved[k] = 2 endif endif endfor idx = where(phase_true_sfc ge 0 and phase_retrieved ge 0, count_all) idx = where(phase_true_sfc eq 0 and phase_retrieved eq 0, count_agree_clear) idx = where(phase_true_sfc gt 0 and phase_retrieved gt 0, count_agree_cloud) idx = where(phase_true_sfc eq 1 and phase_retrieved eq 1, count_agree_water) idx = where(phase_true_sfc eq 2 and phase_retrieved eq 2, count_agree_ice) if (count_all gt 0) then begin pod_clear = float(count_agree_cloud + count_agree_clear) / count_all endif else begin pod_clear = missing endelse if (count_agree_cloud gt 0) then begin pod_phase = float(count_agree_water + count_agree_ice) / count_agree_cloud endif else begin pod_phase = missing endelse print, 'stats for isfc = ', isfc, pod_clear, pod_phase, count_all ;------------------------------------------------------------------------------------ ; make images ;------------------------------------------------------------------------------------ if (~keyword_set(noplot)) then begin p0 = plot(x_x,post_prob_clear, name = 'clear_prob_lidar', line='solid', color='black',thick = 2,$ xtitle = xtitle, xrange = [min(x_x),max(x_x)], $ ytitle = 'frequency', yrange = [0.0, 1.0]) p1 = plot(x_x, post_prob_ice, name = 'ice_prob',line='solid', color='blue', thick = 2, /overplot) p2 = plot(x_x, post_prob_water, name = 'water_prob', line='solid', color='red',thick = 2, /overplot) p3 = plot(x_x, obs_prob_lidar/max(obs_prob_lidar), name = 'obs_prob', line='solid', color='grey',thick = 2, /overplot) leg0 = legend(target=[p0,p1,p2,p3], position = [0.3,0.95], /normal, font_size = 8,/AUTO_TEXT_COLOR) t = text(0.3, 0.8, ['SFC'+string(isfc+1,format="(I02)")], color = "black",/data) p0.save, "poster_prob_1d_"+classifier_name+"_sfc"+string(isfc,format="(I02)")+"_"+day_string+"_"+sensor_name+".png", $ resolution = 100, border=10 endif obs_table[*,isfc] = obs_prob_lidar ice_table[*,isfc] = cond_ratio_ice water_table[*,isfc] = cond_ratio_water clear_table[*,isfc] = cond_ratio_clear cycle_sfc: endfor ;sfc type loop ;------------------------------------------------------- ; constrain ;------------------------------------------------------- idx = where(obs_table lt 0.0 and obs_table ne missing,cc) & if (cc gt 0) then obs_table[idx] = 0.00 idx = where(ice_table lt 0.0 and ice_table ne missing,cc) & if (cc gt 0) then ice_table[idx] = 0.00 idx = where(water_table lt 0.0 and water_table ne missing,cc) & if (cc gt 0) then water_table[idx] = 0.00 idx = where(clear_table lt 0.0 and clear_table ne missing,cc) & if (cc gt 0) then clear_table[idx] = 0.00 ;------------------------------------------------- ; make a structure to hold data and attributes ;------------------------------------------------- ;---- assign attributes @attributes_define.pro attributes.classifier_name = classifier_name attributes.length_classifier_name = strlen(classifier_name) attributes.description = nc_header attributes.sensor = sensor_name attributes.training_source = training_source attributes.timestamp = timestamp() attributes.cod_clear_water_thresh = cod_clear_water_thresh attributes.cod_clear_ice_thresh = cod_clear_ice_thresh attributes.rank = nb_rank attributes.nsfc = nsfc attributes.nchan_used = nchan_used attributes.x_name = xtitle attributes.nbins_x = nbins_x attributes.x_min = x_min attributes.x_bin = x_bin attributes.y_name = ytitle attributes.nbins_y = nbins_y attributes.y_min = y_min attributes.y_bin = y_bin attributes.z_name = ztitle attributes.nbins_z = nbins_z attributes.z_min = z_min attributes.z_bin = z_bin attributes.lat_min = lat_min attributes.lat_max = lat_max attributes.zen_min = zen_min attributes.zen_max = zen_max attributes.solzen_min = solzen_min attributes.solzen_max = solzen_max attributes.solglintzen_min = solglintzen_min attributes.solglintzen_max = solglintzen_max attributes.solscatang_min = solscatang_min attributes.solscatang_max = solscatang_max attributes.lunzen_min = lunzen_min attributes.lunzen_max = lunzen_max attributes.lunglintzen_min = lunglintzen_min attributes.lunglintzen_max = lunglintzen_max attributes.lunscatang_min = lunscatang_min attributes.lunscatang_max = lunscatang_max attributes.tpw_min = slant_tpw_min attributes.tpw_max = slant_tpw_max attributes.tsfc_min = tsfc_min attributes.tsfc_max = tsfc_max attributes.zsfc_min = zsfc_min attributes.zsfc_max = zsfc_max attributes.zsfc_std_min = zsfc_std_min attributes.zsfc_std_max = zsfc_std_max attributes.snow_class_min = snow_class_min attributes.snow_class_max = snow_class_max attributes.solglint_mask_min = solglint_mask_min attributes.solglint_mask_max = solglint_mask_max attributes.lunglint_mask_min = lunglint_mask_min attributes.lunglint_mask_max = lunglint_mask_max attributes.coast_mask_min = coast_mask_min attributes.coast_mask_max = coast_mask_max attributes.city_mask_min = city_mask_min attributes.city_mask_max = city_mask_max attributes.moon_illum_frac_min = moon_illum_frac_min attributes.moon_illum_frac_max = moon_illum_frac_max attributes.rut_solzen_thresh = rut_solzen_thresh ;---- assign table @lut_1d_define.pro lut.clear = clear_table lut.ice = ice_table lut.water = water_table lut.obs = obs_table lut.observation_count = observation_count lut.cloud_fraction = cloud_fraction lut.ice_fraction = ice_fraction lut.water_fraction = water_fraction lut.on_flag = sfc_type_on_flag channels_used.wvl = channel_wvl_used ;---------------------------------------------- ; write to netcdf file ;---------------------------------------------- ;nc_filename = path_dir + sensor_name + "_" + day_string+"_"+classifier_name+".nc" nc_filename = path_dir + sensor_name + "_" + day_string+"_"+test_name+".nc" sd_id = ncdf_create(nc_filename,/CLOBBER,/NETCDF4_FORMAT) write_1d_table, sd_id, lut, attributes, channels_used, success_flag, /class_att ncdf_close, sd_id print, 'write success = ', success_flag skip_all: end