;$Id$ ;----------------------------------------------------------------------------- ; Routines to derive Bayesian Cloud Mask Classifier Files ; ; Requires idl save file from make_xxx.pro ; ; output is a text file with is used by CLAVR-x/PATMOS-x, figures and ; performance metrics ; ;----------------------------------------------------------------------------- @legend_local @tag_exist pro ecm_v1_5_lut, time_diff_max = time_diff_max, SFC_TYPE_INPUT = sfc_type_input, $ SOLZEN_MIN_IN=solzen_min_in, SOLZEN_MAX_IN=solzen_max_in, $ AVHRR1=avhrr1, TRAINING_FILE=training_file, $ Log10_Emiss_Tropo = Log10_Emiss_Tropo, $ Emiss_Tropo = Emiss_Tropo, $ T_11 = T_11, $ T_Std = T_Std, $ LOG10_T_Std = LOG10_T_Std, $ T_Max = T_Max, $ T_Sub_Diff = T_Sub_Diff, $ T_Clear = T_Clear, $ FMFT = FMFT, $ Emiss_375_All = Emiss_375_All, $ Emiss_375_Day = Emiss_375_Day, $ Emiss_375_Night = Emiss_375_Night, $ Emiss_375_Emiss_Tropo_Night = Emiss_375_Emiss_Tropo_Night, $ Log10_Ref_063_Day = Log10_Ref_063_Day, $ Ref_063_Day = Ref_063_Day, $ Ref_DNB_Night = Ref_DNB_Night, $ Ref_063_Min_3x3_Day = Ref_063_Min_3x3_Day, $ Ref_Std_Day = Ref_Std_Day, $ Ref_086_Day = Ref_086_Day, $ Ref_Ratio_Day = Ref_Ratio_Day, $ Ref_138_Day = Ref_138_Day, $ Ref_160_Day = Ref_160_Day, $ Ref_375_Day = Ref_375_Day, $ Ndsi_Day = Ndsi_Day, $ Btd_11_85 = Btd_11_85, $ Btd_375_11_Night = Btd_375_11_Night, $ Btd_375_11_Day = Btd_375_11_Day, $ Btd_375_11_All = Btd_375_11_All, $ Bt_11_67_Covar = Bt_11_67_Covar, $ Btd_11_67 = Btd_11_67, $ Btd_85_73 = Btd_85_73, $ Opd_063_Day = Opd_063_Day, $ Log10_Opd_063_Day = Log10_Opd_063_Day, $ Zcld_Opa = Zcld_Opa, $ output_name = output_name, $ cod_min = cod_min, $ png=png ;----------------------------------------------------------------------------- ; Create and test a naive bayesian classifier for PATMOS-x ; ; 1) Read in and determine bins for the PATMOS-x variables ; ; current list of variables: ; ;time_diff num_cld_layers num_cld_layers_min ; ; 0 = Deep Ocean ; 1 = Shallow Ocean ; 2 = Land ; 3 = Arctic ; 4 = Antarctica ; 5 = Desert ;----------------------------------------------------------------------------- FORWARD_FUNCTION CREATE_BIN, MAKE_T_RATIO if (~keyword_set(output_name)) then output_name = 'temp' if (~keyword_set(solzen_min_in)) then solzen_min_in = 0.00 if (~keyword_set(solzen_max_in)) then solzen_max_in = 180.0 glintzen_vis_max = 180.0 glintzen_vis_min = 40.0 solzen_vis_min = 0.0 solzen_vis_max = 80.0 solzen_swir_day_min = 0.0 solzen_swir_day_max = 90.0 solzen_swir_night_min = 90.0 solzen_swir_night_max = 180.0 lunzen_vis_min = 0.0 lunzen_vis_max = 80.0 lunglintzen_vis_max = 180.0 lunglintzen_vis_min = 40.0 ;---- some parameters missing = -999.0 ncolors = 250 iplot = 1 ; 0 = no plots ioutput = 1 ; 0 = no output n_bounds_reg = 101 total_number_of_sfc_types = 7 skip_sfc_type_flag = make_array(total_number_of_sfc_types+1,/INT, VALUE=0) prior_yes = make_array(total_number_of_sfc_types+1,/FLOAT, VALUE=missing) prior_no = make_array(total_number_of_sfc_types+1,/FLOAT, VALUE=missing) ;----- set flags n_class = keyword_set(Emiss_Tropo) + keyword_set(Log10_Emiss_Tropo) + keyword_set(T_Clear) + $ keyword_set(T_11) + keyword_set(T_Std) + keyword_set(Log10_T_Std) + keyword_set(T_Max) + $ keyword_set(T_Sub_Diff) + $ keyword_set(FMFT) + keyword_set(Emiss_375_All) + keyword_set(Emiss_375_Day) + $ keyword_set(Emiss_375_Night) + keyword_set(Emiss_375_Emiss_Tropo_Night) + $ keyword_set(Opd_063_Day) + keyword_set(Log10_Opd_063_Day) + $ keyword_set(Ref_063_Day) + keyword_set(Ref_063_Min_3x3_Day) + $ keyword_set(Ref_DNB_Night) + $ keyword_set(Ref_Std_Day) + $ keyword_set(Ref_086_Day) + $ keyword_set(Ref_Ratio_Day) + $ keyword_set(Ref_138_Day) + $ keyword_set(Ndsi_Day) + $ keyword_set(Ref_160_Day) + $ keyword_set(Ref_375_Day) + keyword_set(Btd_11_85) + keyword_set(Btd_85_73) + $ keyword_set(Bt_11_67_Covar) + keyword_set(Btd_11_67) + $ keyword_set(Btd_375_11_All) + keyword_set(Btd_375_11_Day) + keyword_set(Btd_375_11_Night) + $ keyword_set(Zcld_Opa) print, "N_Class = ", n_class ;--------------------------------------------------------- ; set missing arguments to missing ;--------------------------------------------------------- if (keyword_set(TRAINING_FILE) eq 0) then begin training_file = 'phase_training.sav' endif if (keyword_set(time_diff_max) eq 0) then begin time_diff_max = 10.0 endif print, 'Maximum allowed time difference (min) = ', time_diff_max if (keyword_set(SFC_TYPE_INPUT) eq 0) then begin sfc_type_input_start = 1 sfc_type_input_end = total_number_of_sfc_types sfc_type_input_name = ['Deep_Water','Shallow_Water','Unfrozen_Land' $ ,'Frozen_Land','Arctic','Antarctic','Desert'] endif else begin sfc_type_input_start = sfc_type_input sfc_type_input_end = sfc_type_input sfc_type_input_name = [' '] endelse ;------------------------------------------------------------------------- ;---- set fidelity of bins ;------------------------------------------------------------------------- min_number_cdf = 100 n_bins_reg = n_bounds_reg - 1 n_bounds = 51 n_bins = n_bounds - 1 X_bounds = make_array(n_bounds,n_class,/FLOAT,VALUE=missing) first_valid_bounds = make_array(n_class,/INT,VALUE=missing) last_valid_bounds = make_array(n_class,/INT,VALUE=missing) X_bin_centers = make_array(n_bins,n_class,/FLOAT,VALUE=missing) X_bin_centers_reg = fltarr(n_bins_reg,n_class)-999.0 X_Name = strarr(n_class) Posterior_By_Class_bin = make_array(n_bins_reg,n_class,/FLOAT,VALUE=missing) delta_bin = make_array(n_class,/FLOAT,VALUE=missing) bin_start = make_array(n_class,/FLOAT,VALUE=missing) bin_end = make_array(n_class,/FLOAT,VALUE=missing) class_cond_yes = make_array(n_bins,n_class,/FLOAT,VALUE=missing) class_cond_no = make_array(n_bins,n_class,/FLOAT,VALUE=missing) X_bounds_reg = make_array(n_bounds_reg,n_class,/FLOAT,VALUE=missing) class_cond_yes_reg = make_array(n_bins_reg,n_class,/FLOAT,VALUE=missing) class_cond_no_reg = make_array(n_bins_reg,n_class,/FLOAT,VALUE=missing) class_cond_ratio_reg = make_array(n_bins_reg,n_class,/FLOAT,VALUE=missing) X_bounds_reg_temp = make_array(n_bounds_reg,/FLOAT,VALUE=missing) class_cond_yes_reg_temp = make_array(n_bins_reg,/FLOAT,VALUE=missing) class_cond_no_reg_temp = make_array(n_bins_reg,/FLOAT,VALUE=missing) ;--- define and set default for limits table_rank = make_array(n_class,/INT,VALUE=1) tsfc_min = make_array(n_class,/FLOAT,VALUE=180.0) tsfc_max = make_array(n_class,/FLOAT,VALUE=340.0) zsfc_min = make_array(n_class,/FLOAT,VALUE=-500.0) zsfc_max = make_array(n_class,/FLOAT,VALUE=9000.0) solzen_min = make_array(n_class,/FLOAT,VALUE=0.0) solzen_max = make_array(n_class,/FLOAT,VALUE=180.0) glintzen_min = make_array(n_class,/FLOAT,VALUE=0.0) glintzen_max = make_array(n_class,/FLOAT,VALUE=180.0) lunzen_min = make_array(n_class,/FLOAT,VALUE=0.0) lunzen_max = make_array(n_class,/FLOAT,VALUE=180.0) lunglintzen_min = make_array(n_class,/FLOAT,VALUE=0.0) lunglintzen_max = make_array(n_class,/FLOAT,VALUE=180.0) coast_min = make_array(n_class,/INT,VALUE=0) coast_max = make_array(n_class,/INT,VALUE=1) sfc_type_flag = make_array(n_class,total_number_of_sfc_types,/FLOAT,VALUE=1) n_sfc_type = total_number_of_sfc_types false_alarm_rate = fltarr(n_sfc_type+1) false_alarm_old_rate = fltarr(n_sfc_type+1) false_alarm_modis_rate = fltarr(n_sfc_type+1) missed_cloud_rate_global = fltarr(n_sfc_type+1) pod_rate = fltarr(n_sfc_type+1) pod_old_rate = fltarr(n_sfc_type+1) pod_modis_rate = fltarr(n_sfc_type+1) skill_rate = fltarr(n_sfc_type+1) skill_old_rate = fltarr(n_sfc_type+1) skill_modis_rate = fltarr(n_sfc_type+1) n_ec_zc_zone = 5 missed_cloud_rate = fltarr(n_ec_zc_zone,n_sfc_type+1) missed_cloud_old_rate = fltarr(n_sfc_type+1) missed_cloud_modis_rate = fltarr(n_sfc_type+1) optimal_pod = fltarr(n_sfc_type+1) optimal_posterior_prob = fltarr(n_sfc_type+1) optimal_skill = fltarr(n_sfc_type+1) global_cloud_frac_calipso = fltarr(n_sfc_type+1) global_cloud_frac_new = fltarr(n_sfc_type+1) global_cloud_frac_old = fltarr(n_sfc_type+1) global_cloud_frac_modis = fltarr(n_sfc_type+1) ;--- set domain domain_string = 'global' lat_north = 90.0 lat_south = -90.0 lon_west = -180.0 lon_east = 180.0 ;--- set thresholds for what is clear and cloudy cld_frac_clear_thresh = 0.0 cld_frac_cloud_thresh = 1.0 if (~keyword_set(cod_min)) then cod_min = 0.01 ;----------------------------------------------------------------------------- ; Read in the IDL save file ;----------------------------------------------------------------------------- data_file = training_file restore, data_file xxx = training_data training_data = !NULL ;-----------------------------------------------------------------------------------------' ; make classifiers needed for mask ;-----------------------------------------------------------------------------------------' ndata = n_elements(xxx.closest_calipso_time_difference) missing_data = make_array(ndata,/FLOAT,VALUE=missing) ;--- check for NaN in median index = where(xxx.refl_0_65um_nom_clear_sky le 0, count) if (count gt 0) then xxx[index].refl_0_65um_nom_clear_sky = 5.0 make_t_ratio, xxx.temp_11_0um_nom, xxx.temp_11_0um_nom_clear_sky, missing, t_ratio make_rfmft, xxx.temp_11_0um_nom, xxx.temp_12_0um_nom, xxx.diff_ch31_ch32_bt_ch31_max_3x3, missing, rfmft_metric ;--- limit t_11um_Std index = where(xxx.temp_11_0um_nom_stddev_3x3 gt 5.0, count) if (count gt 0) then xxx[index].temp_11_0um_nom_stddev_3x3 = 5.0 ;----- make fmft classifier make_fmft, xxx.temp_11_0um_nom, xxx.temp_12_0um_nom, xxx.temp_11_0um_nom_clear_sky, $ xxx.temp_12_0um_nom_clear_sky, missing, fmft_metric ;----- make fmft classifier btd_375_11_metric = (xxx.temp_3_75um_nom - xxx.temp_11_0um_nom) - $ (xxx.temp_3_75um_nom_clear_sky - xxx.temp_11_0um_nom_clear_sky) ;--- T11_max - T11 t_11um_rel = xxx.temp_11_0um_nom_max_3x3 - xxx.temp_11_0um_nom index = where(xxx.temp_11_0um_nom eq missing, count) if (count gt 0) then t_11um_rel[index] = missing index = where(t_11um_rel gt 20.0, count) if (count gt 0) then t_11um_rel[index] = 20 ;--- T11_max_sub - T11_min_sub t_11um_rel_sub = t_11um_rel if (tag_exist(xxx,'temp_11_0um_nom_max_sub') and tag_exist(xxx,'temp_11_0um_nom_min_sub')) then begin t_11um_rel_sub = xxx.temp_11_0um_nom_max_sub - xxx.temp_11_0um_nom_min_sub index = where(xxx.temp_11_0um_nom_max_sub eq missing or xxx.temp_11_0um_nom_min_sub eq missing, count) if (count gt 0) then t_11um_rel_sub[index] = missing index = where(t_11um_rel_sub gt 20.0, count) if (count gt 0) then t_11um_rel_sub[index] = 20 endif ;--- Ref_063_Day metric rgct = xxx.refl_0_65um_nom - xxx.refl_0_65um_nom_clear_sky index = where(xxx.refl_0_65um_nom le 0.0 or xxx.refl_0_65um_nom_clear_sky le 0.0, count) if (count gt 0) then rgct[index] = missing ;--- Ref_DNB_Night metric dnb_refl_test = make_array(ndata, /FLOAT, VALUE=missing) if (tag_exist(xxx,'refl_lunar_dnb_nom')) then begin dnb_refl_test = xxx.refl_lunar_dnb_nom index = where(xxx.refl_lunar_dnb_nom le 0.0, count) if (count gt 0) then dnb_refl_test[index] = missing endif rvct = xxx.refl_0_65um_nom - xxx.refl_0_65um_nom_min_3x3 index = where(xxx.refl_0_65um_nom le 0.0 or xxx.refl_0_65um_nom_min_3x3 le 0.0, count) if (count gt 0) then rvct[index] = missing ;--- reflectance ratio rrct = xxx.refl_0_86um_nom / xxx.refl_0_65um_nom index = where(xxx.refl_0_86um_nom le 0.0 or xxx.refl_0_65um_nom le 0.0, count) if (count gt 0) then rrct[index] = missing ;--- ndsi ndsi = (xxx.refl_0_65um_nom - xxx.refl_1_60um_nom) / (xxx.refl_0_65um_nom + xxx.refl_1_60um_nom) index = where(xxx.refl_1_60um_nom le 0.0 or xxx.refl_0_65um_nom le 0.0, count) if (count gt 0) then ndsi[index] = missing ;--- Ref_375_Day metric ref_375_metric = xxx.refl_3_75um_nom index = where(xxx.refl_3_75um_nom le 0.0, count) if (count gt 0) then ref_375_metric[index] = missing ;---------------------------------------------------------------------- ; Determine Surface Type ;---------------------------------------------------------------------- ;--- comment out the following to not use CLAVR-x values in file sfc_type_mask = xxx.bayes_mask_sfc_type goto, skip_sfc_type_creation make_sfc_type_mask, xxx.land_class, xxx.snow_class, xxx.latitude, xxx.longitude, $ xxx.emiss_3_75um_nom, xxx.sst_background_uni_3x3,xxx.surface_type, $ sfc_type_mask skip_sfc_type_creation: ;============== ;window, 0, xsize = 800, ysize = 800 ;load_color_table, 33, 1.0, 0, 249 ;map_set, 0,0 ;erase, color = 0 ;map_continents, color = 4 ;for i = 0L, n_elements(xxx.land_class) -1 do begin ;if (sfc_type_mask[i] eq 6) then begin ; print, xxx[i].lon, xxx[i].lat, xxx[i].cld_frac_lidar ; plots, xxx[i].lon, xxx[i].lat, psym = 3, color = 3 ;endif ;endfor ;stop ;============== ;------------------------------------------------------------------------------------------------- ; print out distribution of pixels in each surface type ;------------------------------------------------------------------------------------------------- count_total = n_elements(xxx.closest_calipso_time_difference) for isfc = 0,total_number_of_sfc_types-1 do begin print, 'number of pixels = ', n_elements(where(sfc_type_mask eq isfc+1)), ' for sfc type = ', isfc+1 print, '% of pixels = ', 100.0*n_elements(where(sfc_type_mask eq isfc+1))/float(count_total), ' for sfc type = ', isfc+1 endfor ;------------------------------------------------------------------------------------------------- ;----- make a binary calipso cld mask (0=clear, 3=cloud) ;------------------------------------------------------------------------------------------------- cld_mask = xxx.cloud_mask cld_mask_modis = xxx.cloud_mask ; --- disabled this for now, because cloud_mask_aux exist but all values set to missing (Denis B.) ;if (tag_exist(xxx,'cloud_mask_aux')) then cld_mask_modis = xxx.cloud_mask_aux cld_mask_calipso = cld_mask cld_mask_calipso[*] = missing if (tag_exist(xxx,'closest_calipso_cloud_fraction_5km')) then calipso_cloud_fraction = xxx.closest_calipso_cloud_fraction_5km if (tag_exist(xxx,'closest_calipso_cloud_fraction_x5')) then calipso_cloud_fraction = xxx.closest_calipso_cloud_fraction_x5 surface_temperature = xxx.TEMP_11_0UM_NOM_CLEAR_SKY if (tag_exist(xxx,'surface_temperature_nwp')) then surface_temperature = xxx.surface_temperature_nwp index = where(calipso_cloud_fraction le cld_frac_clear_thresh and $ calipso_cloud_fraction ne missing, count) if (count gt 0) then cld_mask_calipso[index] = 0 index = where(calipso_cloud_fraction ge cld_frac_cloud_thresh and $ calipso_cloud_fraction ne missing, count) if (count gt 0) then cld_mask_calipso[index] = 3 ;------------------------------------------------------------------------------------------------- ;--- force very thin cirrus to be called clear ;------------------------------------------------------------------------------------------------- index = where(xxx.closest_calipso_cod ne missing and xxx.closest_calipso_cod lt cod_min, count) if (count gt 0) then begin cld_mask_calipso[index] = 0 calipso_cloud_fraction[index] = 0.0 endif ;---------------------------------------------------------------------------------------- ; 3 point smoothing of calipso cloud mask ;---------------------------------------------------------------------------------------- n = n_elements(cld_mask_calipso) temp = cld_mask_calipso for i = 1L, n - 2L do begin cld_mask_calipso[i] = mean(temp[i-1:i+1]) endfor ;------------------------------------------------------------------------- ; skip pixels that don't seem to match - obvious cloud missed in calipso ; set cld_mask_calipso to 1, not 0 or 3 which are used. ;-------------------------------------------------------------------------- index_skip = where(cld_mask_calipso lt 2 and $ xxx.emiss_tropo_11_0um_nom gt 0.4 and $ sfc_type_mask ne 4 and $ sfc_type_mask ne 5 and $ sfc_type_mask ne 6,count_skip) if (count_skip gt 0) then cld_mask_calipso[index_skip] = 1 index_skip = where(cld_mask_calipso lt 2 and $ xxx.refl_0_65um_nom gt 50.0 and $ sfc_type_mask ne 4 and $ sfc_type_mask ne 5 and $ sfc_type_mask ne 6,count_skip) if (count_skip gt 0) then cld_mask_calipso[index_skip] = 1 ;--------------------------------------------------------------------- ; open a file for writing the results ;--------------------------------------------------------------------- if (ioutput eq 1) then begin openw, lun, output_name+'.txt', /get_lun printf, lun, data_file printf, lun, format = "(a10, 6f10.2)", $ domain_string, lat_north, lat_south, $ lon_west, lon_east, solzen_min_in, solzen_max_in printf, lun, time_diff_max, n_class, n_bounds_reg, n_sfc_type svn_id_string = "$Id$" cdfid = ncdf_create(output_name+'.nc',/NETCDF4_FORMAT,/CLOBBER);,/NETCDF4_FORMAT) ncdf_attput, cdfid, /GLOBAL, 'data_file', data_file ncdf_attput, cdfid, /GLOBAL, 'timestamp', timestamp(), /CHAR ncdf_attput, cdfid, /GLOBAL, 'svn_id', svn_id_string, /CHAR ncdf_attput, cdfid, /GLOBAL, 'domain_string', domain_string ;ncdf_attput, cdfid, /GLOBAL, 'training_data_timestamp', training_data_creation_time_string, /CHAR ;ncdf_attput, cdfid, /GLOBAL, 'training_data_timestamp', training_creation_time, /CHAR ;ncdf_attput, cdfid, /GLOBAL, 'training_data_code_svn_id', training_data_code_id_string, /CHAR ncdf_attput, cdfid, /GLOBAL, 'lat_north', lat_north ncdf_attput, cdfid, /GLOBAL, 'lat_south', lat_south ncdf_attput, cdfid, /GLOBAL, 'lon_west', lon_west ncdf_attput, cdfid, /GLOBAL, 'lon_east', lon_east ncdf_attput, cdfid, /GLOBAL, 'solzen_min', solzen_min_in ncdf_attput, cdfid, /GLOBAL, 'solzen_max', solzen_max_in ncdf_attput, cdfid, /GLOBAL, 'time_diff_max', time_diff_max ncdf_attput, cdfid, /GLOBAL, 'n_class', n_class ncdf_attput, cdfid, /GLOBAL, 'n_bounds_reg', n_bounds_reg ncdf_attput, cdfid, /GLOBAL, 'n_bins_reg', n_bins_reg ncdf_attput, cdfid, /GLOBAL, 'n_sfc_type', n_sfc_type endif ;------------------------------------------------------------------- ; loop through sfc type generate the needed table values ;------------------------------------------------------------------ count_valid_sum = 0L for sfc_type_input = sfc_type_input_start,sfc_type_input_end do begin ;--- assign header name based on sfc_type if (sfc_type_input eq 1) then header = 'Deep_Water' if (sfc_type_input eq 2) then header = 'Shallow_Water' if (sfc_type_input eq 3) then header = 'Unfrozen_Land' if (sfc_type_input eq 4) then header = 'Frozen_Land' if (sfc_type_input eq 5) then header = 'Arctic' if (sfc_type_input eq 6) then header = 'Antarctic' if (sfc_type_input eq 7) then header = 'Desert' ;------------------------------------------------------------------- ; select data for analysis ;------------------------------------------------------------------- index_valid = where(sfc_type_mask eq sfc_type_input $ and abs(xxx.closest_calipso_time_difference) lt time_diff_max and $ xxx.latitude le lat_north and $ xxx.latitude ge lat_south and $ xxx.longitude ge lon_west and $ xxx.longitude le lon_east and $ surface_temperature ne missing and $ xxx.emiss_tropo_11_0um_nom ne missing and $ (rgct ne missing or xxx.solar_zenith_angle gt solzen_vis_max) and $ xxx.emiss_3_75um_nom ne missing and xxx.emiss_3_75um_nom_clear ne missing and $ t_11um_rel ne missing and $ xxx.solar_zenith_angle ge solzen_min_in and xxx.solar_zenith_angle le solzen_max_in and $ ;---> xxx.glint_zenith_angle ge glintzen_vis_min and xxx.glint_zenith_angle le glintzen_vis_max and $ (cld_mask_calipso eq 0 or cld_mask_calipso eq 3) $ , count_valid) print, 'number of pixels after filtering = ', count_valid count_valid_sum = count_valid_sum + count_valid print, 'accumulated number of pixels after filtering = ', count_valid_sum skip_sfc_type_flag[sfc_type_input] = 0 if (count_valid eq 0) then begin print, 'Not enough valid points' prior_yes[sfc_type_input] = missing prior_no[sfc_type_input] = missing skip_sfc_type_flag[sfc_type_input] = 1 goto, skip_sfc_type endif ;------------------------------------------------------------------------------------ ; Assign the classes ;------------------------------------------------------------------------------------ yes_ = where(cld_mask_calipso[index_valid] eq 3,n_yes_) no_ = where(cld_mask_calipso[index_valid] eq 0,n_no_) X = fltarr(count_valid,n_class) Y = fltarr(count_valid,n_class) ;---- set defaults on limits solzen_min[*] = 0.0 solzen_max[*] = 180.0 glintzen_min[*] = 0.0 glintzen_max[*] = 180.0 sfc_type_flag[*,*] = 1 ;these are the internal surface types tsfc_min[*] = 180.0 tsfc_max[*] = 340.0 zsfc_min[*] = -500.0 zsfc_max[*] = 10000.0 coast_min[*] = 0 coast_max[*] = 1 iclass = 0 ;--- note, limits have default values. Only set those that are needed if keyword_set(Emiss_Tropo) then begin X[*,iclass] = xxx[index_valid].emiss_tropo_11_0um_nom & X_Name[iclass] = 'Emiss_Tropo' tsfc_min[iclass] = 230.0 & tsfc_max[iclass] = 340.0 iclass = iclass + 1 endif if keyword_set(Log10_Emiss_Tropo) then begin make_log, xxx[index_valid].emiss_tropo_11_0um_nom, 0.001, 0.0, missing, xtemp X[*,iclass] = xtemp & X_Name[iclass] = 'Log10_Emiss_Tropo' tsfc_min[iclass] = 230.0 & tsfc_max[iclass] = 340.0 iclass = iclass + 1 endif if keyword_set(Zcld_Opa) then begin X[*,iclass] = xxx[index_valid].cld_height_opaque & X_Name[iclass] = 'Zcld_Opa' sfc_type_flag[iclass,*]=[1,1,1,1,1,1,1] iclass = iclass + 1 endif if keyword_set(Emiss_375_All) then begin X[*,iclass] = (xxx[index_valid].ems_375um-xxx[index_valid].emiss_3_75um_nom_clear) / $ xxx[index_valid].emiss_3_75um_nom_clear & X_Name[iclass] = 'Emiss_375' iclass = iclass + 1 endif if keyword_set(Emiss_375_Emiss_Tropo_Night) then begin X[*,iclass] = xxx[index_valid].ems_375um * xxx[index_valid].emiss_tropo_11_0um_nom & X_Name[iclass] = 'Emiss_375_Emiss_Tropo_Night' solzen_min[iclass] = solzen_swir_night_min & solzen_max[iclass] = solzen_swir_night_max iclass = iclass + 1 endif if keyword_set(Emiss_375_Night) then begin X[*,iclass] = xxx[index_valid].emiss_3_75um_nom & X_Name[iclass] = 'Emiss_375_Night' ; X[*,iclass] = (ems20[index_valid]-ems20_clr[index_valid])/ems20_clr[index_valid] & X_Name[iclass] = 'Emiss_375_Night' solzen_min[iclass] = solzen_swir_night_min & solzen_max[iclass] = solzen_swir_night_max glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = glintzen_vis_max tsfc_min[iclass] = 200.0 & tsfc_max[iclass] = 340.0 iclass = iclass + 1 endif if keyword_set(Emiss_375_Day) then begin X[*,iclass] = (xxx[index_valid].emiss_3_75um_nom - xxx[index_valid].emiss_3_75um_nom_clear) / $ xxx[index_valid].emiss_3_75um_nom_clear & X_Name[iclass] = 'Emiss_375_Day' solzen_min[iclass] = solzen_swir_day_min & solzen_max[iclass] = solzen_swir_day_max glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = glintzen_vis_max iclass = iclass + 1 endif if keyword_set(Ref_063_Day) then begin X[*,iclass] = rgct[index_valid] & X_Name[iclass] = 'Ref_063_Day' solzen_min[iclass] = solzen_vis_min & solzen_max[iclass] = solzen_vis_max glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = glintzen_vis_max iclass = iclass + 1 endif if keyword_set(Ref_DNB_Night) then begin X[*,iclass] = dnb_refl_test[index_valid] & X_Name[iclass] = 'Ref_DNB_Night' lunzen_min[iclass] = lunzen_vis_min & lunzen_max[iclass] = lunzen_vis_max lunglintzen_min[iclass] = lunglintzen_vis_min & lunglintzen_max[iclass] = lunglintzen_vis_max sfc_type_flag[iclass,*]=[1,1,1,0,0,0,1] iclass = iclass + 1 endif if keyword_set(Ref_Std_Day) then begin X[*,iclass] = xxx[index_valid].refl_0_65um_nom_stddev_3x3 & X_Name[iclass] = 'Ref_Std_Day' solzen_min[iclass] = solzen_vis_min & solzen_max[iclass] = solzen_vis_max glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = glintzen_vis_max coast_min[iclass] = 0 & coast_max[iclass] = 0 iclass = iclass + 1 endif if keyword_set(Ref_063_Min_3x3_Day) then begin X[*,iclass] = rvct[index_valid] & X_Name[iclass] = 'Ref_063_Min_3x3_Day' solzen_min[iclass] = solzen_vis_min & solzen_max[iclass] = solzen_vis_max glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = glintzen_vis_max coast_min[iclass] = 0 & coast_max[iclass] = 0 iclass = iclass + 1 endif if keyword_set(Ref_086_Day) then begin X[*,iclass] = xxx[index_valid].refl_0_86um_nom & X_Name[iclass] = 'Ref_086_Day' solzen_min[iclass] = solzen_vis_min & solzen_max[iclass] = solzen_vis_max glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = glintzen_vis_max iclass = iclass + 1 endif if keyword_set(Ref_Ratio_Day) then begin X[*,iclass] = rrct[index_valid] & X_Name[iclass] = 'Ref_Ratio_Day' solzen_min[iclass] = solzen_vis_min & solzen_max[iclass] = solzen_vis_max glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = glintzen_vis_max iclass = iclass + 1 endif if keyword_set(Ref_138_Day) then begin X[*,iclass] = xxx[index_valid].refl_1_38um_nom & X_Name[iclass] = 'Ref_138_Day' ii = where(X[*,iclass] gt 30.0,count) if (count gt 0) then X[ii,iclass] = 30.0 solzen_min[iclass] = solzen_vis_min & solzen_max[iclass] = solzen_vis_max glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = glintzen_vis_max iclass = iclass + 1 endif if keyword_set(Ref_160_Day) then begin X[*,iclass] = xxx[index_valid].refl_1_60um_nom & X_Name[iclass] = 'Ref_160_Day' solzen_min[iclass] = solzen_vis_min & solzen_max[iclass] = solzen_vis_max glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = glintzen_vis_max iclass = iclass + 1 endif if keyword_set(Ref_375_Day) then begin X[*,iclass] = ref_375_metric[index_valid] & X_Name[iclass] = 'Ref_375_Day' solzen_min[iclass] = solzen_vis_min & solzen_max[iclass] = solzen_vis_max glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = glintzen_vis_max iclass = iclass + 1 endif if keyword_set(Ndsi_Day) then begin X[*,iclass] = ndsi[index_valid] & X_Name[iclass] = 'Ndsi_Day' solzen_min[iclass] = solzen_vis_min & solzen_max[iclass] = solzen_vis_max glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = glintzen_vis_max iclass = iclass + 1 endif if keyword_set(Opd_063_Day) then begin X[*,iclass] = xxx[index_valid].cld_opd_mask_0_65um_nom & X_Name[iclass] = 'Opd_063_Day' solzen_min[iclass] = solzen_vis_min & solzen_max[iclass] = solzen_vis_max glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = glintzen_vis_max iclass = iclass + 1 endif if keyword_set(Log10_Opd_063_Day) then begin make_log, xxx[index_valid].cld_opd_mask_0_65um_nom, 0.01, 21.0, missing, xtemp X[*,iclass] = xtemp & X_Name[iclass] = 'Log10_Opd_063_Day' solzen_min[iclass] = solzen_vis_min & solzen_max[iclass] = solzen_vis_max glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = glintzen_vis_max iclass = iclass + 1 endif if keyword_set(T_11) then begin X[*,iclass] = xxx[index_valid].temp_11_0um_nom & X_Name[iclass] = 'T_11' iclass = iclass + 1 endif if keyword_set(T_Std) then begin X[*,iclass] = xxx[index_valid].temp_11_0um_nom_stddev_3x3 & X_Name[iclass] = 'T_Std' coast_min[iclass] = 0 & coast_max[iclass] = 0 iclass = iclass + 1 endif if keyword_set(LOG10_T_Std) then begin make_log, xxx[index_valid].temp_11_0um_nom_stddev_3x3, 0.001, 0.0, missing, xtemp X[*,iclass] = xtemp & X_Name[iclass] = 'Log10_T_Std' coast_min[iclass] = 0 & coast_max[iclass] = 0 iclass = iclass + 1 endif if keyword_set(T_max) then begin X[*,iclass] = t_11um_rel[index_valid] & X_Name[iclass] = 'T_Max-T' iclass = iclass + 1 endif if keyword_set(T_Sub_Diff) then begin X[*,iclass] = t_11um_rel_sub[index_valid] & X_Name[iclass] = 'T_Sub_Diff' iclass = iclass + 1 endif if keyword_set(T_Clear) then begin X[*,iclass] = xxx[index_valid].temp_11_0um_nom - xxx[index_valid].temp_11_0um_nom_clear_sky & X_Name[iclass] = 'T-T_Clear' iclass = iclass + 1 endif if keyword_set(FMFT) then begin X[*,iclass] = fmft_metric[index_valid] & X_Name[iclass] = 'FMFT' iclass = iclass + 1 endif if keyword_set(Btd_11_85) then begin X[*,iclass] = xxx[index_valid].temp_11_0um_nom - xxx[index_valid].temp_8_5um_nom & X_Name[iclass] = 'Btd_11_85' iclass = iclass + 1 endif if keyword_set(Btd_85_73) then begin X[*,iclass] = xxx[index_valid].temp_8_5um_nom - xxx[index_valid].temp_7_3um_nom & X_Name[iclass] = 'Btd_85_73' iclass = iclass + 1 endif if keyword_set(Btd_375_11_All) then begin X[*,iclass] = btd_375_11_metric[index_valid] & X_Name[iclass] = 'Btd_375_11_All' glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = 180.0 iclass = iclass + 1 endif if keyword_set(Btd_375_11_Day) then begin X[*,iclass] = btd_375_11_metric[index_valid] & X_Name[iclass] = 'Btd_375_11_Day' solzen_min[iclass] = solzen_swir_day_min & solzen_max[iclass] = solzen_swir_day_max glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = glintzen_vis_max iclass = iclass + 1 endif if keyword_set(Btd_375_11_Night) then begin X[*,iclass] = btd_375_11_metric[index_valid] & X_Name[iclass] = 'Btd_375_11_Night' solzen_min[iclass] = solzen_swir_night_min & solzen_max[iclass] = solzen_swir_night_max glintzen_min[iclass] = glintzen_vis_min & glintzen_max[iclass] = glintzen_vis_max iclass = iclass + 1 endif if keyword_set(Bt_11_67_Covar) then begin X[*,iclass] = xxx[index_valid].temp_11um_vs_67um_covar_5x5 & X_Name[iclass] = 'Bt_11_67_Covar' iclass = iclass + 1 endif if keyword_set(Btd_11_67) then begin X[*,iclass] = xxx[index_valid].temp_11_0um_nom - xxx[index_valid].temp_6_7um_nom & X_Name[iclass] = 'Btd_11_67' iclass = iclass + 1 endif ;-------------------------------------------------------------------------------------------- ; apply each classifiers limits by setting out of limit values as missing ; if filter is missing, ignore the filtering ;-------------------------------------------------------------------------------------------- for iclass = 0, n_class -1 do begin ;--- solzen if (tag_exist(xxx,'solar_zenith_angle')) then begin filter_metric = xxx[index_valid].solar_zenith_angle ii = where(filter_metric ne missing and filter_metric gt solzen_max[iclass],cc) & if (cc gt 0) then X[ii,iclass] = missing ii = where(filter_metric ne missing and filter_metric lt solzen_min[iclass],cc) & if (cc gt 0) then X[ii,iclass] = missing endif ;--- glintzen if (tag_exist(xxx,'glint_zenith_angle')) then begin filter_metric = xxx[index_valid].glint_zenith_angle ii = where(filter_metric ne missing and filter_metric gt glintzen_max[iclass],cc) & if (cc gt 0) then X[ii,iclass] = missing ii = where(filter_metric ne missing and filter_metric lt glintzen_min[iclass],cc) & if (cc gt 0) then X[ii,iclass] = missing endif ;--- lunzen if (tag_exist(xxx,'lunar_zenith_angle')) then begin filter_metric = xxx[index_valid].lunar_zenith_angle ii = where(filter_metric ne missing and filter_metric gt lunzen_max[iclass],cc) & if (cc gt 0) then X[ii,iclass] = missing ii = where(filter_metric ne missing and filter_metric lt lunzen_min[iclass],cc) & if (cc gt 0) then X[ii,iclass] = missing endif ;--- lunar glintzen if (tag_exist(xxx,'lunar_glint_zenith_angle')) then begin filter_metric = xxx[index_valid].lunar_glint_zenith_angle ii = where(filter_metric ne missing and filter_metric gt lunglintzen_max[iclass],cc) & if (cc gt 0) then X[ii,iclass] = missing ii = where(filter_metric ne missing and filter_metric lt lunglintzen_min[iclass],cc) & if (cc gt 0) then X[ii,iclass] = missing endif ;--- tsfc if (tag_exist(xxx,'surface_temperature_nwp')) then begin filter_metric = xxx[index_valid].surface_temperature_nwp ii = where(filter_metric ne missing and filter_metric gt tsfc_max[iclass],cc) & if (cc gt 0) then X[ii,iclass] = missing ii = where(filter_metric ne missing and filter_metric lt tsfc_min[iclass],cc) & if (cc gt 0) then X[ii,iclass] = missing endif ;--- zsfc if (tag_exist(xxx,'surface_elevation')) then begin filter_metric = xxx[index_valid].surface_elevation ii = where(filter_metric ne missing and filter_metric gt zsfc_max[iclass],cc) & if (cc gt 0) then X[ii,iclass] = missing ii = where(filter_metric ne missing and filter_metric lt zsfc_min[iclass],cc) & if (cc gt 0) then X[ii,iclass] = missing endif ;--- coast if (tag_exist(xxx,'coast_mask')) then begin filter_metric = xxx[index_valid].coast_mask ii = where(filter_metric ne missing and filter_metric gt coast_max[iclass],cc) & if (cc gt 0) then X[ii,iclass] = missing ii = where(filter_metric ne missing and filter_metric lt coast_min[iclass],cc) & if (cc gt 0) then X[ii,iclass] = missing endif endfor ;------------------------------------------------------------------------------------------------ ; make a string of classifiers for netcdf file ;------------------------------------------------------------------------------------------------ Classifier_String = '' for iclass = 0, n_class -1 do begin if (iclass gt 0) then begin Classifier_String = Classifier_String + ":" + X_Name[iclass] endif else begin Classifier_String = X_Name[iclass] endelse endfor print,typename(Classifier_String) print,Classifier_String ;------------------------------------------------------------------------------------------------- ; take out correlation ;------------------------------------------------------------------------------------------------- i_class_ref = -1 ncorr = 10000 for i_class = 0, n_class -1 do begin for j_class = 0, n_class -1 do begin XX = X[*,i_class] YY = X[*,j_class] idx = where(XX ne missing and YY ne missing, npoint) if (i_class eq i_class_ref and i_class ne j_class) then begin if (npoint gt 0) then begin nstride = max([1,npoint / ncorr]) coef = linfit(xx[0:npoint-1:nstride],yy[0:npoint-1:nstride]) X[*,j_class] = YY - (coef[0] + coef[1]*XX) endif endif endfor endfor ;------------------------------------------------------------------------------------------------- ; correlation of classifiers ;------------------------------------------------------------------------------------------------- Correlation_Matrix = make_array(n_class,n_class,/FLOAT,value=missing) ncorr = 10000 for i_class = 0, n_class -1 do begin for j_class = 0, n_class -1 do begin XX = X[*,i_class] YY = X[*,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 print, format="(13a15)", " ", X_Name for j_class = 0, n_class -1 do begin print, format="(a15, 12F15.2)", X_Name[j_class], Correlation_Matrix[*,j_class] endfor ;------------------------------------------------------------------------------------ ; Determine the bins for each PATMOS-x variable ;------------------------------------------------------------------------------------ ;min_number_cdf = 100 ;n_bins_reg = n_bounds_reg - 1 ;n_bounds = 51 ;n_bins = n_bounds - 1 ;X_bounds = fltarr(n_bounds,n_class)-999.0 ;first_valid_bounds = intarr(n_class) ;last_valid_bounds = intarr(n_class) ;X_bin_centers = fltarr(n_bins,n_class)-999.0 ;X_bin_centers_reg = fltarr(n_bins_reg,n_class)-999.0 for iclass = 0, n_class-1 do begin ii = where(X[*,iclass] ne missing,cc) if (cc lt min_number_cdf) then begin print,'Number of good pixels = ',cc print, 'Insufficient points for cdf' X_bounds[*,iclass] = missing goto, skip_cdf endif ; determine 1% and 99% values nbins_cdf = 1000 determine_cdf_bounds, X[ii,iclass], nbins_cdf, data_cdf_min, data_cdf_max, cdf_fail if (cdf_fail eq 1) then begin data_cdf_min = min(X[ii,iclass]) data_cdf_max = max(X[ii,iclass]) endif ii = where(X[*,iclass] ne missing and X[*,iclass] gt data_cdf_min and X[*,iclass] lt data_cdf_max,cc) X_bounds[*,iclass] = create_bin(X[ii,iclass],n_bins,missing) skip_cdf: endfor for i = 0, n_class-1 do begin for j = 0,n_bins -1 do begin X_bin_centers[j,i] = 0.5*(X_bounds[j,i]+X_bounds[j+1,i]) endfor endfor ;------------------------------------------------------------------------------------ ; The Prior Probability (avg cloudiness) ;------------------------------------------------------------------------------------ prior_yes[sfc_type_input] = n_elements(yes_)/float(count_valid) ;prior_yes[sfc_type_input] = 0.5 prior_no[sfc_type_input] = 1.0 - prior_yes[sfc_type_input] print, "Prior Yes = ", prior_yes[sfc_type_input] ;prior_yes[sfc_type_input] = 0.5 ;prior_no[sfc_type_input] = 0.5 ;------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------ ; Class-conditional probability - YES ;------------------------------------------------------------------------------------ class_cond_yes[*,*] = 0 for iclass = 0, n_class-1 do begin ii = where(X[yes_,iclass] ne missing,cc) for j = 0, n_bins-1 do begin this_class = $ where(X[yes_[ii],iclass] ge X_bounds[j,iclass] and X[yes_[ii],iclass] lt X_bounds[j+1,iclass],nC) if(j eq n_bins-1) then begin this_class = $ where(X[yes_[ii],iclass] ge X_bounds[j,iclass] and $ X[yes_[ii],iclass] le X_bounds[j+1,iclass],nC) endif if(nC gt 0)then class_cond_yes[j,iclass] = float(nC)/float(n_elements(yes_[ii])) endfor ;class_cond_yes[*,iclass] = ts_smooth(class_cond_yes[*,iclass],3) endfor ;------------------------------------------------------------------------------------ ; Class-conditional probability - NO ;------------------------------------------------------------------------------------ class_cond_no[*,*] = 0 for iclass = 0, n_class-1 do begin ii = where(X[no_,iclass] ne missing,cc) for j = 0, n_bins-1 do begin this_class = $ where(X[no_[ii],iclass] ge X_bounds[j,iclass] and X[no_[ii],iclass] lt X_bounds[j+1,iclass],nC) if (j eq n_bins-1) then begin this_class = $ where(X[no_[ii],iclass] ge X_bounds[j,iclass] and $ X[no_[ii],iclass] le X_bounds[j+1,iclass],nC) endif if(nC gt 0) then class_cond_no[j,iclass] = float(nC)/float(n_elements(no_[ii])) endfor ;class_cond_no[*,iclass] = ts_smooth(class_cond_no[*,iclass],3) endfor ;--------------------------------------------------------------------- ; Convert to Regular Spacing ;--------------------------------------------------------------------- class_cond_yes_reg[*,*] = 0.0 class_cond_no_reg[*,*] = 0.0 class_cond_ratio_reg[*,*] = 1.0 for iclass = 0, n_class-1 do begin ii = where(class_cond_yes[*,iclass] ge 0.0,cc) if (cc ge 1) then begin make_regularly_spaced_data,n_bounds,n_bounds_reg, $ X_bounds[*,iclass], class_cond_yes[*,iclass], class_cond_no[*,iclass], $ X_bounds_reg_temp, class_cond_yes_reg_temp, class_cond_no_reg_temp X_bounds_reg[*,iclass] = X_bounds_reg_temp class_cond_yes_reg[*,iclass] = class_cond_yes_reg_temp class_cond_no_reg[*,iclass] = class_cond_no_reg_temp for j = 0,n_bins_reg -1 do begin X_bin_centers_reg[j,iclass] = 0.5*(X_bounds_reg[j,iclass]+X_bounds_reg[j+1,iclass]) endfor endif endfor ;------------------------------------------------------------------------------------- ; make ratio = no / yes ;------------------------------------------------------------------------------------- class_cond_ratio_reg[*,*] = 1.0 idx = where(class_cond_no_reg eq 0.0,cc) if (cc gt 0) then class_cond_ratio_reg[idx] = 0.0 idx = where(class_cond_no_reg ne 0.0 and class_cond_yes_reg eq 0.0,cc) if (cc gt 0) then class_cond_ratio_reg[idx] = 100.0 idx = where(class_cond_no_reg ne 0.0 and class_cond_yes_reg ne 0.0,cc) if (cc gt 0) then class_cond_ratio_reg[idx] = $ class_cond_no_reg[idx] / class_cond_yes_reg[idx] ;------------------------------------------------------------------------------------- ; store first and last valid bounds ;------------------------------------------------------------------------------------- for i = 0, n_class-1 do begin index_min = where(X_bounds_reg[*,i] eq min(X_bounds_reg[*,i]),count_min) first_valid_bounds[i] = count_min - 1 index_max = where(X_bounds_reg[*,i] eq max(X_bounds_reg[*,i]),count_max) last_valid_bounds[i] = n_bounds_reg - count_max endfor ;------------------------------------------------------------------------------------ ; ;------------------------------------------------------------------------------------- for j = 0, n_class -1 do begin delta_bin[j] = X_bounds_reg[1,j] - X_bounds_reg[0,j] bin_start[j] = X_bounds_reg[0,j] bin_end[j] = X_bounds_reg[n_bounds_reg-1,j] if (delta_bin[j] le 0.0) then begin class_cond_yes_reg[*,j] = missing class_cond_no_reg[*,j] = missing endif endfor ;------------------------------------------------------------------------------------ ; Posterior probability -- (prio.y*cond.y)/(prio.y*cond.y+prio.n*cond.n) ;------------------------------------------------------------------------------------ Posterior = fltarr(count_valid) Posterior_By_Class = fltarr(count_valid,n_class) ;--- the bin # for each observations cond_yes = fltarr(count_valid,n_class) cond_no = fltarr(count_valid,n_class) cond_ratio = fltarr(count_valid,n_class) for i = 0l, count_valid-1l do begin for iclass = 0, n_class-1 do begin if (X[i,iclass] ne missing) then begin ibin = fix((X[i,iclass] - bin_start[iclass])/delta_bin[iclass]) ibin = max([0,min([ibin,n_bins_reg-1])]) cond_yes[i,iclass] = class_cond_yes_reg[ibin,iclass] cond_no[i,iclass] = class_cond_no_reg[ibin,iclass] cond_ratio[i,iclass] = class_cond_ratio_reg[ibin,iclass] ; if (cond_no[i,iclass] eq 0.0) then cond_ratio[i,iclass] = 0.0 else $ ; if (cond_yes[i,iclass] eq 0.0) then cond_ratio[i,iclass] = 100.0 else $ ; cond_ratio[i,iclass] = cond_no[i,iclass] / cond_yes[i,iclass] ; print, 'test ', cond_yes[i,iclass],cond_no[i,iclass], cond_ratio[i,iclass] endif else begin cond_yes[i,iclass] = 1.0 cond_no[i,iclass] = 1.0 cond_ratio[i,iclass] = 1.0 endelse ;-- turn off test if AVHRR1 if (keyword_set(AVHRR1) ne 0) then begin if (X_Name[iclass] eq 'FMFT') then begin cond_yes[i,iclass] = 1.0 cond_no[i,iclass] = 1.0 endif if (X_Name[iclass] eq 'RFMFT') then begin cond_yes[i,iclass] = 1.0 cond_no[i,iclass] = 1.0 endif endif endfor endfor ;--- determine the post. prob. for i = 0l, count_valid-1l do begin ;--- all classes combined using naive formulation r = product(cond_ratio[i,*]) Posterior[i] = 1.0 / ( 1.0 + r/prior_yes[sfc_type_input] - r) ;--- do it by class for j = 0, n_class -1 do begin r = cond_ratio[i,j] Posterior_By_Class[i,j] = 1.0 / ( 1.0 + r/prior_yes[sfc_type_input] - r) endfor endfor ;------------------------------------------------------------------------------------ ; make a posterior of each bin for plotting below ;------------------------------------------------------------------------------------ for iclass = 0, n_class -1 do begin for ibin = 0,n_bins_reg-1 do begin cond_yes_bin = class_cond_yes_reg[ibin, iclass] cond_no_bin = class_cond_no_reg[ibin,iclass] r = class_cond_ratio_reg[ibin,iclass] Posterior_By_Class_Bin[ibin,iclass] = 1.0 / ( 1.0 + r/prior_yes[sfc_type_input] - r) endfor endfor ;----------------------------------------------------------------------------- ; rank importance ;----------------------------------------------------------------------------- Prob_Corr = make_array(n_class) for iclass = 0, n_class - 1 do begin f = Posterior_By_Class[*,iclass] idx = where(Posterior ge 0.0 and f ge 0.0, cc) Prob_Corr[iclass] = correlate(Posterior[idx], f[idx]) print, X_Name[iclass], Prob_Corr[iclass] endfor ;----------------------------------------------------------------------------- if (iplot eq 1) then begin load_color_table, 13, 0, 1.0, ncolors xsize = 800 ysize = 400 for iclass = 0,n_class-1 do begin window, iclass, xsize = xsize, ysize=ysize erase, color = 3 plot, X_bin_centers_reg[*,iclass], Posterior_By_Class_bin[*,iclass], color=0,background=3, $ xtitle = X_Name[iclass], $ ytitle = 'Probability', yrange = [0,1.0], ystyle = 1, $ charsize=1.5, xticklen=0.05,yticklen=0.05, thick=2, psym = -4 oplot, X_bin_centers_reg[*,iclass], (class_cond_yes_reg[*,iclass]/max(class_cond_yes_reg[*,iclass])), color=255, thick=2 oplot, X_bin_centers_reg[*,iclass], (class_cond_no_reg[*,iclass]/max(class_cond_no_reg[*,iclass])), color=80, thick=2 ; legend_local, ['Posterior (This Class Alone)', 'Class Conditional Cloudy (Scaled)', 'Class Conditional Clear (Scaled)'], $ ; color=[0,255,80],line=[0,0,0],thick=[2,2,2], /center,/right, charsize=1.5 legend_local, [header], charsize=1.5, /top, /left if (keyword_set(png) ne 0) then saveimage, 'posterior_class'+strtrim(iclass,2)+'_sfc'+strtrim(sfc_type_input,2)+'.tif',/tif endfor endif ;------------------------------------------------------------------------------------ ; Skill Metrics, Precision - Recall ;------------------------------------------------------------------------------------ precision = fltarr(100) recall = fltarr(100) skill = fltarr(100) pod = fltarr(100) posterior_threshold = fltarr(100) new_cld = intarr(count_valid) old_cld = calipso_cloud_fraction[index_valid] for i = 0, 99 do begin new_cld[*] = 0 posterior_threshold[i] = float(i)/100.0 index = where(Posterior ge posterior_threshold[i], gc) if(gc gt 0)then new_cld[index] = 3 ;skill score a = float(n_elements(where(new_cld eq 3 and cld_mask_calipso[index_valid] eq 3 ))) b = float(n_elements(where(new_cld eq 3 and cld_mask_calipso[index_valid] eq 0 ))) c = float(n_elements(where(new_cld eq 0 and cld_mask_calipso[index_valid] eq 3 ))) d = float(n_elements(where(new_cld eq 0 and cld_mask_calipso[index_valid] eq 0 ))) skill[i] = (a*d - b*c) / ( (a+c)*(b+d)) pod[i] = (a+d) / (a+b+c+d) endfor ;---- determine optimal posterior ii = where(skill eq max(skill)) ;ii = where(pod eq max(pod)) optimal_posterior_prob[sfc_type_input] = posterior_threshold[min(ii)] optimal_skill[sfc_type_input] = skill[min(ii)] optimal_pod[sfc_type_input] = pod[min(ii)] if (iplot eq 1) then begin window, 10, xsize=600, ysize=400 load_color_table, 13, 0, 1.0, ncolors erase, color=3 plot,posterior_threshold,skill,yrange=[0.0,1],xrange=[0,1], color = 0, background=3, $ xtitle = 'Posterior Probability', ytitle = 'Hansen-Kuiper Skill Score', charsize=1 oplot,posterior_threshold,pod,linestyle=2 legend_local, ['Skill', 'POD'], lines=[0,2], /center, box=1, charsize=1 legend_local, ['CLASSIFIERS:', X_Name], /left, /bottom, charsize=1 legend_local, ['optimal posterior probability = '+ $ string(optimal_posterior_prob[sfc_type_input],format='(f5.3)'), $ 'maximum skill score = ' + string(optimal_skill[sfc_type_input], format='(f5.3)'), $ 'maximum pod score = ' + string(optimal_pod[sfc_type_input], format='(f5.3)')], $ /right, /bottom, charsize = 1.1 legend_local, [header], charsize=1.5, /top, /left if (keyword_set(png) ne 0) then saveimage, 'posterior_sfc'+strtrim(sfc_type_input,2)+'.tif',/tif endif ;--------------------------------------------------------- ; output classifiers ;--------------------------------------------------------- skip_sfc_type: if (skip_sfc_type_flag[sfc_type_input] eq 1) then begin prior_yes[sfc_type_input] = missing prior_no[sfc_type_input] = missing optimal_posterior_prob[sfc_type_input] = missing first_valid_bounds[*] = missing last_valid_bounds[*] = missing bin_start[*] = missing bin_end[*] = missing delta_bin[*] = missing class_cond_yes_reg[*,*] = missing class_cond_no_reg[*,*] = missing endif X_Name_output = X_Name X_Name_length = make_array(n_class, /INT, VALUE=0) string_length = 30 blank_strings = make_array(string_length,/STRING,VALUE=' ') blank_string = strjoin(blank_strings) for i_class = 0, n_class - 1 do begin X_Name_length[i_class] = strlen(strtrim(X_Name[i_class])) string_temp = strtrim(X_Name[i_class]) X_Name_output[i_class] = string_temp + strmid(blank_string,0,string_length - strlen(string_temp)) endfor if (ioutput eq 1) then begin ;------------------------------------------------- ; netcdf output ;------------------------------------------------- if (sfc_type_input eq 1) then begin ncdf_attput, cdfid, /GLOBAL, 'classifier_string', Classifier_String add_class_flags, cdfid, N_Class, X_Name ;--- set dimentions nclassid = NCDF_DIMDEF(cdfid, 'nclass', n_class) ; Define the iclass dimension. stringid = NCDF_DIMDEF(cdfid, 'str_length', 30) sfc_type_id = NCDF_DIMDEF(cdfid, 'num_sfc_type', n_sfc_type) n_bounds_id = NCDF_DIMDEF(cdfid, 'n_bounds', n_bounds_reg-1) endif ;------------------------------------------------- ; txt output ;------------------------------------------------- printf, lun, sfc_type_input, ' ',header printf, lun, prior_yes[sfc_type_input], prior_no[sfc_type_input] printf, lun, optimal_posterior_prob[sfc_type_input] for iclass = 0, n_class -1 do begin printf, lun, iclass+1, first_valid_bounds[iclass]+1,last_valid_bounds[iclass]+1, $ ' ',X_Name[iclass] ; note the +1 n first and last valid bounds printf, lun, bin_start[iclass],bin_end[iclass], delta_bin[iclass] printf, lun, class_cond_yes_reg[*,iclass] printf, lun, class_cond_no_reg[*,iclass] endfor endif ;------------------------------------------------------------------------------ ; compute some statistics ;------------------------------------------------------------------------------ if (count_valid gt 0) then begin ;--- make a cloud mask from the posterior probabilities cld_mask_new = bytarr(count_valid) cld_mask_new[*] = missing posterior_thresh = optimal_posterior_prob[sfc_type_input] posterior_thresh = 0.5 index = where(Posterior ge posterior_thresh and Posterior ne missing, count) if (count gt 0) then cld_mask_new[index] = 3 index = where(Posterior lt posterior_thresh and Posterior ne missing, count) if (count gt 0) then cld_mask_new[index] = 0 ;--- bulk stats cld_mask_x = cld_mask_calipso[index_valid] cld_mask_y = cld_mask_new a = float(n_elements(where(cld_mask_y eq 3 and cld_mask_x eq 3 ))) b = float(n_elements(where(cld_mask_y eq 3 and cld_mask_x eq 0 ))) c = float(n_elements(where(cld_mask_y eq 0 and cld_mask_x eq 3 ))) d = float(n_elements(where(cld_mask_y eq 0 and cld_mask_x eq 0 ))) false_alarm_rate[sfc_type_input] = b / (a + b + c + d) missed_cloud_rate_global[sfc_type_input] = c / (a + b + c + d) skill_rate[sfc_type_input] = (a*d - b*c) / ( (a+c)*(b+d)) pod_rate[sfc_type_input] = (a+d) / (a+b+c+d) posterior_yes = float(n_elements(where(cld_mask_y eq 3))) / (a + b + c + d) global_cloud_frac_calipso[sfc_type_input] = prior_yes[sfc_type_input] global_cloud_frac_new[sfc_type_input] = posterior_yes ;------------------------------------------------- ; compare against old mask ; allow cm = 0 or 1 for old to be clear ; allow cm = 2 or 3 for old to be cloudy ;------------------------------------------------- cld_mask_x = cld_mask_calipso[index_valid] cld_mask_y = cld_mask[index_valid] a = float(n_elements(where(cld_mask_y ge 2 and cld_mask_x eq 3 ))) b = float(n_elements(where(cld_mask_y ge 2 and cld_mask_x eq 0 ))) c = float(n_elements(where(cld_mask_y le 1 and cld_mask_x eq 3 ))) d = float(n_elements(where(cld_mask_y le 1 and cld_mask_x eq 0 ))) skill_old_rate[sfc_type_input] = (a*d - b*c) / ( (a+c)*(b+d)) pod_old_rate[sfc_type_input] = (a+d) / (a+b+c+d) posterior_yes = float(n_elements(where(cld_mask_y ge 2))) / (a + b + c + d) global_cloud_frac_old[sfc_type_input] = posterior_yes false_alarm_old_rate[sfc_type_input] = b / (a + b + c + d) missed_cloud_old_rate[sfc_type_input] = c / (a + b + c + d) ;------------------------------------------------- ; compare against modis mask ; allow cm = 0 or 1 for old to be clear ; allow cm = 2 or 3 for old to be cloudy ;------------------------------------------------- cld_mask_x = cld_mask_calipso[index_valid] cld_mask_y = cld_mask_modis[index_valid] a = float(n_elements(where(cld_mask_y ge 2 and cld_mask_x eq 3 ))) b = float(n_elements(where(cld_mask_y ge 2 and cld_mask_x eq 0 ))) c = float(n_elements(where(cld_mask_y le 1 and cld_mask_x eq 3 ))) d = float(n_elements(where(cld_mask_y le 1 and cld_mask_x eq 0 ))) skill_modis_rate[sfc_type_input] = (a*d - b*c) / ( (a+c)*(b+d)) pod_modis_rate[sfc_type_input] = (a+d) / (a+b+c+d) posterior_yes = float(n_elements(where(cld_mask_y ge 2))) / (a + b + c + d) global_cloud_frac_modis[sfc_type_input] = posterior_yes false_alarm_modis_rate[sfc_type_input] = b / (a + b + c + d) missed_cloud_modis_rate[sfc_type_input] = c / (a + b + c + d) endif endfor ; loop over surfaces if (ioutput eq 1) then begin close, lun free_lun, lun ;--- set variables vid0 = NCDF_VARDEF(cdfid, 'sfc_type_names', [stringid,sfc_type_id], /CHAR) vid1 = NCDF_VARDEF(cdfid, 'classifier_names', [stringid,nclassid], /CHAR) vid2 = NCDF_VARDEF(cdfid, 'classifier_names_length', [nclassid], /LONG) vid3 = NCDF_VARDEF(cdfid, 'prior_yes', [sfc_type_id], /FLOAT) vid4 = NCDF_VARDEF(cdfid, 'prior_no', [sfc_type_id], /FLOAT) vid5 = NCDF_VARDEF(cdfid, 'optimal_posterior_prob', [sfc_type_id], /FLOAT) vid6 = NCDF_VARDEF(cdfid, 'first_valid_bounds', [nclassid,sfc_type_id], /SHORT) vid7 = NCDF_VARDEF(cdfid, 'last_valid_bounds', [nclassid,sfc_type_id], /SHORT) vid8 = NCDF_VARDEF(cdfid, 'bin_start', [nclassid,sfc_type_id], /FLOAT) vid9 = NCDF_VARDEF(cdfid, 'bin_end', [nclassid,sfc_type_id], /FLOAT) vid10 = NCDF_VARDEF(cdfid, 'delta_bin', [nclassid,sfc_type_id], /FLOAT) vid11 = NCDF_VARDEF(cdfid, 'class_cond_yes_reg', [n_bounds_id,nclassid,sfc_type_id], /FLOAT) vid12 = NCDF_VARDEF(cdfid, 'class_cond_no_reg', [n_bounds_id,nclassid,sfc_type_id], /FLOAT) vid13 = NCDF_VARDEF(cdfid, 'class_cond_ratio_reg', [n_bounds_id,nclassid,sfc_type_id], /FLOAT) vid14 = NCDF_VARDEF(cdfid, 'pc', [sfc_type_id], /FLOAT) vid15 = NCDF_VARDEF(cdfid, 'skill', [sfc_type_id], /FLOAT) vid16 = NCDF_VARDEF(cdfid, 'missed', [sfc_type_id], /FLOAT) vid17 = NCDF_VARDEF(cdfid, 'false', [sfc_type_id], /FLOAT) vid18 = NCDF_VARDEF(cdfid, 'solzen_min', [nclassid], /FLOAT) vid19 = NCDF_VARDEF(cdfid, 'solzen_max', [nclassid], /FLOAT) vid20 = NCDF_VARDEF(cdfid, 'glintzen_min', [nclassid], /FLOAT) vid21 = NCDF_VARDEF(cdfid, 'glintzen_max', [nclassid], /FLOAT) vid22 = NCDF_VARDEF(cdfid, 'lunzen_min', [nclassid], /FLOAT) vid23 = NCDF_VARDEF(cdfid, 'lunzen_max', [nclassid], /FLOAT) vid24 = NCDF_VARDEF(cdfid, 'lunglintzen_min', [nclassid], /FLOAT) vid25 = NCDF_VARDEF(cdfid, 'lunglintzen_max', [nclassid], /FLOAT) vid26 = NCDF_VARDEF(cdfid, 'tsfc_min', [nclassid], /FLOAT) vid27 = NCDF_VARDEF(cdfid, 'tsfc_max', [nclassid], /FLOAT) vid28 = NCDF_VARDEF(cdfid, 'zsfc_min', [nclassid], /FLOAT) vid29 = NCDF_VARDEF(cdfid, 'zsfc_max', [nclassid], /FLOAT) vid30 = NCDF_VARDEF(cdfid, 'sfc_type_flag', [nclassid,sfc_type_id], /FLOAT) vid31 = NCDF_VARDEF(cdfid, 'coast_min', [nclassid], /SHORT) vid32 = NCDF_VARDEF(cdfid, 'coast_max', [nclassid], /SHORT) vid33 = NCDF_VARDEF(cdfid, 'table_rank', [nclassid], /SHORT) ;--- write data NCDF_CONTROL, cdfid, /ENDEF ; Put the file into data mode NCDF_VARPUT, cdfid, vid0, sfc_type_input_name NCDF_VARPUT, cdfid, vid1, X_Name_output NCDF_VARPUT, cdfid, vid2, X_Name_length NCDF_VARPUT, cdfid, vid3, prior_yes NCDF_VARPUT, cdfid, vid4, prior_no NCDF_VARPUT, cdfid, vid5, optimal_posterior_prob NCDF_VARPUT, cdfid, vid6, first_valid_bounds NCDF_VARPUT, cdfid, vid7, last_valid_bounds NCDF_VARPUT, cdfid, vid8, bin_start NCDF_VARPUT, cdfid, vid9, bin_end NCDF_VARPUT, cdfid, vid10, delta_bin NCDF_VARPUT, cdfid, vid11, class_cond_yes_reg NCDF_VARPUT, cdfid, vid12, class_cond_no_reg NCDF_VARPUT, cdfid, vid13, class_cond_ratio_reg NCDF_VARPUT, cdfid, vid14, pod_rate[1:total_number_of_sfc_types] NCDF_VARPUT, cdfid, vid15, skill_rate[1:total_number_of_sfc_types] NCDF_VARPUT, cdfid, vid16, missed_cloud_rate_global[1:total_number_of_sfc_types] NCDF_VARPUT, cdfid, vid17, false_alarm_rate[1:total_number_of_sfc_types] NCDF_VARPUT, cdfid, vid18, solzen_min NCDF_VARPUT, cdfid, vid19, solzen_max NCDF_VARPUT, cdfid, vid20, glintzen_min NCDF_VARPUT, cdfid, vid21, glintzen_max NCDF_VARPUT, cdfid, vid22, lunzen_min NCDF_VARPUT, cdfid, vid23, lunzen_max NCDF_VARPUT, cdfid, vid24, lunglintzen_min NCDF_VARPUT, cdfid, vid25, lunglintzen_max NCDF_VARPUT, cdfid, vid26, tsfc_min NCDF_VARPUT, cdfid, vid27, tsfc_max NCDF_VARPUT, cdfid, vid28, zsfc_min NCDF_VARPUT, cdfid, vid29, zsfc_max NCDF_VARPUT, cdfid, vid30, sfc_type_flag NCDF_VARPUT, cdfid, vid31, coast_min NCDF_VARPUT, cdfid, vid32, coast_max NCDF_VARPUT, cdfid, vid33, table_rank ; --- close netcdf file ncdf_close, cdfid endif print, format = "(a20,10f8.3)",'CALIPSO Cloud Frac = ', global_cloud_frac_calipso(1:total_number_of_sfc_types) print, format = "(a20,10f8.3)",'New Cloud Frac = ', global_cloud_frac_new(1:total_number_of_sfc_types) print, format = "(a20,10f8.3)",'Old Cloud Frac = ', global_cloud_frac_old(1:total_number_of_sfc_types) print, format = "(a20,10f8.3)",'MODIS Cloud Frac = ', global_cloud_frac_modis(1:total_number_of_sfc_types) print, '-------------------------------------------------------------' print, format = "(a20,10f8.3)",'optimal post prob = ', optimal_posterior_prob(1:total_number_of_sfc_types) print, '-------------------------------------------------------------' print, format = "(a20,10f8.3)", 'POD New = ', pod_rate(1:total_number_of_sfc_types) print, format = "(a20,10f8.3)", 'Skill Score New = ', skill_rate[1:total_number_of_sfc_types] print, format = "(a20,10f8.3)", 'False Alarm New = ', false_alarm_rate[1:total_number_of_sfc_types] print, format = "(a20,10f8.3)", 'Missed Cloud New = ', missed_cloud_rate_global[1:total_number_of_sfc_types] print, '-------------------------------------------------------------' print, format = "(a20,10f8.3)", 'POD Old = ', pod_old_rate(1:total_number_of_sfc_types) print, format = "(a20,10f8.3)", 'Skill Score Old = ', skill_old_rate(1:total_number_of_sfc_types) print, format = "(a20,10f8.3)", 'False Alarm Old = ', false_alarm_old_rate(1:total_number_of_sfc_types) print, format = "(a20,10f8.3)", 'Missed Cloud Old = ', missed_cloud_old_rate(1:total_number_of_sfc_types) print, '-------------------------------------------------------------' print, format = "(a20,10f8.3)", 'POD MODIS = ', pod_modis_rate(1:total_number_of_sfc_types) print, format = "(a20,10f8.3)", 'Skill Score Modis = ', skill_modis_rate(1:total_number_of_sfc_types) print, format = "(a20,10f8.3)", 'False Alarm Modis = ', false_alarm_modis_rate(1:total_number_of_sfc_types) print, format = "(a20,10f8.3)", 'Missed Cloud Modis = ', missed_cloud_modis_rate(1:total_number_of_sfc_types) end ;--------------------------------------------------------- ; resample class conditional distributions to a constant grid ;--------------------------------------------------------- pro make_regularly_spaced_data,n_bounds,n_bounds_reg, $ X_bounds, class_cond_yes, class_cond_no, $ X_bounds_reg, class_cond_yes_reg, class_cond_no_reg class_cond_yes_cdf_reg = fltarr(n_bounds_reg) class_cond_no_cdf_reg = fltarr(n_bounds_reg) class_cond_yes_cdf = fltarr(n_bounds) class_cond_no_cdf = fltarr(n_bounds) xmin = X_bounds[0] xmax = X_bounds[n_bounds-1] dx = (xmax -xmin)/(n_bounds_reg-1) X_bounds_reg = xmin + dx*findgen(n_bounds_reg) X_bin_centers_reg = xmin + dx*findgen(n_bounds_reg-1) + dx/2.0 ;for a cdf of raw distribution class_cond_yes_cdf[0] = 0.0 class_cond_no_cdf[0] = 0.0 for j = 1,n_bounds-1 do begin class_cond_yes_cdf[j] = class_cond_yes_cdf[j-1] + class_cond_yes[j-1] class_cond_no_cdf[j] = class_cond_no_cdf[j-1] + class_cond_no[j-1] endfor ;form a regularly spaced cdf class_cond_yes_cdf_reg[0] = class_cond_yes_cdf[0] class_cond_no_cdf_reg[0] = class_cond_no_cdf[0] class_cond_yes_cdf_reg[n_bounds_reg-1] = class_cond_yes_cdf[n_bounds-1] class_cond_no_cdf_reg[n_bounds_reg-1] = class_cond_no_cdf[n_bounds-1] for j = 1,n_bounds_reg - 2 do begin class_cond_yes_cdf_reg[j] = interpol(class_cond_yes_cdf[*],X_bounds[*],X_bounds_reg[j]) class_cond_no_cdf_reg[j] = interpol(class_cond_no_cdf[*],X_bounds[*],X_bounds_reg[j]) endfor ;convert cdfs into regularly spaced class conditionals for j = 0,n_bounds_reg - 2 do begin class_cond_yes_reg[j] = class_cond_yes_cdf_reg[j+1] - class_cond_yes_cdf_reg[j] class_cond_no_reg[j] = class_cond_no_cdf_reg[j+1] - class_cond_no_cdf_reg[j] endfor end ;--------------------------------------------------------------------- function create_bin, var, nv, missing ;----------------------------------------------------------------------------- ; find bin limits for variable [var] with every nv-th percentile of the pdf ;----------------------------------------------------------------------------- ;nv = 50 index = where(var ne missing, count) min_v = min(var[index]) max_v = max(var[index]) binsize = (max_v-min_v)/float(nv) hist_v = histogram(var,min=min_v,max=max_v-binsize,binsize=binsize,locations=locations) hist_v = hist_v/float(total(hist_v)) cdf_v = fltarr(n_elements(hist_v)) for i = 0, n_elements(hist_v)-1 do cdf_v[i] = total(hist_v[0:i]) boundaries = [interpol(locations,cdf_v,findgen(nv)*(1.0/nv)),max_v] neg = where(boundaries lt min_v, nN) if(nN gt 0)then boundaries[neg] = min_v ;--- check for zeros on the bottom ;index = where(boundaries eq min_v,count) ;if (count gt 1) then begin ;nbottom = return,boundaries end ;-------------------------------------------------------------------------------- ; generic netcdf write routine - assumes no scaling ;-------------------------------------------------------------------------------- pro write_nc, cdfid, varname, data, long_name ;--- read data data_type = size(data,/type) ;--- define variable case data_type of 1: varid = ncdf_vardef(cdfid,varname,/byte) 2: varid = ncdf_vardef(cdfid,varname,/short) 3: varid = ncdf_vardef(cdfid,varname,/long) 4: varid = ncdf_vardef(cdfid,varname,/float) 7: varid = ncdf_vardef(cdfid,varname,/char) endcase ;--- write long name ncdf_attput, cdfid, varid, 'long_name', long_name ;--- write data ncdf_control, cdfid, /ENDEF ncdf_varput, cdfid, varid, data ncdf_control, cdfid, /REDEF end ;-------------------------------------------------------------------------------- ; make bayes mask surface type ;-------------------------------------------------------------------------------- pro make_sfc_type_mask, land_class, snow_class, latitude, longitude, emiss_3_75um_nom, $ sst_background_uni_3x3, sfc_type, sfc_type_mask ;--- land class definitions SHALLOW_OCEAN = 0 LAND = 1 COASTLINE = 2 SHALLOW_INLAND_WATER = 3 EPHEMERAL_WATER = 4 DEEP_INLAND_WATER = 5 MODERATE_OCEAN = 6 DEEP_OCEAN = 7 ; snow class definitions NO_SNOW = 1 SEA_ICE = 2 SNOW = 3 sfc_type_mask = make_array(n_elements(land_class),/BYTE,VALUE=0) for i = 0L, n_elements(land_class) -1 do begin ;--- default LAND/OCEAN values if (land_class[i] eq LAND) then begin sfc_type_mask[i] = 3 endif else begin sfc_type_mask[i] = 1 endelse ;--- #2 - Shallow Ocean if (land_class[i] eq MODERATE_OCEAN or $ land_class[i] eq DEEP_INLAND_WATER or $ land_class[i] eq SHALLOW_INLAND_WATER or $ land_class[i] eq SHALLOW_OCEAN) then begin sfc_type_mask[i] = 2 endif ;--- include into Shallow Ocean SST fronts if (land_class[i] ne LAND and $ sst_background_uni_3x3[i] gt 0.5) then begin sfc_type_mask[i] = 1 endif ;--- #3 Unfrozen Land if (land_class[i] eq LAND or $ land_class[i] eq COASTLINE or $ land_class[i] eq EPHEMERAL_WATER) then begin sfc_type_mask[i] = 3 endif ;--- #4 Frozen Land if (snow_class[i] eq SNOW and latitude[i] gt -60.0) then begin sfc_type_mask[i] = 4 endif ;--- #5 - Arctic if (latitude[i] gt 0 and snow_class[i] eq SEA_ICE) then begin sfc_type_mask[i] = 5 endif ;---- #6 - Antarctica if (latitude[i] lt 0 and snow_class[i] eq SEA_ICE) then begin sfc_type_mask[i] = 6 endif if (latitude[i] lt -60 and snow_class[i] eq SNOW) then begin sfc_type_mask[i] = 6 endif ;--- include greenland (add zsfc?) if (latitude[i] gt 60 and longitude[i] gt -75.0 and longitude[i] lt -10 and $ (land_class[i] eq LAND or land_class[i] eq COASTLINE) and $ snow_class[i] eq SNOW) then begin sfc_type_mask[i] = 6 endif ;--- #7 Desert - subset of unfrozen land if (sfc_type_mask[i] eq 3 and abs(latitude[i]) lt 60.0) then begin if(emiss_3_75um_nom[i] lt 0.93) then begin if(sfc_type[i] eq 8 or sfc_type[i] eq 9 or sfc_type[i] eq 10 or sfc_type[i] eq 12) then begin sfc_type_mask[i] = 7 endif endif endif endfor print, 'ENDING sfc type definition' end ;-------------------------------------------------------------------------------- ; t_ratio classifier ;-------------------------------------------------------------------------------- pro make_t_ratio, t11, t11_clr, missing, t_ratio t_ratio = (t11_clr - t11) / (t11 - 180.0) index = where(t_ratio gt 2, count) if (count gt 0) then t_ratio[index] = 2.0 end ;-------------------------------------------------------------------------------- ; rfmft classifier ;-------------------------------------------------------------------------------- pro make_rfmft, t11, t12, btd_11_12_max_11,missing, rfmft rfmft = (t11 - t12) - (btd_11_12_max_11) index = where(t11 eq missing or $ t12 eq missing or $ btd_11_12_max_11 eq missing or $ abs(t11 - t12) gt 20.0, count) if (count gt 0) then rfmft[index] = missing index = where(rfmft gt 5.0, count) if (count gt 0) then rfmft[index] = 5.0 index = where(rfmft lt -5.0, count) if (count gt 0) then rfmft[index] = -5.0 end ;-------------------------------------------------------------------------------- ; fmft classifier ;-------------------------------------------------------------------------------- pro make_fmft, t11, t12, t11_clr, t12_clr, missing, fmft fmft = (t11_clr - t12_clr) * (t11 - 260.0)/(t11_clr - 260.0) fmft[where(t11_clr le 265)] = 0.0 fmft = (t11 - t12) - fmft index = where(t11 eq missing or t12 eq missing or $ t11_clr eq missing or t12_clr eq missing or $ abs(t11 - t12) gt 20.0, count) if (count gt 0) then fmft[index] = missing end ;-------------------------------------------------------------------------------- ; determine 1% and 99% cdf limits for making classifier bins ;-------------------------------------------------------------------------------- pro determine_cdf_bounds, data_cdf, nbins_cdf, data_cdf_min, data_cdf_max, cdf_fail cdf_fail = 0 bin_min_cdf = min(data_cdf) bin_max_cdf = max(data_cdf) bin_size_cdf = (bin_max_cdf - bin_min_cdf) / nbins_cdf x_cdf = bin_min_cdf + bin_size_cdf*findgen(nbins_cdf) + bin_size_cdf/2.0 histo_pdf = histogram(data_cdf,min=bin_min_cdf,nbins=nbins_cdf) histo_cdf = fltarr(nbins_cdf) histo_cdf[0] = histo_pdf[0] for ibin_cdf = 1,nbins_cdf-1 do begin histo_cdf[ibin_cdf] = histo_cdf[ibin_cdf-1] + histo_pdf[ibin_cdf] endfor histo_cdf = histo_cdf / max(histo_cdf) index = where(histo_cdf lt 0.99,cc) if (cc gt 0) then begin data_cdf_max = x_cdf[max(index)] endif else begin print, 'cdf maximum limit not determined' data_cdf_max = bin_max_cdf cdf_fail = 1 endelse index = where(histo_cdf gt 0.01,cc) if (cc gt 0) then begin data_cdf_min = x_cdf[min(index)] endif else begin print, 'cdf minimum limit not determined' data_cdf_min = bin_min_cdf cdf_fail = 1 endelse end ;-------------------------------------------------------- ; ;-------------------------------------------------------- pro make_log, data, data_floor, data_offset, missing, log_data log_data = data idx = where(data ne missing, cc) if (cc gt 0) then data[idx] = data[idx] + data_offset idx = where(data le 0.0 and data ne missing, cc) if (cc gt 0) then data[idx] = data_floor idx = where(data gt 0.0, cc, complement=idx_bad, ncomplement = n_bad) if (cc gt 0) then begin if (data_offset le 0.0) then begin log_data[idx] = alog10(data[idx]) endif else begin log_data[idx] = alog10(data[idx]) - alog10(data_offset) endelse endif if (n_bad gt 0) then log_data[idx_bad] = missing end ;------------------------------------------------- ; set flags for which tests are on ; helps in table usage, -1 = classifier missing ;------------------------------------------------- pro add_class_flags, cdfid, N_Class, X_Name Log10_Emiss_Tropo_Class_Idx = -1 Emiss_Tropo_Class_Idx = -1 T_11_Class_Idx = -1 T_Std_Class_Idx = -1 LOG10_T_Std_Class_Idx = -1 T_Max_Class_Idx = -1 T_Clear_Class_Idx = -1 FMFT_Class_Idx = -1 Emiss_375_All_Class_Idx = -1 Emiss_375_Day_Class_Idx = -1 Emiss_375_Night_Class_Idx = -1 Emiss_375_Emiss_Tropo_Class_Idx = -1 Log10_Ref_063_Day_Class_Idx = -1 Ref_063_Day_Class_Idx = -1 Ref_063_Min_3x3_Day_Class_Idx = -1 Ref_Std_Day_Class_Idx = -1 Ref_086_Day_Class_Idx = -1 Ref_Ratio_Day_Class_Idx = -1 Ref_138_Day_Class_Idx = -1 Ref_160_Day_Class_Idx = -1 Ref_375_Day_Class_Idx = -1 Ndsi_Day_Class_Idx = -1 Btd_11_85_Class_Idx = -1 Btd_375_11_Night_Class_Idx = -1 Btd_375_11_Day_Class_Idx = -1 Btd_375_11_All_Class_Idx = -1 Bt_11_67_Covar_Class_Idx = -1 Btd_11_67_Class_Idx = -1 Btd_85_73_Class_Idx = -1 Opd_063_Day_Class_Idx = -1 Log10_Opd_063_Day_Class_Idx = -1 Zcld_Opa_Class_Idx = -1 ;--- make classifer mapping - maybe add to Lut? for Class_Idx = 0, N_Class - 1 do begin case X_Name[Class_Idx] of "Log10_Emiss_Tropo": Log10_Emiss_Tropo_Class_Idx = Class_Idx "Emiss_Tropo": Emiss_Tropo_Class_Idx = Class_Idx "T_11": T_11_Class_Idx = Class_Idx "T_Std": T_Std_Class_Idx = Class_Idx "Log10_T_Std": Log10_T_Std_Class_Idx = Class_Idx "T_Max-T": T_Max_Class_Idx = Class_Idx "T_Clear-T": T_Clear_Class_Idx = Class_Idx "FMFT": FMFT_Class_Idx = Class_Idx "Emiss_375_All": Emiss_375_All_Class_Idx = Class_Idx "Emiss_375_Day": Emiss_375_Day_Class_Idx = Class_Idx "Emiss_375_Night": Emiss_375_Night_Class_Idx = Class_Idx "Emiss_375_Emiss_Tropo": Emiss_375_Emiss_Tropo_Class_Idx = Class_Idx "Log10_Ref_063_Day": Log10_Ref_063_Day_Class_Idx = Class_Idx "Ref_063_Day": Ref_063_Day_Class_Idx = Class_Idx "Ref_063_Min_3x3_Day": Ref_063_Min_3x3_Day_Class_Idx = Class_Idx "Ref_Std_Day": Ref_Std_Day_Class_Idx = Class_Idx "Ref_086_Day": Ref_086_Day_Class_Idx = Class_Idx "Ref_Ratio_Day": Ref_Ratio_Day_Class_Idx = Class_Idx "Ref_138_Day": Ref_138_Day_Class_Idx = Class_Idx "Ref_375_Day": Ref_375_Day_Class_Idx = Class_Idx "Ndsi_Day": Ndsi_Day_Class_Idx = Class_Idx "Btd_11_85": Btd_11_85_Class_Idx = Class_Idx "Btd_375_11_Night": Btd_375_11_Night_Class_Idx = Class_Idx "Btd_375_11_Day": Btd_375_11_Day_Class_Idx = Class_Idx "Btd_375_11_All": Btd_375_11_Day_Class_Idx = Class_Idx "Bt_11_67_Covar": Bt_11_67_Covar_Class_Idx = Class_Idx "Btd_11_67": Btd_11_67_Class_Idx = Class_Idx "Btd_85_73": Btd_85_73_Class_Idx = Class_Idx "Opd_063_Day": Opd_063_Day_Class_Idx = Class_Idx "Log10_Opd_063_Day": Log10_Opd_063_Day_Class_Idx = Class_Idx "Zcld_Opd_Day": Zcld_Opa_Class_Idx = Class_Idx endcase endfor ncdf_attput, cdfid, /GLOBAL, 'log10_emiss_tropo_class_idx', Log10_Emiss_Tropo_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'emiss_tropo_class_idx', Emiss_Tropo_Class_Idx ncdf_attput, cdfid, /GLOBAL, 't_11_class_idx', T_11_Class_Idx ncdf_attput, cdfid, /GLOBAL, 't_std_class_idx', T_Std_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'log10_t_std_class_idx', Log10_T_Std_Class_Idx ncdf_attput, cdfid, /GLOBAL, 't_max-t_class_idx', T_Max_Class_Idx ncdf_attput, cdfid, /GLOBAL, 't_clear-t_class_idx', T_Clear_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'fmft_class_idx', FMFT_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'emiss_375_all_class_idx', Emiss_375_All_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'emiss_375_day_class_idx', Emiss_375_Day_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'emiss_375_night_class_idx', Emiss_375_Night_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'emiss_375_emiss_tropo_class_idx', Emiss_375_Emiss_Tropo_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'log10_ref_063_day_class_idx', Log10_Ref_063_Day_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'ref_063_day_class_idx', Ref_063_Day_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'ref_063_min_3x3_day_class_idx', Ref_063_Min_3x3_Day_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'ref_std_day_class_idx', Ref_Std_Day_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'ref_086_day_class_idx', Ref_086_Day_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'ref_ratio_day_class_idx', Ref_Ratio_Day_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'ref_138_day_class_idx', Ref_138_Day_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'ref_375_day_class_idx', Ref_375_Day_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'ndsi_day_class_idx', Ndsi_Day_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'btd_11_85_class_idx', Btd_11_85_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'btd_375_11_night_class_idx', Btd_375_11_Night_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'btd_375_11_day_class_idx', Btd_375_11_Day_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'btd_375_11_all_class_idx', Btd_375_11_All_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'bt_11_67_covar_class_idx', Bt_11_67_Covar_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'btd_11_67_class_idx', Btd_11_67_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'btd_85_73_class_idx', Btd_85_73_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'opd_063_day_class_idx', Opd_063_Day_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'log10_opd_063_day_class_idx', Log10_Opd_063_Day_Class_Idx ncdf_attput, cdfid, /GLOBAL, 'zcld_opd_class_idx', Zcld_Opa_Class_Idx end