Feature/pcastell/cloud filter inputs

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
 

src/Applications/GAAS_App/modis_l2a.pcf

@@ -22,4 +22,11 @@ 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.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

src/Applications/GAAS_App/modis_l2a.py

@@ -47,7 +47,44 @@
 
     if cf('MODIS_L2A_OVERWRITE').upper() == 'YES': Options += " --force"
     if   cf('MODIS_L2A_VERBOSE').upper() == 'YES': Options += " -v"
-
+
+
+    try:
+        cloud_thresh = cf('MODIS_L2A_CLOUD_THRESH')
+        Options += " --cloud_thresh=" + cloud_thresh
+    except:
+        pass
+
+    try:
+        cloudFree = cf('MODIS_L2A_CLOUDFREE')
+        Options += " --cloudFree=" + cloudFree
+    except:
+        pass
+
+    try:
+        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

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,6 +72,12 @@ def makethis_dir(filename):
     coll = '006'
     res = 'c'
     nsyn = 8
+    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
 #   --------------------------
@@ -120,8 +127,32 @@ 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("--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)")
@@ -198,11 +229,21 @@ 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
+    elif options.cloudFree == None:
+        pass
+    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,
+                      aodSTD = options.aodSTD,
+                      aodLength = options.aodLength,
+                      wavs = options.wavs,
                       nsyn=options.nsyn,
                       verbose=options.verbose)
     if modis.nobs < 1:

src/Applications/GAAS_App/mxd04_nnr.py

@@ -118,8 +118,11 @@ 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',
+                 wavs=['440','470','550','660','870'],
                  nsyn=8,
                  verbose=0):
         """
@@ -136,7 +139,10 @@ 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
+        wavs        --- wavelengths to calculate output AOD from the Angstrom Exponent
         nsyn         --- number of synoptic times              
 
         The following attributes are also defined:
@@ -154,7 +160,10 @@ 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
+        self.wavs = wavs
+
         # Initialize superclass
         # ---------------------
         Files = granules(l2_path,prod,syn_time,coll=coll,nsyn=nsyn)
@@ -480,20 +489,25 @@ 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:
                         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 = []
             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
@@ -507,11 +521,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:
@@ -567,32 +583,89 @@ 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":
                 cloudy = (self.cloud_lnd>=self.cloudFree) & (self.cloud>=self.cloudFree)
             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] )
 
+
             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<len(aod550)):
+                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) & (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                
+
 
 #---
 

src/Applications/GAAS_App/viirs_l2a.pcf

@@ -19,4 +19,11 @@ 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.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

src/Applications/GAAS_App/viirs_l2a.py

@@ -46,11 +46,47 @@
                   " --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:
+        cloud_thresh = cf('VIIRS_L2A_CLOUD_THRESH')
+        Options += " --cloud_thresh=" + cloud_thresh
+    except:
+        pass
+
+    try:
+        cloudFree = cf('VIIRS_L2A_CLOUDFREE')
+        Options += " --cloudFree=" + cloudFree
+    except:
+        pass
+
+    try:
+        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