;$Id$ ;----------------------------------------------------------------------------- ; Routines to derive Bayesian Cloud Mask Classifier Files ; ; Requires idl save file from make_training_data.pro ; ; output is a text file with is used by CLAVR-x/PATMOS-x, figures and ; performance metrics ; ;----------------------------------------------------------------------------- @legend_local pro n_bayes_mask, time_diff_max = time_diff_max, SFC_TYPE_INPUT = sfc_type_input, $ SOLZEN_MIN=solzen_min, SOLZEN_MAX=solzen_max, AVHRR1=avhrr1, TRAINING_FILE=training_file, $ Emiss_Tropo = Emiss_Tropo, $ T_11 = T_11, $ T_Std = T_Std, $ T_Max = T_Max, $ 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, $ Ref_063_Day = Ref_063_Day, $ Ref_063_Min_Day = Ref_063_Min_Day, $ Ref_Std = Ref_Std, $ 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_All = Btd_11_85_All, $ 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, $ Vis_Opd_Day = Vis_Opd_Day, $ glint_filter = glint_filter, $ output_name = output_name, $ 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 if (keyword_set(output_name) eq 0) then output_name = 'temp' ;---- 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 Emiss_Tropo_Flag = keyword_set(Emiss_Tropo) T_11_Flag = keyword_set(T_11) T_Std_Flag = keyword_set(T_Std) T_Max_Flag = keyword_set(T_Max) FMFT_Flag = keyword_set(FMFT) Emiss_375_All_Flag = keyword_set(Emiss_375_All) Emiss_375_Day_Flag = keyword_set(Emiss_375_Day) Emiss_375_Night_Flag = keyword_set(Emiss_375_Night) Emiss_375_Emiss_Tropo_Night_Flag = keyword_set(Emiss_375_Emiss_Tropo_Night) Ref_Std_Flag = keyword_set(Ref_Std) Vis_Opd_Day_Flag = keyword_set(Vis_Opd_Day) Ref_063_Min_Day_Flag = keyword_set(Ref_063_Min_Day) Ref_063_Day_Flag = keyword_set(Ref_063_Day) Ref_086_Day_Flag = keyword_set(Ref_086_Day) Ref_Ratio_Day_Flag = keyword_set(Ref_Ratio_Day) Ref_138_Day_Flag = keyword_set(Ref_138_Day) Ref_160_Day_Flag = keyword_set(Ref_160_Day) Ref_375_Day_Flag = keyword_set(Ref_375_Day) Ndsi_Day_Flag = keyword_set(Ndsi_Day) Btd_11_85_Flag = keyword_set(Btd_11_85) Btd_375_11_Night_Flag = keyword_set(Btd_375_11_Night) Btd_375_11_Day_Flag = keyword_set(Btd_375_11_Day) Btd_375_11_All_Flag = keyword_set(Btd_375_11_All) Bt_11_67_Covar_Flag = keyword_set(Bt_11_67_Covar) Btd_11_67_Flag = keyword_set(Btd_11_67) Btd_85_73_Flag = keyword_set(Btd_85_73) n_class = Emiss_Tropo_Flag + T_11_Flag + T_Std_Flag + T_Max_Flag + $ FMFT_FLAG + Emiss_375_All_Flag + Emiss_375_Day_Flag + Emiss_375_Night_Flag + Emiss_375_Emiss_Tropo_Night_Flag + $ Vis_Opd_Day_Flag + Ref_063_Day_Flag + Ref_063_Min_Day_Flag + Ref_Std_Flag + $ Ref_086_Day_Flag + Ref_Ratio_Day_Flag + $ Ref_138_Day_Flag + Ndsi_Day_Flag + Ref_160_Day_Flag + Ref_375_Day_Flag + Btd_11_85_Flag + Btd_85_73_Flag + $ Bt_11_67_Covar_Flag + Btd_11_67_Flag + $ Btd_375_11_All_Flag + Btd_375_11_Day_Flag + Btd_375_11_Night_Flag print, "N_Class = ", n_class ;--------------------------------------------------------- ; set missing arguments to missing ;--------------------------------------------------------- if (keyword_set(TRAINING_FILE) eq 0) then begin ; training_file = 'goes12_list.sav' ; training_file = 'noaa18_list.sav' ; training_file = 'myd02ssh_list.sav' ; training_file = 'modis_list_replica.sav' ; training_file = '/Users/mfoster/Data/patmosx/bayes_training/modis_list_1km_2010_20.sav' training_file = '/home/heidinger/Projects/calipso_colocate/modis_1km_cf3.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_bounds_reg-1,n_class,/FLOAT,VALUE=missing) class_cond_no_reg = make_array(n_bounds_reg-1,n_class,/FLOAT,VALUE=missing) class_cond_ratio_reg = make_array(n_bounds_reg-1,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_bounds_reg-1,/FLOAT,VALUE=missing) class_cond_no_reg_temp = make_array(n_bounds_reg-1,/FLOAT,VALUE=missing) 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 ;--- convert time_diff_max to hours time_diff_max = time_diff_max/60.0 if (keyword_set(solzen_min) eq 0) then solzen_min = 0.0 if (keyword_set(solzen_max) eq 0) then solzen_max = 180.0 glintzen_max = 180.0 glintzen_min = 40.0 ;----------------------------------------------------------------------------- ; Read in the IDL save file ;----------------------------------------------------------------------------- data_file = training_file restore, data_file ;training_data = training_data_temp ;-----------------------------------------------------------------------------------------' ; make metrics needed for mask ;-----------------------------------------------------------------------------------------' ndata = n_elements(training_data.time_diff) missing_data = make_array(ndata,/FLOAT,VALUE=missing) ;--- check for NaN in median index = where(training_data.ref_065um_clr le 0, count) if (count gt 0) then training_data[index].ref_065um_clr = 5.0 t_ratio = (training_data.bt_11um_clr - training_data.bt_11um) / (training_data.bt_11um - 180.0) index = where(t_ratio gt 2, count) if (count gt 0) then t_ratio[index] = 2.0 ;----- make rfmft classifier rfmft = (training_data.bt_11um - training_data.bt_12um) - (training_data.btd_11um_12um_bt_11um_max_3x3) index = where(training_data.bt_11um eq missing or $ training_data.bt_12um eq missing or $ training_data.btd_11um_12um_bt_11um_max_3x3 eq missing or $ abs(training_data.bt_11um - training_data.bt_12um) 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 ;--- limit t_11um_Std index = where(training_data.bt_11um_std_3x3 gt 5.0, count) if (count gt 0) then training_data[index].bt_11um_std_3x3 = 5.0 ;----- make fmft classifier fmft = (training_data.bt_11um_clr - training_data.bt_12um_clr) * $ (training_data.bt_11um - 260.0)/(training_data.bt_11um_clr - 260.0) fmft[where(training_data.bt_11um_clr le 265)] = 0.0 fmft = (training_data.bt_11um - training_data.bt_12um) - fmft index = where(training_data.bt_11um eq missing or training_data.bt_12um eq missing or $ training_data.bt_11um_clr eq missing or training_data.bt_12um_clr eq missing or $ abs(training_data.bt_11um - training_data.bt_12um) gt 20.0, count) if (count gt 0) then fmft[index] = missing ;----- make fmft classifier btd_375_11_metric = (training_data.bt_375um - training_data.bt_11um) - $ (training_data.bt_375um_clr - training_data.bt_11um_clr) ;--- T11_max - T11 t_11um_rel = training_data.bt_11um_max_3x3 - training_data.bt_11um index = where(training_data.bt_11um 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 ;--- Ref_063_Day metric solzen_rgct_max = 80.0 rgct = training_data.ref_065um - training_data.ref_065um_clr index = where(training_data.ref_065um le 0.0 or training_data.ref_065um_clr le 0.0, count) if (count gt 0) then rgct[index] = missing rvct = training_data.ref_065um - training_data.ref_065um_min_3x3 index = where(training_data.ref_065um le 0.0 or training_data.ref_065um_min_3x3 le 0.0, count) if (count gt 0) then rvct[index] = missing ;--- reflectance ratio rrct = training_data.ref_086um / training_data.ref_065um index = where(training_data.ref_086um le 0.0 or training_data.ref_065um le 0.0, count) if (count gt 0) then rrct[index] = missing ;--- ndsi ndsi = (training_data.ref_065um - training_data.ref_160um) / (training_data.ref_065um + training_data.ref_160um) index = where(training_data.ref_160um le 0.0 or training_data.ref_065um le 0.0, count) if (count gt 0) then ndsi[index] = missing ;--- Ref_375_Day metric solzen_rgct_max = 80.0 ref_375_metric = training_data.ref_375um index = where(training_data.ref_375um le 0.0, count) if (count gt 0) then ref_375_metric[index] = missing ;----------------------------------------------------------------------------- ;--- make a surface type ;----------------------------------------------------------------------------- ;--- 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 ;---------------------------------------------------------------------- ; Determine Surface Type ;---------------------------------------------------------------------- ;--- comment out the following to not use CLAVR-x values in file sfc_type_mask = training_data.bayes_mask_sfc_type goto, skip_sfc_type_creation sfc_type_mask = make_array(n_elements(land_class),/BYTE,VALUE=0) print, 'STARTING sfc type definition' for i = 0L, n_elements(training_data.land_class) -1 do begin ;--- default LAND/OCEAN values if (training_data[i].land_class eq LAND) then begin sfc_type_mask[i] = 3 endif else begin sfc_type_mask[i] = 1 endelse ;--- #2 - Shallow Ocean if (training_data[i].land_class eq MODERATE_OCEAN or $ training_data[i].land_class eq DEEP_INLAND_WATER or $ training_data[i].land_class eq SHALLOW_INLAND_WATER or $ training_data[i].land_class eq SHALLOW_OCEAN) then begin sfc_type_mask[i] = 2 endif ;--- include into Shallow Ocean SST fronts if (training_data[i].land_class ne LAND and $ training_data[i].sst_back_uni gt 0.5) then begin sfc_type_mask[i] = 1 endif ;--- #3 Unfrozen Land if (training_data[i].land_class eq LAND or $ training_data[i].land_class eq COASTLINE or $ training_data[i].land_class eq EPHEMERAL_WATER) then begin sfc_type_mask[i] = 3 endif ;--- #4 Frozen Land if (training_data[i].snow_class eq SNOW and training_data[i].lat gt -60.0) then begin sfc_type_mask[i] = 4 endif ;--- #5 - Arctic if (training_data[i].lat gt 0 and training_data[i].snow_class eq SEA_ICE) then begin sfc_type_mask[i] = 5 endif ;---- #6 - Antarctica if (training_data[i].lat lt 0 and training_data[i].snow_class eq SEA_ICE) then begin sfc_type_mask[i] = 6 endif if (training_data[i].lat lt -60 and training_data[i].snow_class eq SNOW) then begin sfc_type_mask[i] = 6 endif if (training_data[i].lat gt 60 and training_data[i].lon gt -75.0 and training_data[i].lon lt -10 and $ (training_data[i].land_class eq LAND or training_data[i].land_class eq COASTLINE) and $ training_data[i].snow_class 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(training_data[i].lat) lt 60.0) then begin if(training_data[i].sfc_emiss_375um 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' 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(training_data.land_class) -1 do begin ;if (sfc_type_mask[i] eq 6) then begin ; print, training_data[i].lon, training_data[i].lat, training_data[i].cld_frac_lidar ; plots, training_data[i].lon, training_data[i].lat, psym = 3, color = 3 ;endif ;endfor ;stop ;============== count_total = n_elements(training_data.time_diff) 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 ;--- remove thin low cloud - BAD IDEA ; index = where(training_data.ec_lidar ne missing and training_data.Zc_lidar lt 4.0 and $ ; training_data.cld_frac eq 1.0 and training_data.ec_lidar lt 0.10, count) ; if (count gt 0) then training_data.cld_frac[index] = 0.0 ;----- make a binary calipso cld mask (0=clear, 3=cloud) cld_mask = training_data.cld_mask cld_mask_modis = training_data.cld_mask_atml2 cld_mask_calipso = cld_mask cld_mask_calipso[*] = missing index = where(training_data.cld_frac_lidar le cld_frac_clear_thresh and training_data.cld_frac ne missing, count) if (count gt 0) then cld_mask_calipso[index] = 0 index = where(training_data.cld_frac_lidar ge cld_frac_cloud_thresh and training_data.cld_frac ne missing, count) if (count gt 0) then cld_mask_calipso[index] = 3 ;--- force very thin cirrus to be called clear index = where(cld_mask_calipso eq 3 and $ training_data.ec_lidar ne missing and training_data.tc_lidar lt 250.0 and training_data.ec_lidar lt 0.03, count) if (count gt 0) then begin cld_mask_calipso[index] = 0 training_data[index].cld_frac_lidar = 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 $ training_data.emiss_tropo 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 $ training_data.ref_065um 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, solzen_max printf, lun, time_diff_max, n_class, n_bounds_reg, n_sfc_type cdfid = ncdf_create(output_name+'.nc',/CLOBBER) ncdf_attput, cdfid, /GLOBAL, 'data_file', data_file ncdf_attput, cdfid, /GLOBAL, 'domain_string', domain_string 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 ncdf_attput, cdfid, /GLOBAL, 'solzen_max', solzen_max 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_sfc_type', n_sfc_type endif ;------------------------------------------------------------------- ; loop through sfc type ;------------------------------------------------------------------ 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' ; if (sfc_type_input eq 8) then header = 'Mountains' ;------------------------------------------------------------------- ; select data for analysis ;------------------------------------------------------------------- index_valid = where(sfc_type_mask eq sfc_type_input and $ abs(training_data.time_diff) lt time_diff_max and $ training_data.lat le lat_north and $ training_data.lat ge lat_south and $ training_data.lon ge lon_west and $ training_data.lon le lon_east and $ training_data.tsfc_model ne missing and $ training_data.emiss_tropo ne missing and $ (rgct ne missing or training_data.solzen gt solzen_rgct_max) and $ training_data.ems_375um ne missing and training_data.ems_375um_clr ne missing and $ t_11um_rel ne missing and $ training_data.solzen ge solzen_min and training_data.solzen le solzen_max and $ training_data.glintzen ge glintzen_min and training_data.glintzen le glintzen_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) ;------------- compute selected classifiers iclass = 0 if (Emiss_Tropo_Flag eq 1) then begin X[*,iclass] = training_data[index_valid].emiss_tropo & X_name[iclass] = 'Emiss_Tropo' iclass = iclass + 1 endif if (Emiss_375_All_Flag eq 1) then begin X[*,iclass] = (training_data[index_valid].ems_375um-training_data[index_valid].ems_375um_clr) / $ training_data[index_valid].ems_375um_clr & X_name[iclass] = 'Emiss_375' iclass = iclass + 1 endif if (Emiss_375_Emiss_Tropo_Night_Flag eq 1) then begin X[*,iclass] = training_data[index_valid].ems_375um * training_data[index_valid].emiss_tropo & X_name[iclass] = 'Emiss_375_Emiss_Tropo_Night' iclass = iclass + 1 endif if (Emiss_375_Night_Flag eq 1) then begin X[*,iclass] = training_data[index_valid].ems_375um & 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' iclass = iclass + 1 endif if (Emiss_375_Day_Flag eq 1) then begin X[*,iclass] = (training_data[index_valid].ems_375um - training_data[index_valid].ems_375um_clr) / $ training_data[index_valid].ems_375um_clr & X_name[iclass] = 'Emiss_375_Day' iclass = iclass + 1 endif if (Ref_063_Day_Flag eq 1) then begin X[*,iclass] = rgct[index_valid] & X_name[iclass] = 'Ref_063_Day' iclass = iclass + 1 endif if (Ref_Std_Flag eq 1) then begin X[*,iclass] = training_data[index_valid].ref_065um_std_3x3 & X_name[iclass] = 'Ref_Std' iclass = iclass + 1 endif if (Ref_063_Min_Day_Flag eq 1) then begin X[*,iclass] = rvct[index_valid] & X_name[iclass] = 'Ref_063_Min_3x3_Day' iclass = iclass + 1 endif if (Ref_086_Day_Flag eq 1) then begin X[*,iclass] = training_data[index_valid].ref_086um & X_name[iclass] = 'Ref_086_Day' iclass = iclass + 1 endif if (Ref_Ratio_Day_Flag eq 1) then begin X[*,iclass] = rrct[index_valid] & X_name[iclass] = 'Ref_Ratio_Day' iclass = iclass + 1 endif if (Ref_138_Day_Flag eq 1) then begin X[*,iclass] = training_data[index_valid].ref_138um & X_name[iclass] = 'Ref_138_Day' ii = where(X[*,iclass] gt 30.0,count) if (count gt 0) then X[ii,iclass] = 30.0 iclass = iclass + 1 endif if (Ref_160_Day_Flag eq 1) then begin X[*,iclass] = training_data[index_valid].ref_160um & X_name[iclass] = 'Ref_160_Day' iclass = iclass + 1 endif if (Ref_375_Day_Flag eq 1) then begin X[*,iclass] = ref_375_metric[index_valid] & X_name[iclass] = 'Ref_375_Day' iclass = iclass + 1 endif if (Ndsi_Day_Flag eq 1) then begin X[*,iclass] = ndsi[index_valid] & X_name[iclass] = 'Ndsi_Day' iclass = iclass + 1 endif if (Vis_Opd_Day_Flag eq 1) then begin X[*,iclass] = training_data[index_valid].opd_mask & X_name[iclass] = 'Vis_Opd_Day' iclass = iclass + 1 endif if (T_11_Flag eq 1) then begin X[*,iclass] = training_data[index_valid].bt_11um & X_name[iclass] = 'T_11' iclass = iclass + 1 endif if (T_Std_Flag eq 1) then begin X[*,iclass] = training_data[index_valid].bt_11um_std_3x3 & X_name[iclass] = 'T_Std' iclass = iclass + 1 endif if (T_max_Flag eq 1) then begin X[*,iclass] = t_11um_rel[index_valid] & X_name[iclass] = 'T_Max-T' iclass = iclass + 1 endif if (FMFT_Flag eq 1) then begin X[*,iclass] = fmft[index_valid] & X_name[iclass] = 'FMFT' iclass = iclass + 1 endif if (Btd_11_85_Flag eq 1) then begin X[*,iclass] = training_data[index_valid].bt_11um - training_data[index_valid].bt_85um & X_name[iclass] = 'Btd_11_85' iclass = iclass + 1 endif if (Btd_85_73_Flag eq 1) then begin X[*,iclass] = training_data[index_valid].bt_85um - training_data[index_valid].bt_73um & X_name[iclass] = 'Btd_85_73' iclass = iclass + 1 endif if (Btd_375_11_All_Flag eq 1) then begin X[*,iclass] = btd_375_11_metric[index_valid] & X_name[iclass] = 'Btd_375_11_All' iclass = iclass + 1 endif if (Btd_375_11_Day_Flag eq 1) then begin X[*,iclass] = btd_375_11_metric[index_valid] & X_name[iclass] = 'Btd_375_11_Day' iclass = iclass + 1 endif if (Btd_375_11_Night_Flag eq 1) then begin X[*,iclass] = btd_375_11_metric[index_valid] & X_name[iclass] = 'Btd_375_11_Night' iclass = iclass + 1 endif if (Bt_11_67_Covar_Flag eq 1) then begin X[*,iclass] = training_data[index_valid].bt_11um_bt_67um_covar & X_name[iclass] = 'Bt_11_67_Covar' iclass = iclass + 1 endif if (Btd_11_67_Flag eq 1) then begin X[*,iclass] = training_data[index_valid].bt_11um - training_data[index_valid].bt_67um & X_name[iclass] = 'Btd_11_67' iclass = iclass + 1 endif ;--- put limits on Emiss_Tropo for iclass = 0, n_class -1 do begin if (X_name[iclass] eq 'Emiss_Tropo') then begin ii = where(training_data[index_valid].tsfc_model lt 230, cc) print, 'emiss tropo removed = ', cc if (cc gt 0) then X[ii,iclass] = missing endif endfor ;--- put limits on Emiss_375_Night for iclass = 0, n_class -1 do begin if (X_name[iclass] eq 'Emiss_375_Night') then begin ii = where(training_data[index_valid].solzen lt 90, cc) if (cc gt 0) then X[ii,iclass] = missing endif if (X_name[iclass] eq 'Btd_375_11_Night') then begin ii = where(training_data[index_valid].solzen lt 90, cc) if (cc gt 0) then X[ii,iclass] = missing endif endfor ;--- put limits on Emiss_375_Day for iclass = 0, n_class -1 do begin if (X_name[iclass] eq 'Emiss_375_Day') then begin ii = where(training_data[index_valid].solzen gt 90, cc) if (cc gt 0) then X[ii,iclass] = missing endif endfor ;--- put limits on Reflectance Tests for iclass = 0, n_class -1 do begin if (X_name[iclass] eq 'Ref_063_Day' or $ X_name[iclass] eq 'Ref_063_Min_3x3_Day' or $ X_name[iclass] eq 'Ref_Std' or $ X_name[iclass] eq 'Ref_086_Day' or $ X_name[iclass] eq 'Ref_Ratio_Day' or $ X_name[iclass] eq 'Ref_138_Day' or $ X_name[iclass] eq 'Ndsi_Day' or $ X_name[iclass] eq 'Ref_160_Day' or $ X_name[iclass] eq 'Btd_375_11_Day' or $ X_name[iclass] eq 'Ref_375_Day' $ ) then begin ii = where(training_data[index_valid].solzen gt solzen_rgct_max, cc) if (cc gt 0) then X[ii,iclass] = missing endif 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 data_cdf = X[ii,iclass] nbins_cdf = 1000 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 goto, skip_cdf 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 goto, skip_cdf endelse ii = where(X[*,iclass] ne missing and X[*,iclass] gt data_cdf_min and X[*,iclass] lt data_cdf_max,cc) print, iclass, cc X_bounds[*,iclass] = create_bin(X[ii,iclass],n_bins,missing) ; determine 1% and 99% values data_cdf = X[ii,iclass] nbins_cdf = 1000 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 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 endelse 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 $ this_class = $ where(X[yes_[ii],iclass] ge X_bounds[j,iclass] and X[yes_[ii],iclass] le X_bounds[j+1,iclass],nC) 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 $ this_class = $ where(X[no_[ii],iclass] ge X_bounds[j,iclass] and X[no_[ii],iclass] le X_bounds[j+1,iclass],nC) 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] ;------------------------------------------------------------------------------------- ; ;------------------------------------------------------------------------------------- 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 ;------------------------------------------------------------------------------------ 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) ; Posterior_by_class_bin[ibin,iclass] = (prior_yes[sfc_type_input]*cond_yes_bin) / $ ; (prior_yes[sfc_type_input]*cond_yes_bin + prior_no[sfc_type_input]*cond_no_bin) endfor 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 if (X_Name[iclass] eq 'Ref_063_Day') then begin savename = 'ref_063_day_'+strtrim(sfc_type_input,2)+'.sav' x = X_bin_centers_reg[*,iclass] y = Posterior_by_class_bin[*,iclass] save, file = savename, x, y endif if (X_Name[iclass] eq 'T_11') then begin savename = 't_11_'+strtrim(sfc_type_input,2)+'.sav' x = X_bin_centers_reg[*,iclass] y = Posterior_by_class_bin[*,iclass] save, file = savename, x, y endif 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 = training_data[index_valid].cld_frac 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 ;--- 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) ;--- 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) ;--- 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 endif NCDF_VARPUT, cdfid, vid3, prior_yes[sfc_type_input], OFFSET=[sfc_type_input-1] NCDF_VARPUT, cdfid, vid4, prior_no[sfc_type_input] , OFFSET=[sfc_type_input-1] NCDF_VARPUT, cdfid, vid5, optimal_posterior_prob[sfc_type_input] , OFFSET=[sfc_type_input-1] for iclass = 0, n_class -1 do begin NCDF_VARPUT, cdfid, vid6, first_valid_bounds[iclass]+1, OFFSET=[iclass,sfc_type_input-1] NCDF_VARPUT, cdfid, vid7, last_valid_bounds[iclass]+1, OFFSET=[iclass,sfc_type_input-1] NCDF_VARPUT, cdfid, vid8, bin_start[iclass], OFFSET=[iclass,sfc_type_input-1] NCDF_VARPUT, cdfid, vid9, bin_end[iclass], OFFSET=[iclass,sfc_type_input-1] NCDF_VARPUT, cdfid, vid10, delta_bin[iclass], OFFSET=[iclass,sfc_type_input-1] for ibin = 0,n_bins_reg-1 do begin NCDF_VARPUT, cdfid, vid11, class_cond_yes_reg[ibin,iclass], OFFSET=[ibin,iclass,sfc_type_input-1] NCDF_VARPUT, cdfid, vid12, class_cond_no_reg[ibin,iclass], OFFSET=[ibin,iclass,sfc_type_input-1] NCDF_VARPUT, cdfid, vid13, class_cond_ratio_reg[ibin,iclass], OFFSET=[ibin,iclass,sfc_type_input-1] endfor endfor ;------------------------------------------------- ; 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 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 ;X_bounds_reg = fltarr(n_bounds_reg) ;class_cond_yes_reg = fltarr(n_bounds_reg-1) ;class_cond_no_reg = fltarr(n_bounds_reg-1) 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