$Id$ @read_1d_table @read_2d_table @read_3d_table @get_prob_1d @get_prob_2d @get_prob_3d @get_prior_1d @get_prior_2d @get_prior_3d @get_prior_map @compute_month_day @leapyr ;-------------------------------------------------------------------------------------- ; test IDL implentation of new naive bayes cloud mask/phase ;-------------------------------------------------------------------------------------- pro test_new_ecm, training_file = training_file, mask_in = mask_in, nskip = nskip, library_name = library_name use_mask_in_flag = 0 if (keyword_set(mask_in)) then use_mask_in_flag = 1 if (~keyword_set(nskip)) then nskip = 10 if (~keyword_set(library_name)) then library_name = "table_library.inc" table_path = "tables/" ;-------------------------------------------------------------------------------------- ;-- constants ;-------------------------------------------------------------------------------------- missing = -999.0 cloud_clear_prob_thresh = 0.5 prob_min_thresh = 0.001 r_max = 1000.0 cod_clear_thresh = 0.100000 cod_ice_water_thresh = 0.100000 lat_min = -90.0 lat_max = 90.0 etrop_cloud_thresh = 0.9 count_min = 10.0 nsfc = 7 solzen_min = 0.0 solzen_max = 180.0 glintzen_min = 0.0 glintzen_max = 180.0 zsfc_min = -100.0 zsfc_max = 20000.0 snow_class_min = 1 snow_class_max = 3 ;-------------------------------------------------------------------------------------- ;--- read prior ;-------------------------------------------------------------------------------------- prior_map_name = table_path + 'nb_cloud_mask_calipso_prior.nc' get_prior_map, prior_map_name ;-------------------------------------------------------------------------------------- ;--- read in observations ;-------------------------------------------------------------------------------------- restore, training_file nx = n_elements(lat) logcod065 = cod065 idx = where(cod065 gt 0.0,cc) & if (cc gt 0) then logcod065[idx] = alog10(cod065[idx]) logcod138 = cod138 idx = where(cod138 gt 0.0,cc) & if (cc gt 0) then logcod138[idx] = alog10(cod138[idx]) ;-------------------------------------------------------------------------------------- ;-- define sensor ;-------------------------------------------------------------------------------------- sensor = 'modis' ;-------------------------------------------------------------------------------------- ;-- define table ;-------------------------------------------------------------------------------------- @table_library_1d.inc ;-------------------------------------------------------------------------------------- ;--- read tables ;------------------------------------------------------------------------------------- table_read_flag[*] = 0 for itable = 0, ntables - 1 do begin temp_table_name = table_path + sensor + "_" + table_name[itable] + ".nc" if (table_rank[itable] eq 1) then begin read_1d_table, temp_table_name, lut_temp, attrs_temp, read_flag *(lut[itable]) = lut_temp *(attrs[itable]) = attrs_temp table_read_flag[itable] = read_flag endif if (table_rank[itable] eq 2) then begin read_2d_table, temp_table_name, lut_temp, attrs_temp, read_flag *(lut[itable]) = lut_temp *(attrs[itable]) = attrs_temp table_read_flag[itable] = read_flag endif if (table_rank[itable] eq 3) then begin read_3d_table, temp_table_name, lut_temp, attrs_temp, read_flag *(lut[itable]) = lut_temp *(attrs[itable]) = attrs_temp table_read_flag[itable] = read_flag endif ;--- turn off tables that were not read in if (table_read_flag[itable] eq 0) then table_on_flag[itable] = 0 endfor ;------------------------------------------------------------------------------------- ; check that the channels are on for each table ;------------------------------------------------------------------------------------- for itable = 0, ntables - 1 do begin for i = 0, nchan_tables - 1 do begin if(chan_on_table[itable,i] gt 0) then begin if(chan_on_flag[chan_on_table[itable,i]] eq 0) then table_on_flag[itable] = 0 endif endfor endfor ;------------------------------------------------------------------------------------- ; define output ;------------------------------------------------------------------------------------- post_prob_cloud = make_array(nx,/FLOAT,VALUE=missing) post_prob_clear = make_array(nx,/FLOAT,VALUE=missing) post_prob_ice = make_array(nx,/FLOAT,VALUE=missing) post_prob_water= make_array(nx,/FLOAT,VALUE=missing) r_clear = make_array(nx,/FLOAT,VALUE=missing) r_water = make_array(nx,/FLOAT,VALUE=missing) r_ice = make_array(nx,/FLOAT,VALUE=missing) cloud_mask_new = make_array(nx,/INT,VALUE=-1) cloud_phase_true = make_array(nx,/INT,VALUE=-1) cloud_phase_new = make_array(nx,/INT,VALUE=-1) cloud_phase_aux = make_array(nx,/INT,VALUE=-1) cloud_phase_old = make_array(nx,/INT,VALUE=-1) ;----------------------------------------------------------- ; determine first valid table index ;----------------------------------------------------------- first_valid_table = -1 itable = 0 while (first_valid_table lt 0 and itable lt ntables) do begin if (table_read_flag[itable] eq 1) then first_valid_table = itable itable++ endwhile if (first_valid_table lt 0) then begin print, 'ERROR, no valid tables found' endif ;----------------------------------------------------------- ; Loop over Pixels and generate probabilities ;----------------------------------------------------------- for i = 0, nx-1, nskip do begin isfc = bayes_mask_type[i] - 1 ; idl is zero based ;--- prior compute_month_day, fix(doy[i]), leapyr(fix(year[i])), month, day imonth = month - 1 idiurnal = 0 get_prior_values, lon[i], lat[i], imonth, idiurnal, prior_prob_clear, prior_prob_ice, prior_prob_water ;--- take prior for this from any of the luts ; f = *(lut[first_valid_table]) ; prior_prob_cloud = f[isfc].cloud_fraction ; prior_prob_clear = 1.0 - prior_prob_cloud ; prior_prob_ice = f[isfc].ice_fraction ; prior_prob_water = f[isfc].water_fraction ; ;-- base prior on etrop ; x = etrop[i,j] ; get_prior, x, isfc, etrop_btd8511_lut, etrop_btd8511_attrs, prior_prob_clear, prior_prob_ice, prior_prob_water ;--- force water and ice priors to be equal (combats over-icing) ; prior_prob_water = 0.5*(1.0-prior_prob_clear) ; prior_prob_ice = prior_prob_water ;---- reset table probabilities to missing clear_cond_ratio[*] = missing water_cond_ratio[*] = missing ice_cond_ratio[*] = missing obs_prob[*] = missing ;---- loop through each table for itable = 0, ntables - 1 do begin if (table_on_flag[itable] eq 0) then goto, skip_table x = missing y = missing z = missing case table_name[itable] OF 'all_etropo_topa_logbt11std': begin x = etrop[i] y = topa[i] z = logt11_std[i] end 'all_etropo': begin x = etrop[i] end 'all_etropo_logbt11std': begin x = etrop[i] y = logt11_std[i] end 'all_etropo_btd8511': begin x = etrop[i] y = t85[i] - t11[i] end 'all_etropo_btd1112': begin x = etrop[i] y = t11[i] - t12[i] end 'all_etropo_bt11bt67covari': begin x = etrop[i] y = bt11bt67covar[i] end 'day_etropo_emiss375': begin x = etrop[i] y = emiss375[i] end 'night_etropo_emiss375': begin x = etrop[i] y = emiss375[i] end 'all_btd1112_btd8511': begin x = t11[i] - t12[i] y = t85[i] - t11[i] end 'day_etropo_logcod065': begin x = etrop[i] y = cod065[i] if (y gt 0.0) then y = alog10(y) end 'day_etropo_logcod138': begin x = etrop[i] y = cod138[i] if (y gt 0.0) then y = alog10(y) end 'day_ref160_ref138': begin x = ref160[i] y = ref138[i] end 'day_etropo_ref160': begin x = etrop[i] y = ref160[i] end 'day_ref160_ref138': begin x = ref160[i] y = ref138[i] end endcase if (table_rank[itable] eq 1) then begin get_prob_1d, x, solzen[i], glintzen[i], zsfc[i], snow_class[i], isfc,$ *(lut[itable]), *(attrs[itable]), prob_min_thresh, r_max, missing, $ cloud_mask_in, use_mask_in_flag, $ z1, z2, z3, z4 endif if (table_rank[itable] eq 2) then begin get_prob_2d, x, y, solzen[i], glintzen[i], zsfc[i], snow_class[i], isfc,$ *(lut[itable]), *(attrs[itable]), prob_min_thresh, r_max, missing, $ cloud_mask_in, use_mask_in_flag, $ z1, z2, z3, z4 endif if (table_rank[itable] eq 3) then begin get_prob_3d, x, y, z, solzen[i], glintzen[i], zsfc[i], snow_class[i], isfc,$ *(lut[itable]), *(attrs[itable]), prob_min_thresh, r_max, missing, $ cloud_mask_in, use_mask_in_flag, $ z1, z2, z3, z4 endif ;----- clear_cond_ratio[itable] = z1 water_cond_ratio[itable] = z2 ice_cond_ratio[itable] = z3 obs_prob[itable] = z4 skip_table: endfor ;-------------------------------------------------------------------------------------------- ; combine probabilities from each tables into a final value using standard naive bayes ;-------------------------------------------------------------------------------------------- ;post_prob_clear[i] = prior_prob_clear ;post_prob_water[i] = prior_prob_water ;post_prob_ice[i] = prior_prob_ice r_clear[i] = 1.0 r_water[i] = 1.0 r_ice[i] = 1.0 for itable = 0, ntables - 1 do begin if (clear_cond_ratio[itable] ne missing) then r_clear[i] = r_clear[i] * clear_cond_ratio[itable] if (water_cond_ratio[itable] ne missing) then r_water[i] = r_water[i] * water_cond_ratio[itable] if (ice_cond_ratio[itable] ne missing) then r_ice[i] = r_ice[i] * ice_cond_ratio[itable] ;---> if (water_prob[itable] ne missing) then post_prob_water[i] = post_prob_water[i] * water_prob[itable] ;---> if (ice_prob[itable] ne missing) then post_prob_ice[i] = post_prob_ice[i] * ice_prob[itable] endfor ;--- turn ratios into probabilities via post_prob = 1.0 / ( 1.0 + r/prior_yes - r) post_prob_clear[i] = 1.0 / (1.0 + r_clear[i]/prior_prob_clear - r_clear[i]) post_prob_water[i] = 1.0 / (1.0 + r_water[i]/prior_prob_water - r_water[i]) post_prob_ice[i] = 1.0 / (1.0 + r_ice[i]/prior_prob_ice - r_ice[i]) ;--- normalize these probabilties so that they sum to 1 post_sum = post_prob_clear[i] + post_prob_water[i] + post_prob_ice[i] post_prob_clear[i] = post_prob_clear[i] / post_sum post_prob_water[i] = post_prob_water[i] / post_sum post_prob_ice[i] = post_prob_ice[i] / post_sum print, 'test = ', i, post_prob_clear[i], post_prob_water[i], post_prob_ice[i], post_prob_clear[i] + post_prob_water[i] + post_prob_ice[i] ;--- make a cloud probability from the complement of the clear probability post_prob_cloud[i] = 1.0 - post_prob_clear[i] ;--- make a 4-level cloud mask if (post_prob_cloud[i] ge 0.00) then cloud_mask_new[i] = 0 if (post_prob_cloud[i] ge 0.10) then cloud_mask_new[i] = 1 if (post_prob_cloud[i] ge 0.50) then cloud_mask_new[i] = 2 if (post_prob_cloud[i] ge 0.90) then cloud_mask_new[i] = 3 ;--- make a new cloud phase if (cloud_mask_new[i] eq 0 or cloud_mask_new[i] eq 1) then cloud_phase_new[i] = 0 if (cloud_mask_new[i] ge 2) then cloud_phase_new[i] = 1 if (cloud_mask_new[i] ge 2 and post_prob_ice[i] gt post_prob_water[i]) then cloud_phase_new[i] = 2 ;--- make old cloud phase if (cloud_type[i] eq 0 or cloud_type[i] eq 1) then cloud_phase_old[i] = 0 if (cloud_type[i] ge 2 and cloud_type[i] le 5) then cloud_phase_old[i] = 1 if (cloud_type[i] ge 6) then cloud_phase_old[i] = 2 ;--- make aux cloud phase if (cloud_type_aux[i] eq 0 or cloud_type_aux[i] eq 1) then cloud_phase_aux[i] = 0 if (cloud_type_aux[i] ge 2 and cloud_type_aux[i] le 5) then cloud_phase_aux[i] = 1 if (cloud_type_aux[i] ge 6) then cloud_phase_aux[i] = 2 ;--- adjust lidar if (cod_lidar[i] ne missing and cod_lidar[i] lt cod_clear_thresh) then phase_lidar[i] = 0 if (cod_lidar[i] gt cod_clear_thresh and codi_lidar[i] ne missing and codi_lidar[i] lt cod_ice_water_thresh) then phase_lidar[i] = 2 ;--- filter lidar (set bad to -1) if (frac_lidar[i] ne 0.0 and frac_lidar[i] ne 1.0) then phase_lidar[i] = -1 ;--- make a true cloud phase cloud_phase_true[i] = phase_lidar[i] if (phase_lidar[i] eq 2) then cloud_phase_true[i] = 1 if (phase_lidar[i] eq 1 or phase_lidar[i] eq 3) then cloud_phase_true[i] = 2 endfor ;----------------------------------------------------------------------- ; destroy pointers ;----------------------------------------------------------------------- lut = !null attrs = !null ;---------------------------------------------------------------------- ; performance ;--------------------------------------------------------------------- for isfc = 0, nsfc - 1 do begin idx_sfc = where(bayes_mask_type eq isfc + 1, cc) idx_valid = where(cloud_phase_true[idx_sfc] ge 0 and cloud_phase_new[idx_sfc] ge 0, count_valid) idx_agree_clear = where(cloud_phase_true[idx_sfc[idx_valid]] eq 0 and cloud_phase_new[idx_sfc[idx_valid]] eq 0, count_agree_clear) idx_agree_cloud = where(cloud_phase_true[idx_sfc[idx_valid]] ge 1 and cloud_phase_new[idx_sfc[idx_valid]] ge 1, count_agree_cloud) idx_agree_water = where(cloud_phase_true[idx_sfc[idx_valid]] eq 1 and cloud_phase_new[idx_sfc[idx_valid]] eq 1, count_agree_water) idx_agree_ice = where(cloud_phase_true[idx_sfc[idx_valid]] eq 2 and cloud_phase_new[idx_sfc[idx_valid]] eq 2, count_agree_ice) idx_cloud = where(cloud_phase_new[idx_sfc[idx_valid]] ge 1, count_cloud) pc_cloud_new = (float(count_agree_clear) + float(count_agree_cloud)) / float(count_valid) pc_phase_new = (float(count_agree_water) + float(count_agree_ice)) / float(count_cloud) print, "new = ", isfc, pc_cloud_new, pc_phase_new idx_valid = where(cloud_phase_true[idx_sfc] ge 0 and cloud_phase_new[idx_sfc] ge 0 and cloud_phase_old[idx_sfc] ge 0, count_valid) idx_agree_clear = where(cloud_phase_true[idx_sfc[idx_valid]] eq 0 and cloud_phase_old[idx_sfc[idx_valid]] eq 0, count_agree_clear) idx_agree_cloud = where(cloud_phase_true[idx_sfc[idx_valid]] ge 1 and cloud_phase_old[idx_sfc[idx_valid]] ge 1, count_agree_cloud) idx_agree_water = where(cloud_phase_true[idx_sfc[idx_valid]] eq 1 and cloud_phase_old[idx_sfc[idx_valid]] eq 1, count_agree_water) idx_agree_ice = where(cloud_phase_true[idx_sfc[idx_valid]] eq 2 and cloud_phase_old[idx_sfc[idx_valid]] eq 2, count_agree_ice) idx_cloud = where(cloud_phase_old[idx_sfc[idx_valid]] ge 1, count_cloud) pc_cloud_old = (float(count_agree_clear) + float(count_agree_cloud)) / float(count_valid) pc_phase_old = (float(count_agree_water) + float(count_agree_ice)) / float(count_cloud) print, "old = ", isfc, pc_cloud_old, pc_phase_old idx_valid = where(cloud_phase_true[idx_sfc] ge 0 and cloud_phase_new[idx_sfc] ge 0 and cloud_phase_aux[idx_sfc] ge 0, count_valid) idx_agree_clear = where(cloud_phase_true[idx_sfc[idx_valid]] eq 0 and cloud_phase_aux[idx_sfc[idx_valid]] eq 0, count_agree_clear) idx_agree_cloud = where(cloud_phase_true[idx_sfc[idx_valid]] ge 1 and cloud_phase_aux[idx_sfc[idx_valid]] ge 1, count_agree_cloud) idx_agree_water = where(cloud_phase_true[idx_sfc[idx_valid]] eq 1 and cloud_phase_aux[idx_sfc[idx_valid]] eq 1, count_agree_water) idx_agree_ice = where(cloud_phase_true[idx_sfc[idx_valid]] eq 2 and cloud_phase_aux[idx_sfc[idx_valid]] eq 2, count_agree_ice) idx_cloud = where(cloud_phase_aux[idx_sfc[idx_valid]] ge 1, count_cloud) pc_cloud_aux = (float(count_agree_clear) + float(count_agree_cloud)) / float(count_valid) pc_phase_aux = (float(count_agree_water) + float(count_agree_ice)) / float(count_cloud) print, "aux = ", isfc, pc_cloud_aux, pc_phase_aux print, "--------------" endfor end