From a73997136268531702a4ce1c8ed69f5854f4f108 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Wed, 15 Nov 2023 14:38:35 -0500 Subject: [PATCH 1/6] add input parameters for cloud filtering --- src/Applications/GAAS_App/modis_l2a.pcf | 4 ++++ src/Applications/GAAS_App/modis_l2a.py | 20 +++++++++++++++++- src/Applications/GAAS_App/mxd04_l2a.py | 28 +++++++++++++++++++++---- src/Applications/GAAS_App/viirs_l2a.pcf | 4 ++++ src/Applications/GAAS_App/viirs_l2a.py | 22 +++++++++++++++++-- src/Applications/GAAS_App/vx04_l2a.py | 28 +++++++++++++++++++++---- 6 files changed, 95 insertions(+), 11 deletions(-) diff --git a/src/Applications/GAAS_App/modis_l2a.pcf b/src/Applications/GAAS_App/modis_l2a.pcf index a5d9e55c..eb729a13 100644 --- a/src/Applications/GAAS_App/modis_l2a.pcf +++ b/src/Applications/GAAS_App/modis_l2a.pcf @@ -22,4 +22,8 @@ MODIS_L2A_AER_X = %s.gaas_bkg.sfc.%y4%m2%d2_%h2z.nc4 MODIS_L2A_NN_FILE = ExtData/g5chem/x/NN/nnr_003.%ident_Tau.net MODIS_L2A_BLANK_ODS = ExtData/g5chem/x/blank_syn8.ods +MODIS_L2A_CLOUD_THRESH = 0.7 +MODIS_L2A_CLOUDFREE = 0.1 +MODIS_L2A_AODMAX = 10 + #END MODIS_L2A processing diff --git a/src/Applications/GAAS_App/modis_l2a.py b/src/Applications/GAAS_App/modis_l2a.py index 3c94ee70..a974e9af 100755 --- a/src/Applications/GAAS_App/modis_l2a.py +++ b/src/Applications/GAAS_App/modis_l2a.py @@ -47,7 +47,25 @@ if cf('MODIS_L2A_OVERWRITE').upper() == 'YES': Options += " --force" if cf('MODIS_L2A_VERBOSE').upper() == 'YES': Options += " -v" - + + try: + aodmax = cf('MODIS_L2A_AODMAX') + Options += " --aodmax=" + aodmax + except: + pass + + try: + cloudFree = cf('MODIS_L2A_CLOUDFREE') + Options += " --cloudFree=" + cloudFree + except: + pass + + try: + cloud_thresh = cf('MODIS_L2A_CLOUD_THRESH') + Options += " --cloud_thresh=" + cloud_thresh + except: + pass + # Generate products # ----------------- i = 0 diff --git a/src/Applications/GAAS_App/mxd04_l2a.py b/src/Applications/GAAS_App/mxd04_l2a.py index 6a1f44c1..34df1e74 100755 --- a/src/Applications/GAAS_App/mxd04_l2a.py +++ b/src/Applications/GAAS_App/mxd04_l2a.py @@ -71,6 +71,9 @@ def makethis_dir(filename): coll = '006' res = 'c' nsyn = 8 + aodmax = 1.0 + cloud_thresh = 0.7 + cloudFree = None # Parse command line options # -------------------------- @@ -120,7 +123,19 @@ def makethis_dir(filename): parser.add_option("-r", "--res", dest="res", default=res, help="Resolution for gridded output (default=%s)"\ - %out_tmpl ) + %res ) + + parser.add_option("--cloud_thresh", dest="cloud_thresh", default=cloud_thresh,type='float', + help="Cloud fractions threshhold for good data (default=%f)"\ + %cloud_thresh ) + + parser.add_option("--cloudFree", dest="cloudFree", default=cloudFree, + help="Extra check for cloudiness when high AOD values are predicted. If not provided, no check is performed. (default=%s)"\ + %cloudFree ) + + parser.add_option("--aodmax", dest="aodmax", default=aodmax,type='float', + help="max AOD value that will be accepted when cloud fraction is greater than cloudFree (default=%f)"\ + %aodmax ) parser.add_option("-u", "--uncompressed", action="store_true", dest="uncompressed",default=False, @@ -198,11 +213,16 @@ def makethis_dir(filename): if options.verbose: print("NNR Retrieving %s %s on "%(prod,algo.upper()),syn_time) + if options.cloudFree == 'None': + options.cloudFree = None + else: + options.cloudFree = float(options.cloudFree) + modis = MxD04_NNR(options.l2_path,prod,algo.upper(),syn_time,aer_x,slv_x, coll=options.coll, - cloud_thresh=0.7, - cloudFree = 0.0, - aodmax = 1.0, + cloud_thresh=options.cloud_thresh, + cloudFree = options.cloudFree, + aodmax = options.aodmax, nsyn=options.nsyn, verbose=options.verbose) if modis.nobs < 1: diff --git a/src/Applications/GAAS_App/viirs_l2a.pcf b/src/Applications/GAAS_App/viirs_l2a.pcf index c6c889ce..74d3df9c 100644 --- a/src/Applications/GAAS_App/viirs_l2a.pcf +++ b/src/Applications/GAAS_App/viirs_l2a.pcf @@ -19,4 +19,8 @@ VIIRS_L2A_AER_X = /nobackup/NNR/Misc/tavg1_2d_aer_Nx VIIRS_L2A_NN_FILE = nnr_001.%ident_TauAE.net VIIRS_L2A_BLANK_ODS = /nobackup/NNR/Misc/blank.ods +VIIRS_L2A_CLOUD_THRESH = 0.7 +VIIRS_L2A_CLOUDFREE = 0.1 +VIIRS_L2A_AODMAX = 10 + #END VIIRS_L2A processing diff --git a/src/Applications/GAAS_App/viirs_l2a.py b/src/Applications/GAAS_App/viirs_l2a.py index f211fcbf..deb18767 100755 --- a/src/Applications/GAAS_App/viirs_l2a.py +++ b/src/Applications/GAAS_App/viirs_l2a.py @@ -46,11 +46,29 @@ " --fname=" + cf('VIIRS_L2A_OUT_TEMPLATE') + \ " --net=" + cf('VIIRS_L2A_NN_FILE') + \ " --aer_x=" + cf('VIIRS_L2A_AER_X') + \ - " --blank_ods=" + cf('VIIRS_L2A_BLANK_ODS') + " --blank_ods=" + cf('VIIRS_L2A_BLANK_ODS') if cf('VIIRS_L2A_OVERWRITE').upper() == 'YES': Options += " --force" if cf('VIIRS_L2A_VERBOSE').upper() == 'YES': Options += " -v" - + + try: + aodmax = cf('VIIRS_L2A_AODMAX') + Options += " --aodmax=" + aodmax + except: + pass + + try: + cloudFree = cf('VIIRS_L2A_CLOUDFREE') + Options += " --cloudFree=" + cloudFree + except: + pass + + try: + cloud_thresh = cf('VIIRS_L2A_CLOUD_THRESH') + Options += " --cloud_thresh=" + cloud_thresh + except: + pass + # Generate products # ----------------- i = 0 diff --git a/src/Applications/GAAS_App/vx04_l2a.py b/src/Applications/GAAS_App/vx04_l2a.py index 7de56401..80441392 100755 --- a/src/Applications/GAAS_App/vx04_l2a.py +++ b/src/Applications/GAAS_App/vx04_l2a.py @@ -67,6 +67,9 @@ def makethis_dir(filename): coll = '002' res = 'c' nsyn = 8 + cloud_thresh = 0.70 + cloudFree = None + aodmax = 1.0 # Parse command line options # -------------------------- @@ -112,7 +115,19 @@ def makethis_dir(filename): parser.add_option("-r", "--res", dest="res", default=res, help="Resolution for gridded output (default=%s)"\ - %out_tmpl ) + %res ) + + parser.add_option("--cloud_thresh", dest="cloud_thresh", default=cloud_thresh,type='float', + help="Cloud fractions threshhold for good data (default=%f)"\ + %cloud_thresh ) + + parser.add_option("--cloudFree", dest="cloudFree", default=cloudFree, + help="Extra check for cloudiness when high AOD values are predicted. If not provided, no check is performed. (default=%s)"\ + %cloudFree ) + + parser.add_option("--aodmax", dest="aodmax", default=aodmax,type='float', + help="max AOD value that will be accepted when cloud fraction is greater than cloudFree (default=%f)"\ + %aodmax ) parser.add_option("-u", "--uncompressed", action="store_true", dest="uncompressed",default=False, @@ -188,11 +203,16 @@ def makethis_dir(filename): if options.verbose: print("NNR Retrieving %s %s on "%(sat,algo.upper()),syn_time) + if options.cloudFree == 'None': + options.cloudFree = None + else: + options.cloudFree = float(options.cloudFree) + viirs = Vx04_NNR(options.l2_path,sat,algo.upper(),syn_time,aer_x, coll=options.coll, - cloud_thresh=0.7, - cloudFree = 0.0, - aodmax = 1.0, + cloud_thresh=options.cloud_thresh, + cloudFree=options.cloudFree, + aodmax=options.aodmax, nsyn=options.nsyn, verbose=options.verbose) if viirs.nobs < 1: From ab439aeafef0368f87c65e703d7f0d5e0ce5e6e3 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Thu, 16 Nov 2023 17:25:06 -0500 Subject: [PATCH 2/6] modified to handle 1 parameter training of AEfit --- src/Applications/GAAS_App/mxd04_nnr.py | 6 ++++++ src/Applications/GAAS_App/vx04_nnr.py | 6 ++++++ 2 files changed, 12 insertions(+) diff --git a/src/Applications/GAAS_App/mxd04_nnr.py b/src/Applications/GAAS_App/mxd04_nnr.py index 9bf76bb5..4ac12251 100644 --- a/src/Applications/GAAS_App/mxd04_nnr.py +++ b/src/Applications/GAAS_App/mxd04_nnr.py @@ -483,11 +483,17 @@ def apply(self,nnFile): wavs = ['440','470','550','660','870'] wav = np.array(wavs).astype(float) nwav = len(wavs) + AEfitb = None for i,targetName in enumerate(self.net.TargetNames): if 'AEfitm' in targetName: AEfitm = targets[:,i] if 'AEfitb' in targetName: AEfitb = targets[:,i] + if 'aTau550' in targetName: + tau550 = targets[:,i] + + if AEfitb == None: + AEfitb = -1.*(tau550 + AEfitm*np.log(550.)) nobs = targets.shape[0] targets_ = np.zeros([nobs,nwav]) targetName = [] diff --git a/src/Applications/GAAS_App/vx04_nnr.py b/src/Applications/GAAS_App/vx04_nnr.py index d8d74ffa..669039b4 100644 --- a/src/Applications/GAAS_App/vx04_nnr.py +++ b/src/Applications/GAAS_App/vx04_nnr.py @@ -338,11 +338,17 @@ def apply(self,nnFile): wavs = ['440','470','550','660','870'] wav = np.array(wavs).astype(float) nwav = len(wavs) + AEfitb = None for i,targetName in enumerate(self.net.TargetNames): if 'AEfitm' in targetName: AEfitm = targets[:,i] if 'AEfitb' in targetName: AEfitb = targets[:,i] + if 'aTau550' in targetName: + tau550 = targets[:,i] + + if AEfitb == None: + AEfitb = -1.*(tau550 + AEfitm*np.log(550.)) nobs = targets.shape[0] targets_ = np.zeros([nobs,nwav]) targetName = [] From fd0738d4e2c6212cccf0a37d69a962bc556cecad Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Fri, 17 Nov 2023 11:48:43 -0500 Subject: [PATCH 3/6] mxd04_nnr changes for cloud and outlier filtering --- src/Applications/GAAS_App/mxd04_nnr.py | 66 ++++++++++++++++++++++++-- 1 file changed, 62 insertions(+), 4 deletions(-) diff --git a/src/Applications/GAAS_App/mxd04_nnr.py b/src/Applications/GAAS_App/mxd04_nnr.py index 4ac12251..fc4709a9 100644 --- a/src/Applications/GAAS_App/mxd04_nnr.py +++ b/src/Applications/GAAS_App/mxd04_nnr.py @@ -118,7 +118,9 @@ def __init__(self,l2_path,prod,algo,syn_time,aer_x,slv_x, glint_thresh=40.0, scat_thresh=170.0, cloudFree=None, - aodmax=1.0, + aodmax=2.0, + aodSTD=3.0, + aodLength=0.5, coll='006', nsyn=8, verbose=0): @@ -136,6 +138,8 @@ def __init__(self,l2_path,prod,algo,syn_time,aer_x,slv_x, cloud_tresh --- cloud fraction treshhold cloudFree --- cloud fraction threshhold for assuring no cloud contaminations when aod is > aodmax if None, no cloud free check is made + aodSTD --- number of standard deviations for checking for outliers + aodLength --- length scale (degrees) to look for outliers coll --- MODIS data collection nsyn --- number of synoptic times @@ -154,7 +158,9 @@ def __init__(self,l2_path,prod,algo,syn_time,aer_x,slv_x, self.algo = algo self.cloudFree = cloudFree self.aodmax = aodmax - + self.aodSTD = aodSTD + self.aodLength = aodLength + # Initialize superclass # --------------------- Files = granules(l2_path,prod,syn_time,coll=coll,nsyn=nsyn) @@ -573,6 +579,7 @@ def apply(self,nnFile): # Do extra cloud filtering if required if self.cloudFree is not None: + # start by checking the cloud masks if self.algo == "LAND": cloudy = (self.cloud_deep>=self.cloudFree) & (self.cloud>=self.cloudFree) elif self.algo == "DEEP": @@ -580,24 +587,75 @@ def apply(self,nnFile): elif self.algo == "OCEAN": cloudy = (self.cloud>=self.cloudFree) + # if cloud fraction exceeds cloudFree and the aod exceeds aodmax, filter out contaminated = np.zeros(np.sum(self.iGood)).astype(bool) for targetName in self.net.TargetNames: name, ch = TranslateTarget[targetName] k = list(self.channels).index(ch) # index of channel result = self.__dict__[name][self.iGood,k] contaminated = contaminated | ( (result > self.aodmax) & cloudy[self.iGood] ) - + + if self.verbose: + print('Filtering out ',np.sum(contaminated),' suspected cloud contaminated pixels') + for targetName in self.net.TargetNames: name, ch = TranslateTarget[targetName] k = list(self.channels).index(ch) # index of channel self.__dict__[name][self.iGood,k][contaminated] = MISSING if doAEfit: - self.ae[self.iGood][contaminated] = MISSING self.ae_[self.iGood][contaminated] = MISSING self.iGood[self.iGood][contaminated] = False + # check for outliers + # start with highest AOD550 value + # find all the pixels within a 1 degree neighborhood + # check if it is outside of mean + N*sigma of the other pixels + # aodSTD parameter is equal to N + # continue until no outliers are found + find_outliers = True + k = list(self.channels).index(550) + aod550 = np.ma.array(self.aod_[self.iGood,k]) + aod550.mask = np.zeros(len(aod550)).astype(bool) + Lon = self.Longitude[self.iGood] + Lat = self.Latitude[self.iGood] + gIndex = np.arange(len(self.iGood))[self.iGood] + iOutliers = [] + while find_outliers: + maxaod = aod550.max() + imax = np.argmax(aod550) + aod550.mask[imax] = True + lon = Lon[imax] + lat = Lat[imax] + + # find the neighborhood of pixels + iHood = (Lon<=lon+self.aodLength) & (Lon>=lon-self.aodLength) & (Lat<=lat+self.aodLength) & (Lat>=lat-self.aodLength) + if np.sum(iHood) <= 1: + #this pixel has no neighbors and is high. Filter it. + iOutliers.append(gIndex[imax]) + else: + aodHood = aod550[iHood] + if maxaod > (aodHood.mean() + self.aodSTD*aodHood.std()): + iOutliers.append(gIndex[imax]) + else: + find_outliers = False # done looking for outliers + if self.verbose: + print("Filtering out ",len(iOutliers)," outlier pixels") + + self.iOutliers = iOutliers + + if len(iOutliers) > 0: + for targetName in self.net.TargetNames: + name, ch = TranslateTarget[targetName] + k = list(self.channels).index(ch) # index of channel + self.__dict__[name][iOutliers,k] = MISSING + + if doAEfit: + self.ae_[iOutliers] = MISSING + + self.iGood[iOutliers] = False + #--- From 8e5c9d2cc8533f11f4ea706c04df9037af39d155 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Fri, 17 Nov 2023 12:00:19 -0500 Subject: [PATCH 4/6] fix to get rid of numpy warning when doing polyfit of modis AE --- src/Applications/GAAS_App/mxd04_nnr.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/Applications/GAAS_App/mxd04_nnr.py b/src/Applications/GAAS_App/mxd04_nnr.py index fc4709a9..bab583fb 100644 --- a/src/Applications/GAAS_App/mxd04_nnr.py +++ b/src/Applications/GAAS_App/mxd04_nnr.py @@ -519,11 +519,13 @@ def apply(self,nnFile): I = np.array(self.channels) < 900 # only visible channels, this is relevant for ocean aechannels = np.array(self.channels)[I] aodT = self.aod[:,I].T - fit = np.polyfit(np.log(aechannels),-1.*np.log(aodT[:,self.iGood]+0.01),1) + iIndex = np.arange(len(self.iGood))[self.iGood] + aodT = aodT[:,iIndex] + 0.01 + mask = aodT.min(axis=0) > 0 + posIndex = iIndex[mask] + fit = np.polyfit(np.log(aechannels),-1.*np.log(aodT[:,mask]+0.01),1) self.ae = MISSING*ones(self.nobs) - self.ae[self.iGood] = fit[0,:] - bad = np.isnan(self.ae) - self.ae[bad] = MISSING + self.ae[posIndex] = fit[0,:] if doAE: From 88ae34425dc35c0f40812834e48205dff3e1b9d9 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Fri, 17 Nov 2023 13:32:58 -0500 Subject: [PATCH 5/6] clean up one thing to use self.nobs instead of len(self.iGood). cleaner --- src/Applications/GAAS_App/mxd04_nnr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/GAAS_App/mxd04_nnr.py b/src/Applications/GAAS_App/mxd04_nnr.py index 69556d35..1d42ecbc 100644 --- a/src/Applications/GAAS_App/mxd04_nnr.py +++ b/src/Applications/GAAS_App/mxd04_nnr.py @@ -626,7 +626,7 @@ def apply(self,nnFile): aod550.mask = np.zeros(len(aod550)).astype(bool) Lon = self.Longitude[self.iGood] Lat = self.Latitude[self.iGood] - gIndex = np.arange(len(self.iGood))[self.iGood] + gIndex = np.arange(self.nobs)[self.iGood] iOutliers = [] while find_outliers: maxaod = aod550.max() From 65449c2d597c3c85a42a71e5889ee56483879a92 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Thu, 30 Nov 2023 12:43:00 -0500 Subject: [PATCH 6/6] update changelog for updates to MODIS and VIIRS NNR updated cloud filter --- CHANGELOG.md | 1 + src/Applications/GAAS_App/modis_l2a.pcf | 7 +- src/Applications/GAAS_App/modis_l2a.py | 27 ++++++-- src/Applications/GAAS_App/mxd04_l2a.py | 25 +++++++- src/Applications/GAAS_App/mxd04_nnr.py | 16 +++-- src/Applications/GAAS_App/viirs_l2a.pcf | 7 +- src/Applications/GAAS_App/viirs_l2a.py | 28 ++++++-- src/Applications/GAAS_App/vx04_l2a.py | 23 ++++++- src/Applications/GAAS_App/vx04_nnr.py | 85 +++++++++++++++++++++---- 9 files changed, 186 insertions(+), 33 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 066f4a2d..ecacff9c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ - VIIRS and MODIS NNR angstrom exponent predictions are now saved in gridded output files (Level 3) - VIIRS and MODIS NNR ods files of AOD can be generated from angstrom exponent fit coefficients - VIIRS NNR ods file generator now takes nsyn (number of synoptic times) as an input. +- VIIRS & MODIS NNR codes have an improved cloud filter and analomaly detection algorithm. There parameters for this can be inputs in the pcf file. ### Fixed diff --git a/src/Applications/GAAS_App/modis_l2a.pcf b/src/Applications/GAAS_App/modis_l2a.pcf index eb729a13..32abf8b7 100644 --- a/src/Applications/GAAS_App/modis_l2a.pcf +++ b/src/Applications/GAAS_App/modis_l2a.pcf @@ -23,7 +23,10 @@ MODIS_L2A_NN_FILE = ExtData/g5chem/x/NN/nnr_003.%ident_Tau.net MODIS_L2A_BLANK_ODS = ExtData/g5chem/x/blank_syn8.ods MODIS_L2A_CLOUD_THRESH = 0.7 -MODIS_L2A_CLOUDFREE = 0.1 -MODIS_L2A_AODMAX = 10 +MODIS_L2A_CLOUDFREE = 0.5 +MODIS_L2A_AODMAX = 2.0 +MODIS_L2A_AODSTD = 3.0 +MODIS_L2A_AODLENGTH = 0.5 +MODIS_L2A_WAVS = 440,470,550,660,870 #END MODIS_L2A processing diff --git a/src/Applications/GAAS_App/modis_l2a.py b/src/Applications/GAAS_App/modis_l2a.py index a974e9af..d75f8931 100755 --- a/src/Applications/GAAS_App/modis_l2a.py +++ b/src/Applications/GAAS_App/modis_l2a.py @@ -48,9 +48,10 @@ if cf('MODIS_L2A_OVERWRITE').upper() == 'YES': Options += " --force" if cf('MODIS_L2A_VERBOSE').upper() == 'YES': Options += " -v" + try: - aodmax = cf('MODIS_L2A_AODMAX') - Options += " --aodmax=" + aodmax + cloud_thresh = cf('MODIS_L2A_CLOUD_THRESH') + Options += " --cloud_thresh=" + cloud_thresh except: pass @@ -61,11 +62,29 @@ pass try: - cloud_thresh = cf('MODIS_L2A_CLOUD_THRESH') - Options += " --cloud_thresh=" + cloud_thresh + aodmax = cf('MODIS_L2A_AODMAX') + Options += " --aodmax=" + aodmax except: pass + try: + aodSTD = cf('MODIS_L2A_AODSTD') + Options += " --aodSTD=" + aodSTD + except: + pass + + try: + aodLenth = cf('MODIS_L2A_AODLENGTH') + Options += " --aodLength=" + aodLength + except: + pass + + try: + wavs = cf('MODIS_L2A_WAVS') + Options += " --wavs=" + wavs + except: + pass + # Generate products # ----------------- i = 0 diff --git a/src/Applications/GAAS_App/mxd04_l2a.py b/src/Applications/GAAS_App/mxd04_l2a.py index 34df1e74..21d46411 100755 --- a/src/Applications/GAAS_App/mxd04_l2a.py +++ b/src/Applications/GAAS_App/mxd04_l2a.py @@ -44,7 +44,8 @@ def makethis_dir(filename): raise IOError("could not create directory "+path) #--------------------------------------------------------------------- - +def wavs_callback(option, opt, value, parser): + setattr(parser.values, option.dest, value.split(',')) if __name__ == "__main__": expid = 'nnr3' @@ -71,9 +72,12 @@ def makethis_dir(filename): coll = '006' res = 'c' nsyn = 8 - aodmax = 1.0 + aodmax = 2.0 cloud_thresh = 0.7 cloudFree = None + aodSTD = 3.0 + aodLength = 0.5 + wavs = '440,470,550,660,870' # Parse command line options # -------------------------- @@ -137,6 +141,18 @@ def makethis_dir(filename): help="max AOD value that will be accepted when cloud fraction is greater than cloudFree (default=%f)"\ %aodmax ) + parser.add_option("--aodSTD", dest="aodSTD", default=aodSTD,type='float', + help="number of standard deviations to check for AOD outliers (default=%f)"\ + %aodSTD ) + + parser.add_option("--aodLength", dest="aodLength", default=aodLength,type='float', + help="length scale (degrees) to check for AOD outliers (default=%f)"\ + %aodLength ) + + parser.add_option("--wavs", dest="wavs", default=wavs,type='string',action='callback',callback=wavs_callback, + help="wavelength to output AOD from predicted Angstrom Exponent (default=%s)"\ + %wavs ) + parser.add_option("-u", "--uncompressed", action="store_true", dest="uncompressed",default=False, help="Do not use n4zip to compress gridded/ODS output file (default=False)") @@ -215,6 +231,8 @@ def makethis_dir(filename): if options.cloudFree == 'None': options.cloudFree = None + elif options.cloudFree == None: + pass else: options.cloudFree = float(options.cloudFree) @@ -223,6 +241,9 @@ def makethis_dir(filename): cloud_thresh=options.cloud_thresh, cloudFree = options.cloudFree, aodmax = options.aodmax, + aodSTD = options.aodSTD, + aodLength = options.aodLength, + wavs = options.wavs, nsyn=options.nsyn, verbose=options.verbose) if modis.nobs < 1: diff --git a/src/Applications/GAAS_App/mxd04_nnr.py b/src/Applications/GAAS_App/mxd04_nnr.py index 1d42ecbc..06cb1efd 100644 --- a/src/Applications/GAAS_App/mxd04_nnr.py +++ b/src/Applications/GAAS_App/mxd04_nnr.py @@ -122,6 +122,7 @@ def __init__(self,l2_path,prod,algo,syn_time,aer_x,slv_x, aodSTD=3.0, aodLength=0.5, coll='006', + wavs=['440','470','550','660','870'], nsyn=8, verbose=0): """ @@ -141,6 +142,7 @@ def __init__(self,l2_path,prod,algo,syn_time,aer_x,slv_x, aodSTD --- number of standard deviations for checking for outliers aodLength --- length scale (degrees) to look for outliers coll --- MODIS data collection + wavs --- wavelengths to calculate output AOD from the Angstrom Exponent nsyn --- number of synoptic times The following attributes are also defined: @@ -160,6 +162,7 @@ def __init__(self,l2_path,prod,algo,syn_time,aer_x,slv_x, self.aodmax = aodmax self.aodSTD = aodSTD self.aodLength = aodLength + self.wavs = wavs # Initialize superclass # --------------------- @@ -486,9 +489,8 @@ def apply(self,nnFile): doAE = True if doAEfit: - wavs = ['440','470','550','660','870'] - wav = np.array(wavs).astype(float) - nwav = len(wavs) + wav = np.array(self.wavs).astype(float) + nwav = len(self.wavs) AEfitb = None for i,targetName in enumerate(self.net.TargetNames): if 'AEfitm' in targetName: @@ -505,7 +507,7 @@ def apply(self,nnFile): targetName = [] for i in range(nwav): targets_[:,i] = -1.*(AEfitm*np.log(wav[i]) + AEfitb) - targetName.append('aTau'+wavs[i]) + targetName.append('aTau'+self.wavs[i]) targets = targets_ self.net.TargetNames = targetName @@ -628,7 +630,8 @@ def apply(self,nnFile): Lat = self.Latitude[self.iGood] gIndex = np.arange(self.nobs)[self.iGood] iOutliers = [] - while find_outliers: + count = 0 + while find_outliers & (count=lon-self.aodLength) & (Lat<=lat+self.aodLength) & (Lat>=lat-self.aodLength) - if np.sum(iHood) <= 1: + if (np.sum(iHood) <= 1) & (maxaod > self.aodmax): #this pixel has no neighbors and is high. Filter it. iOutliers.append(gIndex[imax]) else: @@ -646,6 +649,7 @@ def apply(self,nnFile): iOutliers.append(gIndex[imax]) else: find_outliers = False # done looking for outliers + count +=1 if self.verbose: print("Filtering out ",len(iOutliers)," outlier pixels") diff --git a/src/Applications/GAAS_App/viirs_l2a.pcf b/src/Applications/GAAS_App/viirs_l2a.pcf index 74d3df9c..85ca91ad 100644 --- a/src/Applications/GAAS_App/viirs_l2a.pcf +++ b/src/Applications/GAAS_App/viirs_l2a.pcf @@ -20,7 +20,10 @@ VIIRS_L2A_NN_FILE = nnr_001.%ident_TauAE.net VIIRS_L2A_BLANK_ODS = /nobackup/NNR/Misc/blank.ods VIIRS_L2A_CLOUD_THRESH = 0.7 -VIIRS_L2A_CLOUDFREE = 0.1 -VIIRS_L2A_AODMAX = 10 +VIIRS_L2A_CLOUDFREE = 0.5 +VIIRS_L2A_AODMAX = 2.0 +VIIRS_L2A_AODSTD = 3.0 +VIIRS_L2A_AODLENGTH = 0.5 +VIIRS_L2A_WAVS = 440,470,550,660,870 #END VIIRS_L2A processing diff --git a/src/Applications/GAAS_App/viirs_l2a.py b/src/Applications/GAAS_App/viirs_l2a.py index deb18767..56ae85f4 100755 --- a/src/Applications/GAAS_App/viirs_l2a.py +++ b/src/Applications/GAAS_App/viirs_l2a.py @@ -50,10 +50,10 @@ if cf('VIIRS_L2A_OVERWRITE').upper() == 'YES': Options += " --force" if cf('VIIRS_L2A_VERBOSE').upper() == 'YES': Options += " -v" - + try: - aodmax = cf('VIIRS_L2A_AODMAX') - Options += " --aodmax=" + aodmax + cloud_thresh = cf('VIIRS_L2A_CLOUD_THRESH') + Options += " --cloud_thresh=" + cloud_thresh except: pass @@ -64,11 +64,29 @@ pass try: - cloud_thresh = cf('VIIRS_L2A_CLOUD_THRESH') - Options += " --cloud_thresh=" + cloud_thresh + aodmax = cf('VIIRS_L2A_AODMAX') + Options += " --aodmax=" + aodmax except: pass + try: + aodSTD = cf('VIIRS_L2A_AODSTD') + Options += " --aodSTD=" + aodSTD + except: + pass + + try: + aodLenth = cf('VIIRS_L2A_AODLENGTH') + Options += " --aodLength=" + aodLength + except: + pass + + try: + wavs = cf('VIIRS_L2A_WAVS') + Options += " --wavs=" + wavs + except: + pass + # Generate products # ----------------- i = 0 diff --git a/src/Applications/GAAS_App/vx04_l2a.py b/src/Applications/GAAS_App/vx04_l2a.py index 80441392..acf4170a 100755 --- a/src/Applications/GAAS_App/vx04_l2a.py +++ b/src/Applications/GAAS_App/vx04_l2a.py @@ -42,7 +42,8 @@ def makethis_dir(filename): raise IOError("could not create directory "+path) #--------------------------------------------------------------------- - +def wavs_callback(option, opt, value, parser): + setattr(parser.values, option.dest, value.split(',')) if __name__ == "__main__": expid = 'nnr_001' @@ -70,6 +71,9 @@ def makethis_dir(filename): cloud_thresh = 0.70 cloudFree = None aodmax = 1.0 + aodSTD = 3.0 + aodLength = 0.5 + wavs = '440,470,550,660,870' # Parse command line options # -------------------------- @@ -129,6 +133,18 @@ def makethis_dir(filename): help="max AOD value that will be accepted when cloud fraction is greater than cloudFree (default=%f)"\ %aodmax ) + parser.add_option("--aodSTD", dest="aodSTD", default=aodSTD,type='float', + help="number of standard deviations to check for AOD outliers (default=%f)"\ + %aodSTD ) + + parser.add_option("--aodLength", dest="aodLength", default=aodLength,type='float', + help="length scale (degrees) to check for AOD outliers (default=%f)"\ + %aodLength ) + + parser.add_option("--wavs", dest="wavs", default=wavs,type='string',action='callback',callback=wavs_callback, + help="wavelength to output AOD from predicted Angstrom Exponent (default=%s)"\ + %wavs ) + parser.add_option("-u", "--uncompressed", action="store_true", dest="uncompressed",default=False, help="Do not use n4zip to compress gridded/ODS output file (default=False)") @@ -205,6 +221,8 @@ def makethis_dir(filename): if options.cloudFree == 'None': options.cloudFree = None + elif options.cloudFree == None: + pass else: options.cloudFree = float(options.cloudFree) @@ -213,6 +231,9 @@ def makethis_dir(filename): cloud_thresh=options.cloud_thresh, cloudFree=options.cloudFree, aodmax=options.aodmax, + aodSTD = options.aodSTD, + aodLength = options.aodLength, + wavs = options.wavs, nsyn=options.nsyn, verbose=options.verbose) if viirs.nobs < 1: diff --git a/src/Applications/GAAS_App/vx04_nnr.py b/src/Applications/GAAS_App/vx04_nnr.py index f7191dd0..94def230 100644 --- a/src/Applications/GAAS_App/vx04_nnr.py +++ b/src/Applications/GAAS_App/vx04_nnr.py @@ -61,7 +61,10 @@ def __init__(self,l2_path,sat,algo,syn_time,aer_x, glint_thresh=40.0, scat_thresh=170.0, cloudFree=None, - aodmax=1.0, + aodmax=2.0, + aodSTD=3.0, + aodLength=0.5, + wavs=['440','470','550','660','870'], coll='002', nsyn=8, verbose=0): @@ -82,6 +85,8 @@ def __init__(self,l2_path,sat,algo,syn_time,aer_x, cloud_tresh --- cloud fraction threshhold cloudFree --- cloud fraction threshhold for assuring no cloud contaminations when aod is > aodmax if None, no cloud free check is made + aodSTD --- number of standard deviations for checking for outliers + aodLength --- length scale (degrees) to look for outliers coll --- VIIRS data collection nsyn --- number of synoptic times @@ -98,6 +103,9 @@ def __init__(self,l2_path,sat,algo,syn_time,aer_x, self.algo = algo self.cloudFree = cloudFree self.aodmax = aodmax + self.aodSTD = aodSTD + self.aodLength = aodLength + self.wavs = wavs # Initialize superclass # set anet_wav to True so MODIS wavelengths align with AERONET @@ -335,9 +343,8 @@ def apply(self,nnFile): doAE = True if doAEfit: - wavs = ['440','470','550','660','870'] - wav = np.array(wavs).astype(float) - nwav = len(wavs) + wav = np.array(self.wavs).astype(float) + nwav = len(self.wavs) AEfitb = None for i,targetName in enumerate(self.net.TargetNames): if 'AEfitm' in targetName: @@ -354,7 +361,7 @@ def apply(self,nnFile): targetName = [] for i in range(nwav): targets_[:,i] = -1.*(AEfitm*np.log(wav[i]) + AEfitb) - targetName.append('aTau'+wavs[i]) + targetName.append('aTau'+self.wavs[i]) targets = targets_ self.net.TargetNames = targetName @@ -368,11 +375,14 @@ def apply(self,nnFile): I = np.array(self.channels) < 900 # only visible channels, this is relevant for ocean aechannels = np.array(self.channels)[I] aodT = self.aod[:,I].T - fit = np.polyfit(np.log(aechannels),-1.*np.log(aodT[:,self.iGood]+0.01),1) + iIndex = np.arange(len(self.iGood))[self.iGood] + aodT = aodT[:,iIndex] + 0.01 + mask = aodT.min(axis=0) > 0 + posIndex = iIndex[mask] + fit = np.polyfit(np.log(aechannels),-1.*np.log(aodT[:,mask]+0.01),1) self.ae = MISSING*np.ones(self.nobs) - self.ae[self.iGood] = fit[0,:] - bad = np.isnan(self.ae) - self.ae[bad] = MISSING + self.ae[posIndex] = fit[0,:] + if doAE: for i,targetName in enumerate(self.net.TargetNames): @@ -426,7 +436,8 @@ def apply(self,nnFile): # Do extra cloud filtering if required - if self.cloudFree is not None: + if self.cloudFree is not None: + # start by checking the cloud masks cloudy = (self.cloud>=self.cloudFree) contaminated = np.zeros(np.sum(self.iGood)).astype(bool) @@ -437,18 +448,70 @@ def apply(self,nnFile): contaminated = contaminated | ( (result > self.aodmax) & cloudy[self.iGood] ) icontaminated = np.arange(self.nobs)[self.iGood][contaminated] + + if self.verbose: + print('Filtering out ',np.sum(contaminated),' suspected cloud contaminated pixels') + + for targetName in self.net.TargetNames: name, ch = TranslateTarget(targetName) k = list(self.channels).index(ch) # index of channel self.__dict__[name][icontaminated,k] = MISSING if doAEfit: - self.ae[icontaminated] = MISSING self.ae_[icontaminated] = MISSING self.iGood[icontaminated] = False + # check for outliers + # start with highest AOD550 value + # find all the pixels within a 1 degree neighborhood + # check if it is outside of mean + N*sigma of the other pixels + # aodSTD parameter is equal to N + # continue until no outliers are found + find_outliers = True + k = list(self.channels).index(550) + aod550 = np.ma.array(self.aod_[self.iGood,k]) + aod550.mask = np.zeros(len(aod550)).astype(bool) + Lon = self.Longitude[self.iGood] + Lat = self.Latitude[self.iGood] + gIndex = np.arange(self.nobs)[self.iGood] + iOutliers = [] + count = 0 + while find_outliers & (count=lon-self.aodLength) & (Lat<=lat+self.aodLength) & (Lat>=lat-self.aodLength) + if (np.sum(iHood) <= 1) & (maxaod > self.aodmax): + #this pixel has no neighbors and is high. Filter it. + iOutliers.append(gIndex[imax]) + else: + aodHood = aod550[iHood] + if maxaod > (aodHood.mean() + self.aodSTD*aodHood.std()): + iOutliers.append(gIndex[imax]) + else: + find_outliers = False # done looking for outliers + count +=1 + if self.verbose: + print("Filtering out ",len(iOutliers)," outlier pixels") + + self.iOutliers = iOutliers + + if len(iOutliers) > 0: + for targetName in self.net.TargetNames: + name, ch = TranslateTarget(targetName) + k = list(self.channels).index(ch) # index of channel + self.__dict__[name][iOutliers,k] = MISSING + + if doAEfit: + self.ae_[iOutliers] = MISSING + self.iGood[iOutliers] = False