;$Id$ @apply_filters @compute_class_cond ;-------------------------------------------------------------------------------------------- ; probabilistic phase for 1d naive bayesian ; ;-------------------------------------------------------------------------------------------- pro prob_sc_1d, topa=topa,etropo = etropo, logcod065 = logcod065, logcod138 = logcod138, ndsi = ndsi, $ refl160 = refl160, refl065 = refl065, dref_065 = dref_065, $ emiss_375 = emiss_375, btd8511 = btd8511, bt11bt67covari = bt11bt67covari, fmft = fmft, $ day = day, night=night ;------------------------------- ; read in data ;------------------------------- training_source = 'phase_training.sav' restore, training_source f = training_data training_data = !null f.zsfc = f.zsfc * 1000.0 ; km to m sensor_name='modis' ;------------------------------- ; set any geographical limits ;------------------------------- cod_clear_thresh = 0.100000 cod_ice_water_restore_thresh = 0.0 ;0.100000 cod_ice_water_remove_thresh = 0.01 lat_min = -90.0 lat_max = 90.0 etrop_cloud_thresh = 0.9 count_min = 10.0 missing = -999.0 ;----- set solar limits day_string = 'all' solzen_min = 0.0 solzen_max = 180.0 if (keyword_set(day)) then begin day_string = 'day' solzen_min = 0.0 solzen_max = 80.0 endif if (keyword_set(night)) then begin day_string = 'night' solzen_min = 90.0 solzen_max = 180.0 endif nsfc = 7 nc_header = '1D SC LUT' if (keyword_set(f.tc_opaque_cloud)) then begin classifier_name = 'f.tc_opaque_cloud' xtitle = 'Opaque Cloud Temperature' x = f.tc_opaque_cloud x_min = 180.0 & x_bin = 5.0 & nbins_x = 32 glintzen_min = 0.0 glintzen_max = 180.0 zsfc_min = -100.0 zsfc_max = 20000.0 snow_class_min = 1 snow_class_max = 3 endif if (keyword_set(etropo)) then begin classifier_name = 'etropo' xtitle = 'etrop' x = f.emiss_tropo x_min = -0.4 & x_bin = 0.01 & nbins_x = 150 glintzen_min = 0.0 glintzen_max = 180.0 zsfc_min = -100.0 zsfc_max = 20000.0 snow_class_min = 1 snow_class_max = 3 endif if (keyword_set(btd8511)) then begin classifier_name = 'btd8511' xtitle = 'btd8511' x = f.t85 - f.t11 x_min = -12.0 & x_bin = 0.2 & nbins_x = 100 glintzen_min = 0.0 glintzen_max = 180.0 zsfc_min = -100.0 zsfc_max = 20000.0 snow_class_min = 1 snow_class_max = 3 endif if (keyword_set(bt11bt67covari)) then begin classifier_name = 'bt11bt67covari' xtitle = 'bt11bt67covar' x = f.bt11bt67covar x_min = -3.0 & x_bin = 0.2 & nbins_x = 65 glintzen_min = 0.0 glintzen_max = 180.0 zsfc_min = -100.0 zsfc_max = 20000.0 snow_class_min = 1 snow_class_max = 3 endif if (keyword_set(dref_065)) then begin classifier_name = 'dref_065' xtitle = 'rgct' x = f.ref065 - f.ref065_clr x_min = -40.0 & x_bin = 1.0 & nbins_x = 120 glintzen_min = 40.0 glintzen_max = 180.0 zsfc_min = -100.0 zsfc_max = 20000.0 snow_class_min = 1 snow_class_max = 3 endif if (keyword_set(logcod065)) then begin classifier_name = 'logcod065' xtitle = 'logcod065' idx = where(cod065 gt 0.0,cc) x = f.cod065 x[*] = missing idx = where(f.cod065 gt 0) & if (cc gt 0) then x[idx] = alog10(cod065[idx]) x_min = -2 & x_bin = 0.1 & nbins_x = 45.0 glintzen_min = 40.0 glintzen_max = 180.0 zsfc_min = -100.0 zsfc_max = 20000.0 snow_class_min = 1 snow_class_max = 3 endif if (keyword_set(logcod138)) then begin classifier_name = 'logcod138' xtitle = 'logcod138' idx = where(f.cod138 gt 0.0,cc) x = f.cod138 x[*] = missing idx = where(f.cod138 gt 0) & if (cc gt 0) then x[idx] = alog10(cod138[idx]) x_min = -2 & x_bin = 0.1 & nbins_x = 40.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 endif if (keyword_set(emiss_375)) then begin classifier_name = 'emiss375' xtitle = 'emiss375' x = f.emiss375 x_min = 0.25 & x_bin = 0.05 & nbins_x = 70.0 glintzen_min = 40.0 glintzen_max = 180.0 zsfc_min = -100.0 zsfc_max = 20000.0 snow_class_min = 1 snow_class_max = 3 endif if (keyword_set(ndsi)) then begin classifier_name = 'ndsi' xtitle = 'ndsi' x = (f.ref065 - f.ref160) / (f.ref065 + f.ref160) idx = where(f.ref065 eq missing or f.ref160 eq missing,cc) if (cc gt 0) then x[idx] = missing x_min = -0.5 & x_bin = 0.025 & nbins_x = 75.0 glintzen_min = 40.0 glintzen_max = 180.0 zsfc_min = -100.0 zsfc_max = 20000.0 snow_class_min = 1 snow_class_max = 3 endif if (keyword_set(refl065)) then begin classifier_name = 'refl065' xtitle = 'refl065' x = f.ref065 x_min = -20 & x_bin = 2.0 & nbins_x = 80.0 glintzen_min = 40.0 glintzen_max = 180.0 zsfc_min = -100.0 zsfc_max = 20000.0 snow_class_min = 1 snow_class_max = 3 endif ;------------------------------------------------ ; make bin bounds ;------------------------------------------------ x_x = x_min + findgen(nbins_x)*x_bin + x_bin/2 ;----------------------------------------------- ; make surface type ;---------------------------------------------- sfc_type_mask = f.bayes_mask_sfc_type ;----------------------------------------------- ; apply filters ;---------------------------------------------- idx = where(f.phase_lidar eq 2,cc) print, 'before = ', cc apply_filters, f.phase_lidar, f.cld_frac_lidar, f.cod_lidar, f.codi_lidar, $ cod_clear_thresh, cod_ice_water_restore_thresh, cod_ice_water_remove_thresh, $ f.solzen, solzen_min, solzen_max, $ f.land_class, f.glintzen, glintzen_min, glintzen_max, $ f.snow_class, snow_class_min, snow_class_max, $ f.zsfc, zsfc_min, zsfc_max, $ f.lat, lat_min, lat_max, $ f.emiss_tropo, etrop_cloud_thresh, missing ;------------------------------- ; select valid data within these limits ;------------------------------- idx = where(f.phase_lidar ge 0.0 and x ne missing, cc) print, 'total number after filtering= ', cc ;--- throw away filtered data x = x[idx] phase_lidar = f[idx].phase_lidar cod_lidar = f[idx].cod_lidar tc_lidar = f[idx].tc_lidar zc_opa = f[idx].zc_opaque_cloud sfc_type_mask = sfc_type_mask[idx] nb_rank = 1 laplace_flag = 0 ;---------------------------------------------------------- ; change sc_mask_true = 0 = no, 1 = yes ;---------------------------------------------------------- sc_mask_true = phase_lidar ;---- prob relative to all pixels ;sc_mask_true[*] = -1 ;idx = where(phase_lidar eq 0 or phase_lidar eq 1 or phase_lidar eq 3 or tc_lidar gt 273.15, cc) ;if (cc gt 0) then sc_mask_true[idx] = 0 ;idx = where(phase_lidar eq 2 and tc_lidar gt 0.0 and tc_lidar le 273.15 and cod_lidar ge 1.0, cc) ;if (cc gt 0) then sc_mask_true[idx] = 1 ;-- a test ;idx = where(f.tc_opaque_cloud gt 250 and f.tc_opaque_cloud lt 273, cc) ;if (cc gt 0) then sc_mask_true[idx] = 1 ;---- prob relative to all clouds sc_mask_true[*] = -1 idx = where(phase_lidar eq 1 or phase_lidar eq 3 or tc_lidar gt 273.15 or tc_lidar lt 230.0, cc) if (cc gt 0) then sc_mask_true[idx] = 0 idx = where(phase_lidar eq 2 and tc_lidar gt 0.0 and tc_lidar le 273.15, cc) if (cc gt 0) then sc_mask_true[idx] = 1 ;------------------------------------------------------------- ; make arrays ;------------------------------------------------------------- cond_ratio_table = make_array(nbins_x, nsfc, /FLOAT, VALUE=MISSING) obs_prob_table = make_array(nbins_x, nsfc, /FLOAT, VALUE=0.0) sc_fraction = make_array(nsfc,/FLOAT,VALUE=missing) observation_count = make_array(nsfc,/FLOAT,VALUE=missing) ;-------------------------------------------------------------------- ; loop over surface type ;-------------------------------------------------------------------- ;for isfc = 0, nsfc - 1 do begin ;sfc_loop for isfc = 0, 0 do begin ;sfc_loop idx = where(sfc_type_mask eq isfc + 1 and sc_mask_true ne -1 ,cc) if (cc gt 0) then begin x_sfc = x[idx] sc_mask_true_sfc = sc_mask_true[idx] observation_count[isfc] = cc print, 'isfc, count = ', isfc, cc endif else begin print, 'skipping surface' goto, skip_sfc endelse ;--- make priors idx = where(sc_mask_true_sfc ge 0, count_all) idx = where(sc_mask_true_sfc eq 1, count_sc) sc_fraction[isfc] = float(count_sc) / float(count_all) prior_yes = sc_fraction[isfc] prior_no = 1.0 - prior_yes cond_ratio = make_array(nbins_x, /FLOAT, VALUE=MISSING) post_prob = make_array(nbins_x, /FLOAT, VALUE=MISSING) obs_prob = make_array(nbins_x, /FLOAT, VALUE=0.0) idx = where(sc_mask_true_sfc ge 0 and sc_mask_true_sfc eq 1,count_yes_all) idx = where(sc_mask_true_sfc ge 0 and sc_mask_true_sfc eq 0,count_no_all) cond_ratio[*] = 1.0 max_r = 100.0 for i = 0, nbins_x -1 do begin ;point_loop x_min_bin = x_min + i*x_bin x_max_bin = x_min + (i+1)*x_bin idx_all = where(x_sfc ge x_min_bin and x_sfc le x_max_bin and $ sc_mask_true_sfc ne -1, count_all_bin) idx_yes = where(sc_mask_true_sfc[idx_all] eq 1, count_yes_bin) idx_no = where(sc_mask_true_sfc[idx_all] eq 0, count_no_bin) compute_class_cond, nb_rank, laplace_flag, $ count_min, count_yes_all, count_no_all, count_yes_bin, count_no_bin, $ max_r, prior_yes, cond_yes, cond_no, r, posterior_prob cond_ratio[i] = r post_prob[i] = posterior_prob obs_prob[i] = float(count_all_bin) endfor ;end of pixel loop ;--- normalize obs_prob idx = where(obs_prob gt 0,cc) if (cc gt 0) then obs_prob[idx] = obs_prob[idx] / total(obs_prob) ;----------------- smooth for ix = 1, nbins_x - 2 do begin v = obs_prob[ix-1:ix+1] idx = where(v gt 0,cc) if (cc gt 0) then begin v = cond_ratio[ix-1:ix+1] cond_ratio[ix] = mean(v[idx]) v = post_prob[ix-1:ix+1] post_prob[ix] = mean(v[idx]) endif endfor ;------------------- ; fill in holes ;------------------- idx = where(obs_prob eq 0,cc) & print, 'number with zero = ', cc for i = 0, 5 do begin for ix = 1, nbins_x - 2 do begin if (obs_prob[ix] eq 0) then begin v = obs_prob[ix-1:ix+1] idx = where(v gt 0,cc) if (cc gt 0) then begin v = cond_ratio[ix-1:ix+1] cond_ratio[ix] = mean(v[idx]) v = post_prob[ix-1:ix+1] post_prob[ix] = mean(v[idx]) endif endif endfor endfor ;------------------------------------------------------------------------------------ ; apply table ;------------------------------------------------------------------------------------ n = n_elements(sc_mask_true_sfc) sc_retrieved = make_array(n, /FLOAT, VALUE = missing) for k = 0, n - 1 do begin i = min([nbins_x -1,max([0,fix((x_sfc[k] - x_min) / x_bin)])]) if (post_prob[i] ne missing) then sc_retrieved[k] = 0 if (post_prob[i] ge 0.50) then sc_retrieved[k] = 1 endfor idx = where(sc_mask_true_sfc ge 0 and sc_retrieved ge 0, count_all) idx = where(sc_mask_true_sfc eq 0 and sc_retrieved eq 0, count_agree_nosc) idx = where(sc_mask_true_sfc eq 0 and sc_retrieved eq 1, count_false_sc) idx = where(sc_mask_true_sfc eq 1 and sc_retrieved eq 0, count_missed_sc) idx = where(sc_mask_true_sfc eq 1 and sc_retrieved eq 1, count_agree_sc) idx = where(sc_mask_true_sfc eq 1, count_sc_mask_true) idx = where(sc_retrieved eq 1, count_sc_new) pod_sc = float(count_agree_sc + count_agree_nosc) / count_all missed_sc = float(count_missed_sc) / count_all false_sc = float(count_false_sc) / count_all frac_sc_true = float(count_sc_mask_true) / count_all frac_sc_new = float(count_sc_new) / count_all print, 'frac sc true = ', frac_sc_true print, 'frac sc new = ', frac_sc_new print, 'pod sc = ', pod_sc print, 'false sc = ', false_sc print, 'missed sc = ', missed_sc ;------------------------------------------------------------------------------------ ; make images ;------------------------------------------------------------------------------------ p0 = plot(x_x,post_prob, name = 'posterior_probability', line='solid', color='black',thick = 2,$ xtitle = xtitle, xrange = [min(x_x),max(x_x)], $ ytitle = 'frequency', yrange = [0.0, 1.0]) p1 = plot(x_x, post_prob, name = 'post_prob', line='solid', color='orange',thick = 2, /overplot) p2 = plot(x_x, obs_prob/max(obs_prob), name = 'obs_prob', line='solid', color='grey',thick = 2, /overplot) leg0 = legend(target=[p1,p2], position = [0.3,0.95], /normal, font_size = 8,/AUTO_TEXT_COLOR) t = text(0.3, 0.8, ['SFC'+string(isfc+1,format="(I02)")], color = "black",/data) obs_prob_table[*,isfc] = obs_prob cond_ratio_table[*,isfc] = cond_ratio skip_sfc: endfor ;sfc type loop ;---------------------------------------------- ; write to netcdf file ;---------------------------------------------- nc_filename = sensor_name + "_" + day_string+"_"+classifier_name+"_sc_prob.nc" sd_id = ncdf_create(nc_filename,/CLOBBER) ncdf_attput,sd_id, "description", nc_header, /CHAR, /global ncdf_attput,sd_id, "sensor", sensor_name, /CHAR, /global ncdf_attput,sd_id, "training_source", training_source, /CHAR, /global ncdf_attput,sd_id, "timestamp",timestamp(), /CHAR, /global ncdf_attput,sd_id, "nbins_x", float(nbins_x), /FLOAT, /global ncdf_attput,sd_id, "x_min", float(x_min), /FLOAT, /global ncdf_attput,sd_id, "x_bin", float(x_bin), /FLOAT, /global ncdf_attput,sd_id, "nsfc", float(nsfc), /FLOAT, /global ncdf_attput,sd_id, "latitude_minimum", lat_min, /FLOAT, /global ncdf_attput,sd_id, "latitude_maximum", lat_max, /FLOAT, /global ncdf_attput,sd_id, "calipso_cod_clear_threshold", cod_clear_thresh, /FLOAT, /global ncdf_attput,sd_id, "solar_zenith_minimum", solzen_min, /FLOAT, /global ncdf_attput,sd_id, "solar_zenith_maximum", solzen_max, /FLOAT, /global ncdf_attput,sd_id, "glint_zenith_minimum", glintzen_min, /FLOAT, /global ncdf_attput,sd_id, "glint_zenith_maximum", glintzen_max, /FLOAT, /global ncdf_attput,sd_id, "surface_elevation_minimum", zsfc_min, /FLOAT, /global ncdf_attput,sd_id, "surface_elevation_maximum", zsfc_max, /FLOAT, /global ncdf_attput,sd_id, "snow_class_minimum", snow_class_min, /FLOAT, /global ncdf_attput,sd_id, "snow_class_maximum", snow_class_max, /FLOAT, /global ;--- get dimensions nx_id = NCDF_DIMDEF(sd_id, 'nbins_x', nbins_x) nk_id = NCDF_DIMDEF(sd_id, 'nsfc', nsfc) ;--- create variable sds_id1 = NCDF_VARDEF(sd_id, 'cond_ratio_table', [nx_id,nk_id], /FLOAT) sds_id2 = NCDF_VARDEF(sd_id, 'obs_prob_table', [nx_id,nk_id], /FLOAT) sds_id3 = NCDF_VARDEF(sd_id, 'observation_count', [nk_id], /FLOAT) sds_id4 = NCDF_VARDEF(sd_id, 'prior_yes', [nk_id], /FLOAT) NCDF_CONTROL, sd_id, /ENDEF ; Put the file into data mode NCDF_VARPUT, sd_id, sds_id1, cond_ratio_table NCDF_VARPUT, sd_id, sds_id2, obs_prob_table NCDF_VARPUT, sd_id, sds_id3, observation_count NCDF_VARPUT, sd_id, sds_id4, sc_fraction ;--- close netcdf files ncdf_close, sd_id end