;"$Id:$" ;--------------------------------------------------------------- @read_group_attributes @read_1d_group_lut @read_2d_group_lut @read_3d_group_lut @read_prior @compute_prior_from_map @read_level2 @tut_routine @rut_routine @get_prob_mask_phase @grab_prior_from_lut @view pro test_on_level2, use_log = use_log use_log_flag = 0 if (~keyword_set(use_log)) then use_log_flag = 1 ;------------------------------------------------------------------------------ ; user input ;------------------------------------------------------------------------------ nc_filename= './combined_tables/test.nc' L2_Path = './level2_test/' xstride = 8 ystride = 8 ;acmf_path = './acmf/' acmf_path = 'none' acmf_file = 'none' Prior_Source = 1 ; 0 = lut, 1 = modis prior ;Prior_Name = 'prior/nb_cloud_mask_modis_prior.nc' Prior_Name = 'prior/nb_cloud_mask_calipso_prior.nc' L2_Names = file_search(L2_Path,'*level2.*') Num_L2 = L2_Names.length Nsfc = 7 Missing = -999.0 Prob_Min_Thresh = 0.001 R_Prod_Min = 0.0 R_Prod_Max = 1000.0 R_Log_Min = -10.0 R_Log_Max = 10.0 ;------------------------------------------------------------------------------ ; set channel on/off and channel mapping ; Chan_On = channel on / off binary flags ; Chan_Wvl = channel nominal wavelength in nm ; This is set up using CLAVR-x channel numbers ;------------------------------------------------------------------------------ Nchan = 45 Chan_On = make_array(Nchan,/INT,VALUE=1) Chan_On[31-1] = 1 Chan_Wvl = make_array(Nchan,/INT,VALUE=0) Chan_On[1-1] = 1 & Chan_Wvl[1-1] = 650 Chan_On[2-1] = 1 & Chan_Wvl[2-1] = 860 Chan_On[26-1] = 1 & Chan_Wvl[26-1] = 1380 Chan_On[6-1] = 1 & Chan_Wvl[6-1] = 1600 Chan_On[20-1] = 1 & Chan_Wvl[20-1] = 3750 Chan_On[27-1] = 1 & Chan_Wvl[27-1] = 6700 Chan_On[28-1] = 1 & Chan_Wvl[28-1] = 7300 Chan_On[29-1] = 1 & Chan_Wvl[29-1] = 8500 Chan_On[30-1] = 1 & Chan_Wvl[30-1] = 9700 Chan_On[31-1] = 1 & Chan_Wvl[31-1] = 11000 Chan_On[32-1] = 1 & Chan_Wvl[32-1] = 12000 Chan_On[33-1] = 1 & Chan_Wvl[33-1] = 13300 Chan_On[37-1] = 1 & Chan_Wvl[37-1] = 6200 Chan_On[38-1] = 1 & Chan_Wvl[38-1] = 10400 Chan_On[44-1] = 1 & Chan_Wvl[44-1] = 700 ;DNB ;------------------------------------------------------------------- ; Loop over level2 files ;------------------------------------------------------------------- for L2_Idx = 0, Num_L2-1 do begin L2_Name = L2_Names[L2_Idx] Root_Name = strmid(L2_Name,strpos(L2_Name,'clavrx_')) ;--- determine format netcdf_flag = 0 if (strpos(Root_Name,'.nc') gt 0) then netcdf_flag = 1 if (netcdf_flag eq 0) then Root_Name = strmid(Root_Name, 0, strpos(Root_Name,'.hdf')) if (netcdf_flag eq 1) then Root_Name = strmid(Root_Name, 0, strpos(Root_Name,'.nc')) Output_Name = Root_Name + '.sav' ;--- read acmf if appropriate acmf = missing if (acmf_path ne 'none') then begin temp = strsplit(Root_Name,'_', /extract) temp = temp[4] temp = strsplit(temp,'.',/extract) time_string = temp[0] acmf_file = file_search(acmf_path,'*'+time_string+'*.nc') acmf = read_standard(acmf_file[0],'BCM') acmf = congrid(acmf,542,542) endif ;------------------------------------------------------------------- ; read lookup table ;------------------------------------------------------------------- ;------------------------------------------------------------------ ; note, this does not follow the fortran code ; fortran will have things like lut[0].clear and attributes[0].nbins_x ; I continue to try and get idl to do this ;------------------------------------------------------------------ ;--- open file for output sd_id = ncdf_open(nc_filename) ;--- open file and get the number and names of classifiers ncdf_attget,sd_id, "nclassifiers", nclass, /global ;-- open lut file sds_id = ncdf_varid(sd_id,'classifier_names') ncdf_varget, sd_id, sds_id, classifier_names classifier_names = string(classifier_names) sds_id = ncdf_varid(sd_id, 'conf_clear_prob_clear_thresh') ncdf_varget, sd_id, sds_id, Conf_Clear_Prob_Clear_Thresh sds_id = ncdf_varid(sd_id, 'prob_clear_prob_cloud_thresh') ncdf_varget, sd_id, sds_id, Prob_Clear_Prob_Cloud_Thresh sds_id = ncdf_varid(sd_id, 'prob_cloud_conf_cloud_thresh') ncdf_varget, sd_id, sds_id, Prob_Cloud_Conf_Cloud_Thresh sds_id = ncdf_varid(sd_id, 'rut_clear_prob_clear_thresh') ncdf_varget, sd_id, sds_id, Rut_Clear_Prob_Clear_Thresh sds_id = ncdf_varid(sd_id, 'tut_clear_prob_clear_thresh') ncdf_varget, sd_id, sds_id, Tut_Clear_Prob_Clear_Thresh ;--- define space luts=ptrarr(nclass, /allocate_heap) attrs=ptrarr(nclass, /allocate_heap) channels_used=ptrarr(nclass, /allocate_heap) ;--- get group id numbers - each group is a classifier grp_ids = ncdf_groupsinq(sd_id) Class_Use_Flag = make_array(Nclass,/INT,VALUE=0) Class_Rank = make_array(Nclass,/INT,VALUE=0) Class_Xname = make_array(Nclass,/STRING,VALUE=0) Class_Yname = make_array(Nclass,/STRING,VALUE=0) Class_Zname = make_array(Nclass,/STRING,VALUE=0) for iclass = 0, Nclass -1 do begin grp_id = grp_ids[iclass] classifier_name = ncdf_groupname(grp_id) ;--- read attributes read_group_attributes, grp_id, attrs_temp *(attrs[iclass]) = attrs_temp case attrs_temp.rank of 1: read_1d_group_lut, grp_id, lut_temp, attrs_temp, chans_temp, success_flag, use_log_flag 2: read_2d_group_lut, grp_id, lut_temp, attrs_temp, chans_temp, success_flag, use_log_flag 3: read_3d_group_lut, grp_id, lut_temp, attrs_temp, chans_temp, success_flag, use_log_flag else: begin print, 'unknown rank of classifier in read_combined_table' success_flag = 0 end endcase *(luts[iclass]) = lut_temp *(channels_used[iclass]) = chans_temp Class_Use_Flag[iclass] = success_flag Class_Rank[iclass] = attrs_temp.rank Class_Xname[iclass] = attrs_temp.x_name Class_Yname[iclass] = attrs_temp.y_name Class_Zname[iclass] = attrs_temp.z_name endfor ;--- close file lut file ncdf_close, sd_id print, '==================================================' print, 'Finished Lut Read' print, '==================================================' ;------------------------------------------------------------------- ; turn off tables based on chan_on ;------------------------------------------------------------------- for Class_Idx = 0, Nclass -1 do begin if (Class_Use_Flag[Class_Idx] eq 0) then goto, skip_class_0 x = *(Channels_Used[Class_Idx]) for Chan_Idx = 0, Nchan - 1 do begin if (Chan_On[Chan_Idx] eq 0) then begin idx = where(Chan_Wvl[Chan_Idx] eq X.Wvl and Chan_Wvl[Chan_Idx] ne 0,cc) if (cc gt 0) then begin Class_Use_Flag[Class_Idx] = 0 endif endif endfor skip_class_0: print, Classifier_Names[Class_Idx], Class_Use_Flag[Class_Idx] endfor Number_Valid_Classes = total(Class_Use_Flag) if (Number_Valid_Classes eq 0) then begin print, 'no active tests' stop endif ;------------------------------------------------------------------- ; read prior table ;------------------------------------------------------------------- read_prior, Prior_Name, Prior ;--------------------------------------------------------------------------------------------- ; read level2 ;--------------------------------------------------------------------------------------------- print, 'Begin reading Level2' read_level2, L2_Name, L2, xstride, ystride print, 'Done reading Level2' ;--------------------------------------------------------------------------------------------- ; compute prior from a map ;--------------------------------------------------------------------------------------------- Diurnal = 1 compute_prior_from_map_2d, l2.Lon, l2.Lat, l2.Bad_Pixel_Mask, l2.Nx, l2.Ny, Missing, $ prior.lon_min, prior.dLon, prior.Nlon, $ prior.lat_min, prior.dLat, prior.Nlat, $ Prior, L2.Month, Diurnal, Pixel_Prior_Prob_Cloud, $ Pixel_Prior_Prob_Ice_Cloud, Pixel_Prior_Prob_Water_Cloud ;--------------------------------------------------------------------------------------------- ; compute prior from the values in the lut (the first one) ;--------------------------------------------------------------------------------------------- Sfc_Prior_Prob_Clear = make_array(Nsfc,/FLOAT,Value=Missing) Sfc_Prior_Prob_Water = make_array(Nsfc,/FLOAT,Value=Missing) Sfc_Prior_Prob_Ice = make_array(Nsfc,/FLOAT,Value=Missing) for Sfc_Idx = 1,Nsfc do begin grab_prior_from_lut, *(luts[0]), Sfc_Idx, Prior_Prob_Clear, Prior_Prob_Water, Prior_Prob_Ice Sfc_Prior_Prob_Clear[Sfc_Idx-1] = Prior_Prob_Clear Sfc_Prior_Prob_Water[Sfc_Idx-1] = Prior_Prob_Water Sfc_Prior_Prob_Ice[Sfc_Idx-1] = Prior_Prob_Ice endfor print, "After Compute Prior" ;--------------------------------------------------------------------------------------------- ; run data from ecm ;--------------------------------------------------------------------------------------------- Clear_Cond_Ratio = make_array(Nclass,/FLOAT,Value=Missing) Water_Cond_Ratio = make_array(Nclass,/FLOAT,Value=Missing) Ice_Cond_Ratio = make_array(Nclass,/FLOAT,Value=Missing) Obs_Prob = make_array(Nclass,/FLOAT,Value=Missing) Post_Prob_Cloud = make_array(L2.Nx,L2.Ny,/Float,Value=Missing) Post_Prob_Cloud_By_Class = make_array(L2.Nx,L2.Ny,Nclass,/Float,Value=Missing) Post_Prob_Clear = make_array(L2.Nx,L2.Ny,/Float,Value=Missing) Post_Prob_Water = make_array(L2.Nx,L2.Ny,/Float,Value=Missing) Post_Prob_Ice = make_array(L2.Nx,L2.Ny,/Float,Value=Missing) Cloud_Phase = make_array(L2.Nx,L2.Ny,/Integer,Value=Missing) Cloud_Phase_Uncer = make_array(L2.Nx,L2.Ny,/Float,Value=Missing) Cloud_Mask_Binary = make_array(L2.Nx,L2.Ny,/INT,VALUE=-1) Cloud_Mask = make_array(L2.Nx,L2.Ny,/INT,VALUE=-1) Cloud_Type = make_array(L2.Nx,L2.Ny,/Integer,Value=Missing) TUT = make_array(L2.Nx,L2.Ny,/INT,VALUE=-1) RUT = make_array(L2.Nx,L2.Ny,/INT,VALUE=-1) for Elem_Idx = 0, L2.Nx - 1 do begin print, 'processing Elem_Idx = ', Elem_Idx, ' of ', L2.Nx-1 for Line_Idx = 0, L2.Ny - 1 do begin ;--- check for missing data like space if (L2.zen[Elem_Idx,Line_Idx] lt 0.0) then goto, skip_pixel ;--- select surface type index Sfc_Idx = fix(l2.sfc_type_mask[Elem_Idx,Line_Idx]) if (Sfc_Idx lt 1) then goto, skip_pixel Sfc_Idx = Sfc_Idx - 1 ; idl is zero based ;--- TEST FILTER ; if (Sfc_Idx ne 0) then goto, skip_pixel ; if (L2.etropo11[Elem_Idx,Line_Idx] gt 0.2) then goto, skip_pixel ; if (L2.ndsi[Elem_Idx,Line_Idx] gt 0.2) then goto, skip_pixel ; if (L2.drefl065clr[Elem_Idx,Line_Idx] gt 10.0) then goto, skip_pixel ;--- TEST FILTER ;---- reset table probabilities to missing Clear_Cond_Ratio[*] = Missing Water_Cond_Ratio[*] = Missing Ice_Cond_Ratio[*] = Missing Obs_Prob[*] = Missing R_Clear = 1.0 R_Water = 1.0 R_Ice = 1.0 for Class_Idx = 0, Nclass - 1 do begin if (Class_Use_Flag[Class_Idx] eq 0) then goto, skip_class ;--- set up input x = Missing y = Missing z = Missing for i = 0, Class_Rank[Class_Idx] -1 do begin case i of 0: dim_name = Class_Xname[Class_Idx] 1: dim_name = Class_Yname[Class_Idx] 2: dim_name = Class_Zname[Class_Idx] endcase case dim_name of 'etropo11': value = L2.etropo11[Elem_Idx,Line_Idx] 'etropo10': value = L2.etropo10[Elem_Idx,Line_Idx] 'topa': value = L2.Topa[Elem_Idx,Line_Idx] 'zopa': value = L2.Zopa[Elem_Idx,Line_Idx] 'tsfc': value = L2.tsfc[Elem_Idx,Line_Idx] 'tpw': value = L2.tpw[Elem_Idx,Line_Idx] 'logzopa': value = L2.LogZopa[Elem_Idx,Line_Idx] 'dzopasfc': value = L2.dZopaSfc[Elem_Idx,Line_Idx] 'dtsfcbt11': value = L2.dtsfcbt11[Elem_Idx,Line_Idx] 'dtsfcopa': value = L2.dtsfcopa[Elem_Idx,Line_Idx] 'logdtsfcopa': value = L2.logdtsfcopa[Elem_Idx,Line_Idx] 'bt62': value = L2.bt62[Elem_Idx,Line_Idx] 'bt11': value = L2.bt11[Elem_Idx,Line_Idx] 'bt133': value = L2.bt133[Elem_Idx,Line_Idx] 'btdclr10': value = L2.bt10_clr[Elem_Idx,Line_Idx] - L2.bt10[Elem_Idx,Line_Idx] 'btdclr11': value = L2.bt11_clr[Elem_Idx,Line_Idx] - L2.bt11[Elem_Idx,Line_Idx] 'btdclr12': value = L2.bt12_clr[Elem_Idx,Line_Idx] - L2.bt12[Elem_Idx,Line_Idx] 'bt11minsub': value = L2.bt11minsub[Elem_Idx,Line_Idx] 'bt11maxsub': value = L2.bt11maxsub[Elem_Idx,Line_Idx] 'bt10': value = L2.bt10[Elem_Idx,Line_Idx] 'dbtd3811clr': value = L2.dbtd3811clr[Elem_Idx,Line_Idx] 'btd3811': value = L2.btd3811[Elem_Idx,Line_Idx] 'btd3810': value = L2.btd3810[Elem_Idx,Line_Idx] 'emiss3811': value = L2.emiss3811[Elem_Idx,Line_Idx] 'emiss3810': value = L2.emiss3810[Elem_Idx,Line_Idx] 'btd7362': value = L2.btd7362[Elem_Idx,Line_Idx] 'btd1173': value = L2.btd1173[Elem_Idx,Line_Idx] 'btd11133': value = L2.btd11133[Elem_Idx,Line_Idx] 'btd13373': value = L2.btd13373[Elem_Idx,Line_Idx] 'btd8573': value = L2.bt85[Elem_Idx,Line_Idx] - L2.bt73[Elem_Idx,Line_Idx] 'btd8511': value = L2.btd8511[Elem_Idx,Line_Idx] 'btd1112': value = L2.btd1112[Elem_Idx,Line_Idx] 'btd1110': value = L2.btd1110[Elem_Idx,Line_Idx] 'fmft': value = L2.fmft[Elem_Idx,Line_Idx] 'bt11std': value = L2.bt11std[Elem_Idx,Line_Idx] 'logbt11std': value = L2.logbt11std[Elem_Idx,Line_Idx] 'bt10std': value = L2.bt10std[Elem_Idx,Line_Idx] 'logbt10std': value = L2.logbt10std[Elem_Idx,Line_Idx] 'dbt11max3x3': value = L2.dbt11max3x3[Elem_Idx,Line_Idx] 'dbt11maxsub': value = L2.dbt11maxsub[Elem_Idx,Line_Idx] 'dbt10max3x3': value = L2.dbt10max3x3[Elem_Idx,Line_Idx] 'logdbt11maxsub': value = L2.logdbt11maxsub[Elem_Idx,Line_Idx] 'logdbt11maxminsub': value = L2.logdbt11maxminsub[Elem_Idx,Line_Idx] 'logdbt11max3x3': value = L2.logdbt11max3x3[Elem_Idx,Line_Idx] 'logcod065': value = L2.logcod065[Elem_Idx,Line_Idx] 'logcod138': value = L2.logcod138[Elem_Idx,Line_Idx] 'refl047': value = L2.refl047[Elem_Idx,Line_Idx] 'refl065': value = L2.refl065[Elem_Idx,Line_Idx] 'refl138': value = L2.refl138[Elem_Idx,Line_Idx] 'refl160': value = L2.refl160[Elem_Idx,Line_Idx] 'refl38': value = L2.refl38[Elem_Idx,Line_Idx] 'refldnb': value = L2.refldnb[Elem_Idx,Line_Idx] 'logrefl065': value = L2.logrefl065[Elem_Idx,Line_Idx] 'logrefl138': value = L2.logrefl138[Elem_Idx,Line_Idx] 'logrefl160': value = L2.logrefl160[Elem_Idx,Line_Idx] 'drefl065clr': value = L2.drefl065clr[Elem_Idx,Line_Idx] 'drefl065maxsubclr': value = L2.drefl065maxsubclr[Elem_Idx,Line_Idx] 'drefl065min3x3': value = L2.drefl065min3x3[Elem_Idx,Line_Idx] 'logdrefl065minsub': value = L2.logdrefl065minsub[Elem_Idx,Line_Idx] 'logdrefl065maxminsub': value = L2.logdrefl065maxminsub[Elem_Idx,Line_Idx] 'logrefl065cv': value = L2.refl065cv[Elem_Idx,Line_Idx] 'drefl065minsub': value = L2.drefl065minsub[Elem_Idx,Line_Idx] 'refl065maxsub': value = L2.refl065maxsub[Elem_Idx,Line_Idx] 'drefl065maxminsub': value = L2.drefl065maxminsub[Elem_Idx,Line_Idx] 'refl065cv': value = L2.refl065cv[Elem_Idx,Line_Idx] 'refrat086065': value = L2.refrat086065[Elem_Idx,Line_Idx] 'refrat138065': value = L2.refrat138065[Elem_Idx,Line_Idx] 'refl065std': value = L2.refl065std[Elem_Idx,Line_Idx] 'logrefl065std': value = L2.logrefl065std[Elem_Idx,Line_Idx] 'ndsi': value = L2.ndsi[Elem_Idx,Line_Idx] 'bt6711covar' : value = L2.bt6711covar[Elem_Idx,Line_Idx] 'btd1167' : value = L2.btd1167[Elem_Idx,Line_Idx] 'bt11ratio': value = btratio(L2.bt11[Elem_Idx,Line_Idx],L2.tsfc[Elem_Idx,Line_Idx], $ L2.ttropo[Elem_Idx,Line_Idx], Missing) else: begin print, "unknown classifier dimension name: ", dim_name print, classifier_names[Class_Idx] Class_Use_Flag[Class_Idx] = 0 stop goto, skip_class end endcase case i of 0: x = value 1: y = value 2: z = value endcase endfor ;print, dim_name, value ;--- set prior from values in tables ;--- set prior if (Prior_Source eq 1) then begin Prior_Prob_Clear = 1.0 - Pixel_Prior_Prob_Cloud[Elem_Idx,Line_Idx] if (Pixel_Prior_Prob_Ice_Cloud[Elem_Idx] ne missing) then begin Prior_Prob_Ice = Pixel_Prior_Prob_Ice_Cloud[Elem_Idx,Line_Idx] Prior_Prob_Water = Pixel_Prior_Prob_Cloud[Elem_Idx,Line_Idx] - Pixel_Prior_Prob_Ice_Cloud[Elem_Idx,Line_Idx] endif else begin Prior_Prob_Water = 0.5*Pixel_Prior_Prob_Cloud[Elem_Idx,Line_Idx] Prior_Prob_Ice = 0.5*Pixel_Prior_Prob_Cloud[Elem_Idx,Line_Idx] endelse endif else begin Prior_Prob_Clear = Sfc_Prior_Prob_Clear[Sfc_Idx] ;note sfc_idx starts at zero Prior_Prob_Water = Sfc_Prior_Prob_Water[Sfc_Idx] Prior_Prob_Ice = Sfc_Prior_Prob_Ice[Sfc_Idx] Pixel_Prior_Prob_Cloud[Elem_Idx,Line_Idx] = 1.0 - Prior_Prob_Clear endelse ;--- constrain Prior_Prob_Clear = min([1.0,max([0,Prior_Prob_Clear])]) Prior_Prob_Water = min([1.0,max([0,Prior_Prob_Water])]) Prior_Prob_Ice = min([1.0,max([0,Prior_Prob_Ice])]) ;---- compute for this class if (Sfc_Idx lt 0 or Sfc_Idx gt 7) then goto, skip_pixel get_prob_mask_phase, x, $ y, $ z, $ l2.zen[Elem_Idx,Line_Idx],$ l2.solzen[Elem_Idx,Line_Idx], $ l2.lunzen[Elem_Idx,Line_Idx], $ l2.solglintzen[Elem_Idx,Line_Idx], $ l2.lunglintzen[Elem_Idx,Line_Idx], $ l2.solscatang[Elem_Idx,Line_Idx], $ l2.tpw[Elem_Idx,Line_Idx], $ l2.tsfc[Elem_Idx,Line_Idx], $ l2.zsfc[Elem_Idx,Line_Idx], $ l2.zsfcstd[Elem_Idx,Line_Idx], $ l2.land_class[Elem_Idx,Line_Idx], $ l2.snow_class[Elem_Idx,Line_Idx], $ l2.solglint_mask[Elem_Idx,Line_Idx], $ l2.lunglint_mask[Elem_Idx,Line_Idx], $ l2.coast_mask[Elem_Idx,Line_Idx], $ l2.city_mask[Elem_Idx,Line_Idx], $ l2.moon_illum_frac[Elem_Idx,Line_Idx], $ Sfc_Idx,$ *(luts[Class_Idx]), $ *(attrs[Class_Idx]), $ Missing, $ z1, $ z2, $ z3, $ z4 Clear_Cond_Ratio[Class_Idx] = z1 Water_Cond_Ratio[Class_Idx] = z2 Ice_Cond_Ratio[Class_Idx] = z3 Obs_Prob[Class_Idx] = z4 ;remove this? akh Weight = 1.0 ; if (z1 lt 0.0) then print, 'clear cond ratio negative ', z1 ; if (z2 lt 0.0) then print, 'water cond ratio negative ', z2 ; if (z3 lt 0.0) then print, 'water cond ratio negative ', z3 ;--- combine this class into the running product of all classes if (use_log_flag eq 0) then begin if (Clear_Cond_Ratio[Class_Idx] ne Missing) then begin R_Clear = max([R_Prod_min,min([R_Prod_Max,R_Clear * (Clear_Cond_Ratio[Class_Idx])])]) endif if (Water_Cond_Ratio[Class_Idx] ne Missing) then begin R_Water = max([R_Prod_min,min([R_Prod_Max,R_Water * Water_Cond_Ratio[Class_Idx]])]) endif if (Ice_Cond_Ratio[Class_Idx] ne Missing) then begin R_Ice = max([R_Prod_min,min([R_Prod_Max,R_Ice * Ice_Cond_Ratio[Class_Idx]])]) endif endif if (use_log_flag eq 1) then begin if (Clear_Cond_Ratio[Class_Idx] ne Missing) then begin R_Clear = max([R_Log_Min,min([R_Log_Max,R_Clear + (Clear_Cond_Ratio[Class_Idx])])]) endif if (Water_Cond_Ratio[Class_Idx] ne Missing) then begin R_Water = max([R_Log_Min,min([R_Log_Max,R_Water + Water_Cond_Ratio[Class_Idx]])]) endif if (Ice_Cond_Ratio[Class_Idx] ne Missing) then begin R_Ice = max([R_Log_Min,min([R_Log_Max,R_Ice + Ice_Cond_Ratio[Class_Idx]])]) endif endif ; if (z1 eq Missing) then begin ; print, 'Missing Clear_Cond_Ratio ', Class_Idx, z1,z2,z3,z4 ; endif ;--- store class value if (use_log_flag eq 0) then begin Clear_Cond_Ratio_Temp = Clear_Cond_Ratio[Class_Idx] endif if (use_log_flag eq 1) then begin Clear_Cond_Ratio_Temp = exp(Clear_Cond_Ratio[Class_Idx]) endif if (Clear_Cond_Ratio[Class_Idx] ne Missing) then begin Clear_Cond_Ratio_Temp = Clear_Cond_Ratio[Class_Idx] if (use_log_flag eq 1) then begin Clear_Cond_Ratio_Temp = exp(Clear_Cond_Ratio[Class_Idx]) endif Post_Prob_Clear_by_Class = $ 1.0 / (1.0 + Clear_Cond_Ratio_Temp/Prior_Prob_Clear - Clear_Cond_Ratio_Temp) ;-- convert from clear to cloud prob Post_Prob_Cloud_by_Class[Elem_Idx,Line_Idx,Class_Idx] = 1.0 -Post_Prob_Clear_by_Class endif ;print, 'test ', x, z1, R_Clear, R_Max, Post_Prob_Clear_by_Class, Prior_Prob_Clear, Post_Prob_Cloud_by_Class[Elem_Idx,Line_Idx,Class_Idx] skip_class: endfor ;end of Class_Idx Loop ;-------------------------------------------------------------------------------------- ; compute posterior probs for this pixel ;-------------------------------------------------------------------------------------- ;--- turn ratios into probabilities via post_prob = 1.0 / ( 1.0 + r/prior_yes - r) Post_Prob_Clear[Elem_Idx,Line_Idx] = 1.0 / (1.0 + R_Clear/Prior_Prob_Clear - R_Clear) Post_Prob_Water[Elem_Idx,Line_Idx] = 1.0 / (1.0 + R_Water/Prior_Prob_Water - R_Water) Post_Prob_Ice[Elem_Idx,Line_Idx] = 1.0 / (1.0 + R_Ice/Prior_Prob_Ice - R_Ice) ;--- normalize these probabilities so that they sum to 1 - won't happen if prior is not from lut Post_Sum = Post_Prob_clear[Elem_Idx,Line_Idx] + Post_Prob_Water[Elem_Idx,Line_Idx] + Post_Prob_Ice[Elem_Idx,Line_Idx] Post_Prob_Clear[Elem_Idx,Line_Idx] = Post_Prob_Clear[Elem_Idx,Line_Idx] / Post_Sum Post_Prob_Water[Elem_Idx,Line_Idx] = Post_Prob_Water[Elem_Idx,Line_Idx] / Post_Sum Post_Prob_Ice[Elem_Idx,Line_Idx] = Post_Prob_Ice[Elem_Idx,Line_Idx] / Post_Sum ;--- make a cloud probability from the complement of the clear probability Post_Prob_Cloud[Elem_Idx,Line_Idx] = 1.0 - Post_Prob_Clear[Elem_Idx,Line_Idx] ;--- make cloud phase and its uncertainty Post_Prob_Water_Norm = Post_Prob_Water[Elem_Idx,Line_Idx] / Post_Prob_Cloud[Elem_Idx,Line_Idx] Post_Prob_Ice_Norm = Post_Prob_Ice[Elem_Idx,Line_Idx] / Post_Prob_Cloud[Elem_Idx,Line_Idx] Cloud_Phase[Elem_Idx,Line_Idx] = 0 if (Post_Prob_Clear[Elem_Idx,Line_Idx] lt 0.50) then begin Cloud_Phase[Elem_Idx,Line_Idx] = 1 Cloud_Phase_Uncer[Elem_Idx,Line_Idx] = 1.0 - Post_Prob_Water_Norm if (Post_Prob_Ice_Norm gt Post_Prob_Water_Norm) then begin Cloud_Phase[Elem_Idx,Line_Idx] = 2 Cloud_Phase_Uncer[Elem_Idx,Line_Idx] = 1.0 - Post_Prob_Ice_Norm endif endif Cloud_Type = Cloud_Phase idx = where(Cloud_Phase eq 2,cc) if (cc gt 0) then Cloud_Type[idx] = 4 idx = where(Cloud_Phase eq 1 and L2.Topa lt 273.15,cc) if (cc gt 0) then Cloud_Type[idx] = 2 ;--- binary mask if (Post_Prob_Cloud[Elem_Idx,Line_Idx] lt Prob_Clear_Prob_Cloud_Thresh[Sfc_Idx] and $ Post_Prob_Cloud[Elem_Idx,Line_Idx] ne Missing) then begin Cloud_Mask_Binary[Elem_Idx,Line_Idx] = 0 endif if (Post_Prob_Cloud[Elem_Idx,Line_Idx] ge Prob_Clear_Prob_Cloud_Thresh[Sfc_Idx]) then begin Cloud_Mask_Binary[Elem_Idx,Line_Idx] = 1 endif ;--- 4-level mask if (Cloud_Mask_Binary[Elem_Idx,Line_Idx] eq 0) then begin Cloud_Mask[Elem_Idx,Line_Idx] = 0 if (Post_Prob_Cloud[Elem_Idx,Line_Idx] ge Conf_Clear_Prob_Clear_Thresh[Sfc_Idx]) then Cloud_Mask[Elem_Idx,Line_Idx] = 1 endif if (Cloud_Mask_Binary[Elem_Idx,Line_Idx] eq 1) then begin Cloud_Mask[Elem_Idx,Line_Idx] = 2 if (Post_Prob_Cloud[Elem_Idx,Line_Idx] ge Prob_Cloud_Conf_Cloud_Thresh[Sfc_Idx]) then Cloud_Mask[Elem_Idx,Line_Idx] = 3 endif ;--- Thermal Uniformity Test Filter (Needs to Use 10.4 when necessary) TUT[ELem_Idx,Line_Idx] = tut_routine(L2.Coast_Mask[Elem_Idx,Line_Idx], $ L2.bt11std[Elem_Idx,Line_Idx], $ L2.zsfcstd[Elem_Idx,Line_Idx], $ TUT_Clear_Prob_Clear_Thresh[Sfc_Idx]) ;--- Reflectance Uniformity Test Filter if (L2.Solzen[Elem_Idx,Line_Idx] lt attrs_temp.rut_solzen_thresh) then begin RUT_THRESH = RUT_CLEAR_PROB_CLEAR_THRESH[Sfc_Idx] Is_Land = 1 if (l2.sfc_type_mask[Elem_Idx,Line_Idx] eq 0 or l2.sfc_type_mask[Elem_Idx,Line_Idx] eq 1) then Is_Land = 0 RUT[Elem_Idx,Line_Idx] = RUT_ROUTINE(L2.Coast_Mask[Elem_Idx,Line_Idx], $ Is_Land, $ L2.Refl065[Elem_Idx,Line_Idx], $ L2.Refl065std[Elem_Idx,Line_Idx], $ RUT_THRESH) endif ;--- apply uniformity - bad idea ;if (Cloud_Mask[Elem_Idx,Line_Idx] eq 0 and $ ; (TUT[Elem_Idx,Line_Idx] eq 1 or RUT[Elem_Idx,Line_Idx] eq 1)) then Cloud_Mask[Elem_Idx,Line_Idx] = 1 skip_pixel: endfor endfor ;------------------------------------------------------------------------------------------------------- ; Imaging ;------------------------------------------------------------------------------------------------------- ;goto, skip_plot xsize0 = 1000 ysize0 = 300 load_color_table, 74, 1, 1.0, 249 window, 0, xsize = xsize0, ysize = ysize0 erase, color = 0 view, L2.refl065, /refl_image, /nocbar, /current,position=[0.0,0.0,0.334,1.0], /UP_DOWN, text_size = 0.2 ;view, L2.bt11, /t11_image, /nocbar, /current,position=[0.0,0.0,0.334,1.0], /UP_DOWN, text_size = 0.2 load_color_table, 74, 1, 1.0, 249 view, Post_Prob_Cloud, sds_min=-101, sds_max=1, position=[0.336,0.0,0.664,1.0], /current, /nocbar, /UP_DOWN, text_size = 0.2 ;view, L2.land_class, /land_image, position=[0.336,0.0,0.664,1.0], /current, /nocbar, /UP_DOWN load_color_table, 74, 1, 1.0, 249 view, l2.cld_prob, sds_min=-101, sds_max=1, position=[0.666,0.0,1.00,1.0], /current, /nocbar, /UP_DOWN, text_size = 0.2 load_color_table, 74, 1, 1.0, 249 legend_local, [Root_Name], position=[0.0,1.0],/norm, textcolor=5, bxcolor=1, /clear, charsize=1.5 saveimage, Root_Name+"_cld_prob_montage.png",/png wdelete xsize0 = 1000 ysize0 = 300 load_color_table, 74, 1, 1.0, 249 window, 1, xsize = xsize0, ysize = ysize0 erase, color = 0 view, L2.bt11, /t11_image, /nocbar, /current,position=[0.0,0.0,0.334,1.0], /UP_DOWN load_color_table, 1, 0, 1.0, 249 view, Cloud_Type, /type_image, position=[0.336,0.0,0.664,1.0], /current, /nocbar, /UP_DOWN ;view, Cloud_Phase, sds_min=-101, sds_max=2, position=[0.336,0.0,0.664,1.0], /current, /nocbar, /UP_DOWN view, l2.ctype, /type_image, position=[0.666,0.0,1.00,1.0], /current, /nocbar, /UP_DOWN load_color_table, 74, 1, 1.0, 249 legend_local, [Root_Name], position=[0.0,1.0],/norm, textcolor=5, bxcolor=1, /clear, charsize=1.5 saveimage, Root_Name+"_cld_phase_montage.png",/png wdelete xsize0 = 600 ysize0 = 600 view, L2.refl065, /refl_image, title = 'Refl_065', xsize=xsize0,ysize=ysize0,/UP_DOWN saveimage, Root_Name+"_refl065.png",/png wdelete ;xsize0 = 600 ;ysize0 = 600 ;view, L2.refldnb, /refl_image, title = 'Refl_DNB', xsize=xsize0,ysize=ysize0,/UP_DOWN ;saveimage, Root_Name+"_refldnb.png",/png ;wdelete view, L2.bt11, /t11_image, title = 'Bt_11', xsize=xsize0,ysize=ysize0,/UP_DOWN saveimage, Root_Name+"_bt11.png",/png wdelete view, L2.btd1112, coltab=0, sds_min=-2,sds_max=4, title = 'Btd_11_12', xsize=xsize0,ysize=ysize0,/UP_DOWN saveimage, Root_Name+"_btd1112.png",/png wdelete view, L2.sfc_type_mask, sds_min=-101,sds_max=7, coltab=74,title = 'mask surface type', xsize=xsize0,ysize=ysize0,/UP_DOWN saveimage, Root_Name+"_mask_sfc_type.png",/png wdelete view, L2.snow_class, /snow_image,title = 'snow class', xsize=xsize0,ysize=ysize0,/UP_DOWN saveimage, Root_Name+"_snow_class.png",/png wdelete view, L2.land_class, /land_image,title = 'snow class', xsize=xsize0,ysize=ysize0,/UP_DOWN saveimage, Root_Name+"_land_class.png",/png wdelete view, Pixel_Prior_Prob_Cloud, sds_min=-101, sds_max=1, coltab=74, rev=1, title = 'Prior_Prob_Cloud', xsize=xsize0,ysize=ysize0,/UP_DOWN saveimage, Root_Name+"_prior_cld_prob.png",/png wdelete view, Post_Prob_Cloud, sds_min=-101, sds_max=1, coltab=74, rev=1, title = 'Post_Prob_Cloud', xsize=xsize0,ysize=ysize0,/UP_DOWN saveimage, Root_Name+"_cld_prob.png",/png wdelete view, l2.cld_prob, sds_min=-101, sds_max=1, coltab=74, rev=1, title = 'Old Post_Prob_Cloud', xsize=xsize0,ysize=ysize0,/UP_DOWN saveimage, Root_Name+"_old_cld_prob.png",/png wdelete view, l2.cld_frac_nwp, sds_min=-101, sds_max=1, coltab=74, rev=1, title = 'NWP Cloud Frac',xsize=xsize0,ysize=ysize0, /UP_DOWN saveimage, Root_Name+"_cld_frac_nwp.png",/png wdelete view, Cloud_Mask, /mask_image, land_class=l2.land_class, title = '4-Level Cloud Mask', xsize=xsize0,ysize=ysize0,/UP_DOWN saveimage, Root_Name+"_cloud_mask.png",/png wdelete view, RUT, /frac_image, title = 'RUT', xsize=xsize0,ysize=ysize0,/UP_DOWN;, text_size = 0.2 saveimage, Root_Name+"_rut.png",/png wdelete view, TUT, /frac_image, title = 'TUT', xsize=xsize0,ysize=ysize0,/UP_DOWN;, text_size = 0.2 saveimage, Root_Name+"_tut.png",/png wdelete for Class_Idx = 0, Nclass - 1 do begin view, Post_Prob_Cloud_by_Class[*,*,Class_Idx], sds_min=-101, sds_max=1, coltab=74, rev=1, title = classifier_names[Class_Idx], xsize=xsize0,ysize=ysize0, /UP_DOWN saveimage, Root_Name+"_"+classifier_names[class_idx]+"_class_cld_prob.png",/png wdelete endfor skip_plot: save, file = Output_Name, l2, Post_Prob_Clear, Post_Prob_Water, Post_Prob_Ice, Post_Prob_Cloud, Post_Prob_Cloud_by_Class, $ Pixel_Prior_Prob_Cloud, Cloud_Mask_Binary, Cloud_Mask, ACMF, TUT, RUT, $ Cloud_Phase, Cloud_Phase_Uncer endfor ;end of L2 loop end