;--- ;------- ;--------------------- ;--------------------------------------------------------------- ;ecm_v1_5_routines ;--------------------------------------------------------------- ;--------------------- ;------- ;--- @compute_month_day @leapyr ;--------------------------------------------------------------- ; Detect Obivious Cloud - this works on arrays ; input: sfc_idx - the bayes surface type ; solzen - solar zenith angle ; glintzen - glint zenith angle ; t11_metric - ir window temperature metric ; ref065_metric - visible reflectance metric ; t11_rel ;--------------------------------------------------------------- function obvious_cloud_detection, sfc_idx, coast_mask, solzen, glintzen, t11_metric, ref065_metric, $, t11_rel_3x3_metric, t11_rel_sub_metric, t11_sub_flag nx = n_elements(solzen[*,0]) ny = n_elements(solzen[0,*]) missing = -999.0 post_prob_obvious = make_array(nx,ny,/FLOAT,VALUE=missing) ;--- apply 3x3 metric (t_max_3x3 - t) t11_rel_3x3_thresh_ocn = 3.0 t11_rel_3x3_thresh_land = 10.0 t11_rel_3x3_thresh_desert = 12.0 idx = where((Sfc_Idx eq 1 or Sfc_Idx eq 2) and coast_mask eq 0 and t11_rel_3x3_metric gt t11_rel_3x3_thresh_ocn, cc) if (cc gt 0) then post_prob_obvious[idx] = 1.0 idx = where(Sfc_Idx eq 3 and coast_mask eq 0 and t11_rel_3x3_metric gt t11_rel_3x3_thresh_land, cc) if (cc gt 0) then post_prob_obvious[idx] = 1.0 idx = where(Sfc_Idx eq 7 and coast_mask eq 0 and t11_rel_3x3_metric gt t11_rel_3x3_thresh_desert, cc) if (cc gt 0) then post_prob_obvious[idx] = 1.0 ;------ sub pixel t11 difference (t_max_sub - t_min_sub) t11_rel_sub_thresh_ocn = 2.0 t11_rel_sub_thresh_land = 5.0 t11_rel_sub_thresh_desert = 7.0 idx = where((Sfc_Idx eq 1 or Sfc_Idx eq 2) and coast_mask eq 0 and t11_rel_sub_metric gt t11_rel_sub_thresh_ocn, cc) if (cc gt 0) then post_prob_obvious[idx] = 1.0 idx = where(Sfc_Idx eq 3 and coast_mask eq 0 and t11_rel_sub_metric gt t11_rel_sub_thresh_land, cc) if (cc gt 0) then post_prob_obvious[idx] = 1.0 idx = where(Sfc_Idx eq 7 and coast_mask eq 0 and t11_rel_sub_metric gt t11_rel_sub_thresh_desert, cc) if (cc gt 0) then post_prob_obvious[idx] = 1.0 ;--- apply t11 metric idx = where((Sfc_Idx eq 1 or Sfc_Idx eq 2) and t11_metric lt 270.0 and t11_metric ne missing, cc) if (cc gt 0) then post_prob_obvious[idx] = 1.0 idx = where(Sfc_Idx eq 7 and t11_metric lt 240.0 and t11_metric ne missing, cc) if (cc gt 0) then post_prob_obvious[idx] = 1.0 idx = where(Sfc_Idx eq 3 and t11_metric lt 250.0 and t11_metric ne missing, cc) if (cc gt 0) then post_prob_obvious[idx] = 1.0 ;--- apply ref065 metric idx = where((Sfc_Idx eq 1 or Sfc_Idx eq 2) and solzen lt 60.0 and glintzen gt 40 and ref065_metric gt 15.0, cc) if (cc gt 0) then post_prob_obvious[idx] = 1.0 idx = where(Sfc_Idx eq 3 and solzen lt 60.0 and ref065_metric gt 40.0, cc) if (cc gt 0) then post_prob_obvious[idx] = 1.0 idx = where(Sfc_Idx eq 7 and solzen lt 60.0 and ref065_metric gt 40.0, cc) if (cc gt 0) then post_prob_obvious[idx] = 1.0 return, post_prob_obvious end ;--------------------------------------------------------------- ; prior prob into a structure ;--------------------------------------------------------------- function read_prior, prior_name ;--- Missing = -999.0 ;--- open cdfid = ncdf_open(prior_name) ;--- attributes ncdf_attget, cdfid, /GLOBAL, 'number_longitudes', Nlon ncdf_attget, cdfid, /GLOBAL, 'number_latitudes', Nlat ncdf_attget, cdfid, /GLOBAL, 'number_months', Nmonths ncdf_attget, cdfid, /GLOBAL, 'number_times', Ndiurnal ;--- structure Prior_Lut = {Nlon: Nlon, $ Nlat: Nlat, $ Nmonths: Nmonths, $ Ndiurnal: Ndiurnal, $ Dlon: 0.0, $ Dlat: 0.0, $ Lon_Min: 0.0, $ Lon_Max: 0.0, $ Lat_Min: 0.0, $ Lat_Max: 0.0, $ Cld_Frac: make_array(Nlon,Nlat,Nmonths,Ndiurnal,/FLOAT,VALUE=Missing)} ncdf_attget, cdfid, /GLOBAL, 'longitude_spacing', x & Prior_Lut.Dlon = 1 ncdf_attget, cdfid, /GLOBAL, 'latitude_spacing',x & Prior_Lut.Dlat = x ncdf_attget, cdfid, /GLOBAL, 'longitude_min',x & Prior_Lut.Lon_Min = x ncdf_attget, cdfid, /GLOBAL, 'longitude_max',x & Prior_Lut.Lon_Max = x ncdf_attget, cdfid, /GLOBAL, 'latitude_min',x & Prior_Lut.Lat_Min = x ncdf_attget, cdfid, /GLOBAL, 'latitude_max', x & Prior_Lut.Lat_Max = x ncdf_varget, cdfid, ncdf_varid(cdfid,'cloud_fraction_table_smoothed'), x Prior_Lut.Cld_Frac = x return, Prior_Lut ;--- close ncdf_close, cdfid end ;--------------------------------------------------------------- ; compute prior value ; ; input: Prior_Lut - prior cloud fraction structure ; Lon - longitude ; Lat - latitude ; Month - calendar month (1-12) ; output: Prior_Cld_Frac = the cloud fraction value from Prior_Lut ;--------------------------------------------------------------- function COMPUTE_PRIOR, Prior_Lut, Lon, Lat, Month Month_Idx = Month - 1 Diurnal_Idx = 0 ; daily averaged Lon_Idx = max([0,min([Prior_Lut.Nlon-1,fix((Lon - Prior_Lut.Lon_Min)/Prior_Lut.Dlon)])]) Lat_Idx = max([0,min([Prior_Lut.Nlat-1,fix((Lat - Prior_Lut.Lat_Min)/Prior_Lut.Dlat)])]) Prior_Cld_Frac = Prior_Lut.Cld_Frac[Lon_Idx,Lat_Idx,Month_Idx,Diurnal_Idx] return, Prior_Cld_Frac end ;--------------------------------------------------------------- ; read needed stuff in an "l2" structure ;--------------------------------------------------------------- function read_level2, level2_name forward_function read_and_check missing = -999.0 x = read_standard_hdf(level2_name,'latitude') Num_Elem_Idx = n_elements(x[*,0]) Num_Line_Idx = n_elements(x[0,*]) l2 = {Num_Elem_Idx: Num_Elem_Idx, $ Num_Line_Idx: Num_Line_Idx, $ year: 0, $ doy: 0, $ month: 0, $ dom: 0, $ lat: x, $ lon: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ solzen: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ glintzen: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ coast_mask: make_array(Num_Elem_Idx,Num_Line_Idx,/INT,VALUE=missing), $ land_class: make_array(Num_Elem_Idx,Num_Line_Idx,/INT,VALUE=missing), $ snow_class: make_array(Num_Elem_Idx,Num_Line_Idx,/INT,VALUE=missing), $ sfc_idx: make_array(Num_Elem_Idx,Num_Line_Idx,/INT,VALUE=missing), $ tsfc: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ zsfc: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ zsfc_stddev_3x3: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ ref065: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ ref065_clr: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ ref065_max_sub: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ ref065_min_sub: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ ref065_max_3x3: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ ref065_min_3x3: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ ref065_stddev_3x3: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ ref086: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ ref138: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ ref160: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ bt375: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ bt375_clr: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ bt85: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ bt85_clr: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ bt11: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ bt11_clr: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ bt11_min_sub: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ bt11_max_sub: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ bt11_max_3x3: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ bt11_stddev_3x3: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ bt12: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ bt12_clr: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ etrop: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ cld_prob: make_array(Num_Elem_Idx,Num_Line_Idx,/FLOAT,VALUE=missing), $ cld_mask: make_array(Num_Elem_Idx,Num_Line_Idx,/INT,VALUE=missing)} l2.lon = read_and_check(level2_name,'longitude',Num_Elem_Idx,Num_Line_Idx,Missing) l2.solzen = read_and_check(level2_name,'solar_zenith_angle',Num_Elem_Idx,Num_Line_Idx,Missing) l2.glintzen = read_and_check(level2_name,'glint_zenith_angle',Num_Elem_Idx,Num_Line_Idx,Missing) l2.tsfc = read_and_check(level2_name,'surface_temperature_nwp',Num_Elem_Idx,Num_Line_Idx,Missing) l2.zsfc = read_and_check(level2_name,'surface_elevation',Num_Elem_Idx,Num_Line_Idx,Missing) l2.zsfc_stddev_3x3 = read_and_check(level2_name,'surface_elevation_stddev_3x3',Num_Elem_Idx,Num_Line_Idx,Missing) l2.sfc_idx = read_and_check(level2_name,'bayes_mask_sfc_type',Num_Elem_Idx,Num_Line_Idx,Missing) l2.coast_mask = read_and_check(level2_name,'coast_mask',Num_Elem_Idx,Num_Line_Idx,Missing) l2.land_class = read_and_check(level2_name,'land_class',Num_Elem_Idx,Num_Line_Idx,Missing) l2.snow_class = read_and_check(level2_name,'snow_class',Num_Elem_Idx,Num_Line_Idx,Missing) l2.bt11 = read_and_check(level2_name,'temp_11_0um_nom',Num_Elem_Idx,Num_Line_Idx,Missing) l2.bt11_min_sub = read_and_check(level2_name,'temp_11_0um_nom_min_sub',Num_Elem_Idx,Num_Line_Idx,Missing) l2.bt11_max_sub = read_and_check(level2_name,'temp_11_0um_nom_max_sub',Num_Elem_Idx,Num_Line_Idx,Missing) l2.bt11_max_3x3 = read_and_check(level2_name,'temp_11_0um_nom_max_3x3',Num_Elem_Idx,Num_Line_Idx,Missing) l2.bt11_stddev_3x3 = read_and_check(level2_name,'temp_11_0um_nom_stddev_3x3',Num_Elem_Idx,Num_Line_Idx,Missing) l2.bt12 = read_and_check(level2_name,'temp_12_0um_nom',Num_Elem_Idx,Num_Line_Idx,Missing) l2.bt85 = read_and_check(level2_name,'temp_8_5um_nom',Num_Elem_Idx,Num_Line_Idx,Missing) l2.bt375 = read_and_check(level2_name,'temp_3_75um_nom',Num_Elem_Idx,Num_Line_Idx,Missing) l2.bt11_clr = read_and_check(level2_name,'temp_11_0um_nom_clear_sky',Num_Elem_Idx,Num_Line_Idx,Missing) l2.bt12_clr = read_and_check(level2_name,'temp_12_0um_nom_clear_sky',Num_Elem_Idx,Num_Line_Idx,Missing) l2.bt85_clr = read_and_check(level2_name,'temp_8_5um_nom_clear_sky',Num_Elem_Idx,Num_Line_Idx,Missing) l2.bt375_clr = read_and_check(level2_name,'temp_3_75um_nom_clear_sky',Num_Elem_Idx,Num_Line_Idx,Missing) l2.etrop = read_and_check(level2_name,'emiss_tropo_11_0um_nom',Num_Elem_Idx,Num_Line_Idx,Missing) l2.ref065 = read_and_check(level2_name,'refl_0_65um_nom',Num_Elem_Idx,Num_Line_Idx,Missing) l2.ref065_clr = read_and_check(level2_name,'refl_0_65um_nom_clear_sky',Num_Elem_Idx,Num_Line_Idx,Missing) l2.ref065_max_sub = read_and_check(level2_name,'refl_0_65um_nom_max_sub',Num_Elem_Idx,Num_Line_Idx,Missing) l2.ref065_min_sub = read_and_check(level2_name,'refl_0_65um_nom_min_sub',Num_Elem_Idx,Num_Line_Idx,Missing) l2.ref065_min_3x3 = read_and_check(level2_name,'refl_0_65um_nom_min_3x3',Num_Elem_Idx,Num_Line_Idx,Missing) l2.ref065_stddev_3x3 = read_and_check(level2_name,'refl_0_65um_nom_stddev_3x3',Num_Elem_Idx,Num_Line_Idx,Missing) l2.ref086 = read_and_check(level2_name,'refl_0_86um_nom',Num_Elem_Idx,Num_Line_Idx,Missing) l2.ref138 = read_and_check(level2_name,'refl_1_38um_nom',Num_Elem_Idx,Num_Line_Idx,Missing) l2.ref160 = read_and_check(level2_name,'refl_1_60um_nom',Num_Elem_Idx,Num_Line_Idx,Missing) l2.cld_prob = read_and_check(level2_name,'cloud_probability',Num_Elem_Idx,Num_Line_Idx,Missing) l2.cld_mask = read_and_check(level2_name,'cloud_mask',Num_Elem_Idx,Num_Line_Idx,Missing) ;--- attributes hdfid = hdf_sd_start(level2_name) hdf_sd_attrinfo, hdfid, hdf_sd_attrfind(hdfid,'START_YEAR'),data = year hdf_sd_attrinfo, hdfid, hdf_sd_attrfind(hdfid,'START_DAY'),data = doy hdf_sd_end, hdfid ;--- compute year, doy, dom, month l2.year = fix(year[0]) l2.doy = fix(doy[0]) compute_month_day, l2.doy, leapyr(l2.year), month, dom l2.month = month l2.dom = dom return, l2 end ;------------------------------------------------------------------------------------------ ; read level2 variables and check if chosen sds is missing and set to it missing ;------------------------------------------------------------------------------------------ function read_and_check, level2_name, sds_name,nx,ny,missing x = read_standard_hdf(level2_name,sds_name) if (n_elements(x) eq 1) then begin x = make_array(nx,ny,/FLOAT,VALUE=Missing) endif return, x end ;------------------------------------------------------------------------------------------------------------ ;--- make a cloud mask ;------------------------------------------------------------------------------------------------------------ pro MAKE_CLOUD_MASK, Sfc_Idx, Post_Prob, Rut_Test, Tut_Test, Cld_Mask, Cld_Mask_Binary ;--- based on lut training with 25% of clear and cloud as probably CONF_CLEAR_PROB_CLEAR_THRESH = [0.03, 0.03, 0.11, 0.28, 0.36, 0.38, 0.04] PROB_CLEAR_PROB_CLOUD_THRESH = [0.50, 0.50, 0.50, 0.50, 0.50, 0.50, 0.50] PROB_CLOUDY_CONF_CLOUD_THRESH = [0.98, 0.98, 0.96, 0.84, 0.82, 0.68, 0.97] YES = 1 NO = 0 ;--- mask values UNKNOWN = -1 CLEAR = 0 PROB_CLEAR = 1 PROB_CLOUDY = 2 CONF_CLOUDY = 3 CLEAR_BINARY = 0 CLOUDY_BINARY = 1 Cld_Mask = UNKNOWN Cld_Mask_Binary = UNKNOWN if (Sfc_Idx ge 1) then begin Cld_Mask = CLEAR Cld_Mask_Binary = CLEAR_BINARY if (Post_Prob ge PROB_CLOUDY_CONF_CLOUD_THRESH[Sfc_Idx-1]) then begin Cld_Mask = CONF_CLOUDY Cld_Mask_Binary = CLOUDY_BINARY endif if ((Post_Prob ge PROB_CLEAR_PROB_CLOUD_THRESH[Sfc_Idx-1]) and $ (Post_Prob lt PROB_CLOUDY_CONF_CLOUD_THRESH[Sfc_Idx-1])) then begin Cld_Mask = PROB_CLOUDY Cld_Mask_Binary = CLOUDY_BINARY endif if ((Post_Prob gt CONF_CLEAR_PROB_CLEAR_THRESH[Sfc_Idx-1]) and $ (Post_Prob lt PROB_CLEAR_PROB_CLOUD_THRESH[Sfc_Idx-1])) then begin Cld_Mask = PROB_CLEAR Cld_Mask_Binary = CLEAR_BINARY endif endif ;----- include RUT and TUT filter on Clear if (Cld_Mask eq CLEAR and (Rut_Test eq YES or Tut_Test eq YES)) then begin Cld_Mask = PROB_CLEAR endif ;---- use spatial analysis filter for probably cloudy end ;==================================================================== ; Function Name: TUT_Routine ; ; Function: ; Thermal Uniformity Routine (TUT) Test ; ; Description: ; Computes the Thermal Uniformity Routine (TUT) test (YES/NO), as described in ; the ABI cloud mask ATBD ; ; ; Calling Sequence: ; Test_Results(Num_Tests,Elem_Idx,Line_Idx) = TUT_Routine ($ ; Is_Land, $ ; Is_Coast, $ ; BT_IRWIN_Stddev_3x3, $ ; Sfc_Hgt_Stddev_3x3) ; ; Inputs: ; Is pixel land (YES/NO) ; Is pixel coast (YES/NO) ; Standard Deviation of the 11 micron BT over a 3x3 box ; Standard Deviation of the surface height over a 3x3 box ; ; Outputs: ; Function returns pass (sym%YES) or fail (sym%NO) result of the test via Test_Result ; ; Dependencies: ; Cloud Mask threshold include file with various needed thresholds. ; ; Restrictions: None ; ;==================================================================== FUNCTION TUT_ROUTINE,Is_Land, $ Is_Coast, $ BT_IRWIN_Stddev_3x3, $ Sfc_Hgt_Stddev_3x3 NO = 0 YES = 1 BT_IRWIN_CLR_UNI_THRESH_LAND = 1.1 BT_IRWIN_CLR_UNI_THRESH_OCN = 0.6 Test_Result = NO IF (Is_Coast eq NO) THEN BEGIN ; ;7K/km is the adiabatic lapse rate ; Test_Threshold = 3.0 * 7.0*Sfc_Hgt_Stddev_3x3/1000.0 IF (Is_Land eq YES) THEN BEGIN IF (BT_IRWIN_Stddev_3x3 gt BT_IRWIN_CLR_UNI_THRESH_LAND $ +Test_Threshold) THEN BEGIN Test_Result = YES ENDIF ENDIF ELSE BEGIN IF (BT_IRWIN_Stddev_3x3 gt BT_IRWIN_CLR_UNI_THRESH_OCN $ +Test_Threshold) THEN BEGIN Test_Result = YES ENDIF ENDELSE ENDIF RETURN, Test_Result END ;==================================================================== ; Function Name: RUT_Routine ; ; Function: ; Reflectance Uniformity Routine (RUT) Test ; ; Description: ; Computes the Reflectance Uniformity Routine (RUT) test (YES/NO), as described in ; the ABI cloud mask ATBD ; ; ; Calling Sequence: ; Test_Results(Num_Tests,Elem_Idx,Line_Idx) = RUT_Routine ($ ; Is_Land, $ ; Is_Snow, $ ; Is_Coast, $ ; Refl_Vis_Clear(Elem_Idx,Line_Idx), $ ; Refl_Vis_Stddev_3x3(Elem_Idx,Line_Idx), $ ; Sol_Zen) ; ; Inputs: ; Is pixel land (YES/NO) ; Is pixel snow (YES/NO) ; Is pixel coast (YES/NO) ; 0.64 miron clear sky reflectance ; Standard Deviation of the 0.64 micron reflectance over a 3x3 box ; Solar zenith angle ; ; Outputs: ; Function returns pass (sym%YES) or fail (sym%NO) result of the test (Test_Result) ; ; Dependencies: ; Cloud Mask threshold include file with various needed thresholds. ; ; Restrictions: None ; ;==================================================================== FUNCTION RUT_ROUTINE,Is_Land, $ Is_Snow, $ Is_Coast, $ Refl_Vis_Clear, $ Refl_Vis_Stddev_3x3, $ Sol_Zen NO = 0 YES = 1 RUT_SOL_ZEN_THRESH = 80.0 REFL_VIS_CLR_UNI_THRESH_LAND = 0.20 REFL_VIS_CLR_UNI_THRESH_OCEAN = 1.0 Test_Result = NO IF (Is_Snow eq NO AND Is_Coast eq NO) THEN BEGIN IF (Sol_Zen lt RUT_SOL_ZEN_THRESH) THEN BEGIN ;--- compute threshold IF (Is_Land eq YES) THEN BEGIN Test_Threshold = MAX([0.5,Refl_Vis_Clear * REFL_VIS_CLR_UNI_THRESH_LAND]) ENDIF ELSE BEGIN Test_Threshold = REFL_VIS_CLR_UNI_THRESH_OCEAN ENDELSE ;--- apply test IF (Refl_Vis_Stddev_3x3 gt Test_Threshold) THEN BEGIN Test_Result = YES ENDIF ENDIF ENDIF RETURN, Test_Result END