;$Id:$ @apply_filters @fmft @write_2d_table @make_true_phase @smooth_nan ;-------------------------------------------------------------------------------------------- ; probabilistic phase for 2d naive bayesian ; ;-------------------------------------------------------------------------------------------- pro prob_overlap_2d, sensor_name = sensor_name, noplot=noplot, test_name=test_name, path_dir = path_dir if (~keyword_set(path_dir)) then path_dir='./' if (~keyword_set(sensor_name)) then sensor_name='myd02ssh' if (sensor_name eq 'myd02ssh') then training_source = 'ecm_training_myd02ssh_train.sav' if (sensor_name eq 'viirs') then training_source = 'ecm_training_viirs_train.sav' if (sensor_name eq 'abhi') then training_source = 'ecm_training_abhi_train.sav' if (sensor_name eq 'avhrr') then training_source = 'ecm_training_avhrr_train.sav' if (~keyword_set(test_name)) then test_name = "topa_beta1112" ;------------------------------- ; read in data ;------------------------------- training_path = './training_data/' restore, training_path + training_source f = training_data training_data = !null nb_rank = 2 nsfc = 1 mod_flag = 0 ;---- some parameters laplace_flag = 0 max_r = 100.0 count_min = 10 missing = -999.0 ;--- set any geographical limits lat_min = -90.0 lat_max = 90.0 ;--- calipso filtering cod_clear_ice_thresh = 0.1 cod_clear_water_thresh = 0.5 ;--- product filtering cod_cloud_thresh = 1.0 etrop_cloud_thresh = 0.9 ;--- default classifier limits - only modify in classifier definition if needed zen_min = 0.0 zen_max = 90.0 solzen_min = 0.0 solzen_max = 180.0 solglintzen_min = 0.0 solglintzen_max = 180.0 solscatang_min = 0.0 solscatang_max = 180.0 solglint_mask_min = 0 solglint_mask_max = 1 lunzen_min = 0.0 lunzen_max = 180.0 lunglintzen_min = 0.0 lunglintzen_max = 180.0 lunscatang_min = 0.0 lunscatang_max = 180.0 lunglint_mask_min = 0 lunglint_mask_max = 1 coast_mask_min = 0 coast_mask_max = 1 slant_tpw_min = 0.0 slant_tpw_max = 100.0 tsfc_min = 180.0 tsfc_max = 340.0 zsfc_min = -100.0 zsfc_max = 20000.0 zsfc_std_min = 0.0 zsfc_std_max = 1000.0 snow_class_min = 1 snow_class_max = 3 sfc_type_on_flag = make_array(nsfc,/INT,VALUE=1) ;----- set solar limits day_string = 'all' solzen_min = 0.0 solzen_max = 180.0 nc_header = '2D Overlap LUT' case test_name of "etropo_beta1185": begin classifier_name = 'etropo_beta1185' xtitle = 'etropo' x = f.emiss_tropo_11_0um_nom x_min = 0.0 & x_bin = 0.02 & nbins_x = 60 ytitle = 'beta1112' y = f.BETA_11UM_85UM_TROPOPAUSE y_min = 0.2 & y_bin = 0.02 & nbins_y = 60.0 nchan_used = 2 & channel_wvl_used = [11000,8500] end "topa_beta1112": begin classifier_name = 'topa_beta1112' xtitle = 'topa' x = f.cld_temp_opaque x_min = 200.0 & x_bin = 2.0 & nbins_x = 60 ytitle = 'beta1112' y = f.BETA_11UM_12UM_TROPOPAUSE y_min = 0.2 & y_bin = 0.02 & nbins_y = 60.0 nchan_used = 2 & channel_wvl_used = [11000,12000] end "etropo_beta1112": begin classifier_name = 'etropo_beta1112' xtitle = 'etropo' x = f.emiss_tropo_11_0um_nom x_min = 0.0 & x_bin = 0.02 & nbins_x = 60 ytitle = 'beta1112' y = f.BETA_11UM_12UM_TROPOPAUSE y_min = 0.2 & y_bin = 0.02 & nbins_y = 60.0 nchan_used = 2 & channel_wvl_used = [11000,12000] end "beta1112_beta11133": begin classifier_name = 'beta1112_beta11133' xtitle = 'beta1112' x = f.BETA_11UM_12UM_TROPOPAUSE x_min = 0.2 & x_bin = 0.02 & nbins_x = 60.0 ytitle = 'beta11133' y = f.BETA_11UM_133UM_TROPOPAUSE y_min = 0.0 & y_bin = 0.02 & nbins_y = 60.0 nchan_used = 3 & channel_wvl_used = [11000,12000,13300] end "topa_beta11133": begin classifier_name = 'topa_beta11133' xtitle = 'topa' x = f.cld_temp_opaque x_min = 200.0 & x_bin = 2.0 & nbins_x = 60 ytitle = 'beta11133' y = f.BETA_11UM_133UM_TROPOPAUSE y_min = 0.0 & y_bin = 0.02 & nbins_y = 60.0 nchan_used = 2 & channel_wvl_used = [11000,13300] end "etropo_beta11133": begin classifier_name = 'etropo_beta11133' xtitle = 'etropo' x = f.emiss_tropo_11_0um_nom x_min = 0.0 & x_bin = 0.02 & nbins_x = 60 ytitle = 'beta11133' y = f.BETA_11UM_133UM_TROPOPAUSE y_min = 0.0 & y_bin = 0.02 & nbins_y = 60.0 nchan_used = 2 & channel_wvl_used = [11000,13300] end "refl065_btd1112": begin classifier_name = 'cod065_btd1112' xtitle = 'btd1112' x = f.temp_11_0um_nom - f.temp_12_0um_nom x_min = -2 & x_bin = 0.25 & nbins_x = 40.0 ytitle = 'rel065' y = f.refl_0_65um_nom y_min = 0.0 & y_bin = 1.0 & nbins_y = 60.0 nchan_used = 2 & channel_wvl_used = [650,11000] solzen_min = 0.0 solzen_max = 80.0 end else: begin print, "error: unknown test" stop end endcase ;------------------------------------------------ ; fake other dims ;------------------------------------------------ ztitle = "blank" nbins_z = 1 z_min = 0.0 z_bin = 0.0 ;------------------------------------------------ ; make bin bounds ;------------------------------------------------ x_x = x_min + findgen(nbins_x)*x_bin + x_bin/2 x_y = y_min + findgen(nbins_y)*y_bin + y_bin/2 ;----------------------------------------------- ; make surface type ;---------------------------------------------- sfc_type_mask = f.bayes_mask_sfc_type ;----------------------------------------------- ; apply filters ;---------------------------------------------- phase_lidar = f.CLOSEST_CALIPSO_CLOUD_PHASE_TOP_LAYER cod_lidar = f.CLOSEST_CALIPSO_COD slant_tpw = f.total_precipitable_water_nwp / cos(f.sensor_zenith_angle * !DTOR) frac_lidar = f.CLOSEST_CALIPSO_CLOUD_FRACTION_15KM tc_lidar = f.CLOSEST_CALIPSO_TOP_MID_TEMPERATURE codi_lidar = f.CLOSEST_CALIPSO_COD_ICE overlap_lidar = f.CLOSEST_CALIPSO_MULTILAYER_FLAG idx_ice = where(phase_lidar eq 1 or phase_lidar eq 3, count_ice) idx_overlap = where(overlap_lidar eq 1, count_overlap) print, 'initial overlap fraction = ', float(count_overlap) / float(count_ice) ;------------------------------------------------------------------------------ ; set non-ice to missing since we only want to process ice phase ;------------------------------------------------------------------------------ make_true_phase_no_filter, phase_lidar, phase_true idx_ice = where(phase_true eq 2, count_ice) idx_overlap = where(overlap_lidar eq 1 and phase_true eq 2, count_overlap) print, 'second overlap fraction = ', float(count_overlap) / float(count_ice) make_true_phase, phase_lidar, frac_lidar, tc_lidar, cod_lidar, codi_lidar, $ cod_clear_water_thresh, cod_clear_ice_thresh, $ f.bayes_mask_sfc_type, f.cld_temp_opaque, f.cld_height_opaque, $ f.solar_zenith_angle, f.glint_mask, f.coast_mask, $ f.surface_elevation, f.surface_temperature_nwp, $ f.refl_0_65um_nom_min_3x3, f.refl_0_65um_nom, $ f.temp_11_0um_nom_max_3x3, f.temp_11_0um_nom, $ f.cld_opd_mask_0_65um_nom, cod_cloud_thresh, $ f.emiss_tropo_11_0um_nom, etrop_cloud_thresh, mod_flag, phase_true, missing idx_ice = where(phase_true eq 2, count_ice) idx_overlap = where(phase_true eq 2 and overlap_lidar eq 1, count_overlap) print, 'final overlap fraction = ', float(count_overlap) / float(count_ice) ;--- remove clouds where overlap is not possible ;idx = where(f.emiss_tropo_11_0um_nom ge 0.8, cc) ;if (cc gt 0) then phase_true[idx] = -1 ;idx = where(f.temp_11_0um_nom lt 250.0, cc) ;if (cc gt 0) then phase_true[idx] = -1 apply_filters, phase_true, $ f.sensor_zenith_angle, zen_min, zen_max, $ f.solar_zenith_angle, solzen_min, solzen_max, $ f.lunar_zenith_angle, lunzen_min, lunzen_max, $ f.land_class, $ f.glint_zenith_angle, solglintzen_min, solglintzen_max, $ f.lunar_glint_zenith_angle, lunglintzen_min, lunglintzen_max, $ f.scattering_angle, solscatang_min, solscatang_max, $ f.snow_class, snow_class_min, snow_class_max, $ f.glint_mask, solglint_mask_min, solglint_mask_max, $ f.coast_mask, coast_mask_min, coast_mask_max, $ slant_tpw, slant_tpw_min, slant_tpw_max, $ f.surface_temperature_nwp, tsfc_min, tsfc_max, $ f.surface_elevation, zsfc_min, zsfc_max, $ f.surface_elevation_stddev_3x3, zsfc_std_min, zsfc_std_max, $ f.latitude, lat_min, lat_max, missing ;------------------------------- ; select valid data within these limits ;------------------------------- ;--- remove all clear and water idx = where(phase_true eq 2 and overlap_lidar ge 0 and x ne missing and y ne missing,cc) print, 'total number after filtering= ', cc ;--- throw away filtered data x = x[idx] y = y[idx] phase_true = phase_true[idx] sfc_type_mask = sfc_type_mask[idx] overlap_lidar = overlap_lidar[idx] idx_save = idx ;------------------------------------------------------------- ; convert phase true to be limited to single and overlap ice ; anything that is not ice should be gone from here ;------------------------------------------------------------- temp = phase_true idx = where(phase_true ne 2,cc) if (cc gt 0) then temp[idx] = -1 idx = where(phase_true eq 2,cc) if (cc gt 0) then temp[idx] = 0 phase_true = temp idx = where(phase_true eq 0 and overlap_lidar eq 1,cc) if (cc gt 0) then phase_true[idx] = 1 ;features = make_array(3,phase_true.length) ;features[0,*] = f[idx_save].emiss_tropo_11_0um_nom ;features[1,*] = f[idx_save].BETA_11UM_12UM_TROPOPAUSE ;features[2,*] = f[idx_save].BETA_11UM_133UM_TROPOPAUSE ;labels = phase_true ;save, file = 'features.sav', features, labels ;stop ;------------------------------------------------------------- ; make arrays ;------------------------------------------------------------- single_table = make_array(nbins_x, nbins_y, /FLOAT, VALUE=MISSING) overlap_table = make_array(nbins_x, nbins_y, /FLOAT, VALUE=MISSING) obs_table = make_array(nbins_x, nbins_y, /FLOAT, VALUE=0.0) single_fraction = missing overlap_fraction = missing observation_count = missing ;-------------------------------------------------------------------- ; loop over surface type ;-------------------------------------------------------------------- x_sfc = x y_sfc = y phase_true_sfc = phase_true observation_count = phase_true_sfc.length ;--- make priors idx = where(phase_true_sfc ge 0, count_all) idx = where(phase_true_sfc eq 0, count_yes_single_all) idx = where(phase_true_sfc eq 1, count_yes_overlap_all) single_fraction = float(count_yes_single_all) / float(count_all) overlap_fraction = float(count_yes_overlap_all) / float(count_all) print, 'fractions = ', single_fraction, overlap_fraction cond_ratio_single = make_array(nbins_x, nbins_y, /FLOAT, VALUE=MISSING) cond_ratio_overlap = make_array(nbins_x, nbins_y, /FLOAT, VALUE=MISSING) post_prob_single = make_array(nbins_x, nbins_y, /FLOAT, VALUE=MISSING) post_prob_overlap = make_array(nbins_x, nbins_y, /FLOAT, VALUE=MISSING) obs_prob_lidar = make_array(nbins_x, nbins_y, /FLOAT, VALUE=0.0) post_prob_single_direct = make_array(nbins_x, nbins_y, /FLOAT, VALUE=MISSING) post_prob_overlap_direct = make_array(nbins_x, nbins_y, /FLOAT, VALUE=MISSING) for i = 0, nbins_x -1 do begin ;i point_loop x_min_bin = x_min + i*x_bin x_max_bin = x_min + (i+1)*x_bin for j = 0, nbins_y -1 do begin ;j-point loop y_min_bin = y_min + j*y_bin y_max_bin = y_min + (j+1)*y_bin idx_all = where(x_sfc ge x_min_bin and x_sfc lt x_max_bin and $ y_sfc ge y_min_bin and y_sfc lt y_max_bin and $ phase_true_sfc ge 0, count_all_bin) if (count_all_bin le 0) then goto, skip_bin if (count_all_bin lt count_min) then goto, skip_bin idx_single = where(phase_true_sfc[idx_all] eq 0, count_yes_single_bin) idx_overlap = where(phase_true_sfc[idx_all] eq 1, count_yes_overlap_bin) ;--- single count_yes_bin = count_yes_single_bin count_yes_all = count_yes_single_all count_no_bin = count_all_bin - count_yes_single_bin count_no_all = count_all - count_yes_single_all prior_yes = single_fraction compute_class_cond, nb_rank, laplace_flag, $ count_yes_all, count_no_all, count_yes_bin, count_no_bin, $ max_r, prior_yes, cond_yes, cond_no, r, posterior_prob cond_ratio_single[i,j] = r post_prob_single[i,j] = posterior_prob post_prob_single_direct[i,j] = float(count_yes_single_bin) / float(count_all_bin) ;--- overlap count_yes_bin = count_yes_overlap_bin count_yes_all = count_yes_overlap_all count_no_bin = count_all_bin - count_yes_overlap_bin count_no_all = count_all - count_yes_overlap_all prior_yes = overlap_fraction compute_class_cond, nb_rank, laplace_flag, $ count_yes_all, count_no_all, count_yes_bin, count_no_bin, $ max_r, prior_yes, cond_yes, cond_no, r, posterior_prob cond_ratio_overlap[i,j] = r post_prob_overlap[i,j] = posterior_prob post_prob_overlap_direct[i,j] = float(count_yes_overlap_bin) / float(count_all_bin) if (count_all gt count_min) then begin obs_prob_lidar[i,j] = float(count_all_bin) endif skip_bin: endfor ;end of j-point loop endfor ;end of i-point loop ;--- normalize obs_prob idx = where(obs_prob_lidar gt 0,cc) if (cc gt 0) then obs_prob_lidar[idx] = obs_prob_lidar[idx] / total(obs_prob_lidar) obs_prob_thresh = 3.0 / nbins_x / nbins_y smooth_width = 3 ;--- new smooth cond_ratio_single = smooth_nan(cond_ratio_single,smooth_width,missing) cond_ratio_overlap = smooth_nan(cond_ratio_overlap,smooth_width,missing) post_prob_single = smooth_nan(post_prob_single,smooth_width,missing) post_prob_overlap = smooth_nan(post_prob_overlap,smooth_width,missing) idx = where(obs_prob_lidar eq 0,cc) & print, 'number with zero = ', cc skip_post_process: ;------------------------------------------------------------------------------------ ; apply table ;------------------------------------------------------------------------------------ n = n_elements(phase_true_sfc) phase_retrieved = make_array(n, /FLOAT, VALUE = missing) for k = 0, n - 1 do begin i = min([nbins_x -1,max([0,fix((x_sfc[k] - x_min) / x_bin)])]) j = min([nbins_y -1,max([0,fix((y_sfc[k] - y_min) / y_bin)])]) ; max_prob = max([post_prob_clear[i,j], post_prob_water[i,j] , post_prob_ice[i,j]]) ; if (post_prob_clear[i,j] eq max_prob) then phase_retrieved[k] = 0 ; if (post_prob_water[i,j] eq max_prob) then phase_retrieved[k] = 1 ; if (post_prob_ice[i,j] eq max_prob) then phase_retrieved[k] = 2 phase_retrieved[k] = 0 if (post_prob_overlap[i,j] gt 0.50) then begin phase_retrieved[k] = 1 endif endfor idx = where(phase_true_sfc ge 0 and phase_retrieved ge 0, count_all) idx = where(phase_true_sfc eq 0 and phase_retrieved eq 0, count_agree_single) idx = where(phase_true_sfc eq 1 and phase_retrieved eq 1, count_agree_overlap) pod_overlap = float(count_agree_single + count_agree_overlap) / count_all print, 'pod overlap = ', pod_overlap ;------------------------------------------------------------------------------------ ; make images ;------------------------------------------------------------------------------------ if (~keyword_set(noplot)) then begin window, 0, xsize = 600, ysize = 600 load_color_table, 33, 0, 1.0, 249 erase, color = 0 view, post_prob_single, coltab=33, /current, position = [0.10,0.06,0.98,0.98], cbar_position=cbar_position, /nocbar, missing_color=1, $ sds_min=-101, sds_max=1, /mask_missing, title = 'SINGLE' plot, [0.0,0.0],[100.0,100.0],color = 3, background = 1 , $ xrange = [min(x_x),max(x_x)],ystyle = 1, /nodata, /noerase, $ yrange = [min(x_y),max(x_y)],xstyle = 1, $ xtitle = xtitle, ytitle = ytitle legend_local, ['SINGLE'], color = 3, /top, /right, charsize = 2, /clear ;saveimage, classifier_name+'_post_prob_single.png',/png window, 1, xsize = 600, ysize = 600 load_color_table, 33, 0, 1.0, 249 erase, color = 0 view, post_prob_overlap, coltab=33, /current, position = [0.10,0.06,0.98,0.98], cbar_position=cbar_position, /nocbar, missing_color=1, $ sds_min=-101, sds_max=1, /mask_missing, title = 'OVERLAP' plot, [0.0,0.0],[100.0,100.0],color = 3, background = 1 , $ xrange = [min(x_x),max(x_x)],ystyle = 1, /nodata, /noerase, $ yrange = [min(x_y),max(x_y)],xstyle = 1, $ xtitle = xtitle, ytitle = ytitle legend_local, ['OVERLAP'], color = 3, /top, /right, charsize = 2, /clear ;saveimage, classifier_name+'_post_prob_overlap.png',/png window, 2, xsize = 600, ysize = 600 load_color_table, 33, 0, 1.0, 249 erase, color = 0 view, obs_prob_lidar, coltab=33, /current, position = [0.10,0.06,0.98,0.98], cbar_position=cbar_position, /nocbar, missing_color=1, $ sds_min=-101, sds_max=max(obs_prob_lidar), /mask_missing, gam_fac=0.1 plot, [0.0,0.0],[100.0,100.0],color = 3, background = 1 , $ xrange = [min(x_x),max(x_x)],ystyle = 1, /nodata, /noerase, $ yrange = [min(x_y),max(x_y)],xstyle = 1, $ xtitle = xtitle, ytitle = ytitle legend_local, ['OBS_PROB'], color = 3, /top, /right, charsize = 2, /clear ;saveimage, classifier_name+'_obs_prob_lidar.png',/png window, 3, xsize = 600, ysize = 600 load_color_table, 33, 0, 1.0, 249 erase, color = 0 view, post_prob_overlap_direct, coltab=33, /current, position = [0.10,0.06,0.98,0.98], cbar_position=cbar_position, /nocbar, missing_color=1, $ sds_min=-101, sds_max=1, /mask_missing, title = 'OVERLAP DIRECT' plot, [0.0,0.0],[100.0,100.0],color = 3, background = 1 , $ xrange = [min(x_x),max(x_x)],ystyle = 1, /nodata, /noerase, $ yrange = [min(x_y),max(x_y)],xstyle = 1, $ xtitle = xtitle, ytitle = ytitle legend_local, ['OVERLAP'], color = 3, /top, /right, charsize = 2, /clear ;saveimage, classifier_name+'_post_prob_overlap_direct.png',/png endif ;-------------------------------------------- ; ;-------------------------------------------- obs_table[*,*] = obs_prob_lidar single_table[*,*] = cond_ratio_single overlap_table[*,*] = cond_ratio_overlap ;------------------------------------------------- ; make a structure to hold data and attributes ;------------------------------------------------- ;-- assign attributes @attributes_define.pro attributes.classifier_name = classifier_name attributes.length_classifier_name = strlen(classifier_name) attributes.description = nc_header attributes.sensor = sensor_name attributes.training_source = training_source attributes.timestamp = timestamp() attributes.cod_clear_ice_thresh = cod_clear_ice_thresh attributes.cod_clear_water_thresh = cod_clear_water_thresh attributes.rank = nb_rank attributes.nsfc = nsfc attributes.nchan_used = nchan_used attributes.x_name = xtitle attributes.nbins_x = nbins_x attributes.x_min = x_min attributes.x_bin = x_bin attributes.y_name = ytitle attributes.nbins_y = nbins_y attributes.y_min = y_min attributes.y_bin = y_bin attributes.z_name = ztitle attributes.nbins_z = nbins_z attributes.z_min = z_min attributes.z_bin = z_bin attributes.lat_min = lat_min attributes.lat_max = lat_max attributes.zen_min = zen_min attributes.zen_max = zen_max attributes.solzen_min = solzen_min attributes.solzen_max = solzen_max attributes.solglintzen_min = solglintzen_min attributes.solglintzen_max = solglintzen_max attributes.solscatang_min = solscatang_min attributes.solscatang_max = solscatang_max attributes.lunzen_min = lunzen_min attributes.lunzen_max = lunzen_max attributes.lunglintzen_min = lunglintzen_min attributes.lunglintzen_max = lunglintzen_max attributes.lunscatang_min = lunscatang_min attributes.lunscatang_max = lunscatang_max attributes.tpw_min = slant_tpw_min attributes.tpw_max = slant_tpw_max attributes.tsfc_min = tsfc_min attributes.tsfc_max = tsfc_max attributes.zsfc_min = zsfc_min attributes.zsfc_max = zsfc_max attributes.snow_class_min = snow_class_min attributes.snow_class_max = snow_class_max attributes.solglint_mask_min = solglint_mask_min attributes.solglint_mask_max = solglint_mask_max attributes.lunglint_mask_min = lunglint_mask_min attributes.lunglint_mask_max = lunglint_mask_max attributes.coast_mask_min = coast_mask_min attributes.coast_mask_max = coast_mask_max attributes.rut_solzen_thresh = -999.0 ;--- assign tables @lut_2d_overlap_define.pro lut.single = single_table lut.overlap = overlap_table lut.obs = obs_table lut.observation_count = observation_count lut.single_fraction = single_fraction lut.overlap_fraction = overlap_fraction lut.on_flag = sfc_type_on_flag channels_used.wvl = channel_wvl_used ;---------------------------------------------- ; write to netcdf file ;---------------------------------------------- nc_filename = path_dir + sensor_name + "_" + day_string+"_"+classifier_name+".nc" sd_id = ncdf_create(nc_filename,/CLOBBER,/NETCDF4_FORMAT) write_2d_overlap_table, sd_id, lut, attributes, channels_used, success_flag, /class_att ncdf_close, sd_id print, 'write success = ', success_flag end