From bca3ab111cb4389bd95c9c41e3599d690a7d0692 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Wed, 28 Sep 2022 19:11:05 -0400 Subject: [PATCH 01/67] PC - working version of VIIRS reader vx04.py --- .../GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 938 ++++++++++++++++++ 1 file changed, 938 insertions(+) create mode 100644 src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py new file mode 100644 index 00000000..d9820f3c --- /dev/null +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -0,0 +1,938 @@ +""" +Reads Level 2 Deep Blue and Dark Target retrievals applied to VIIRS +granules for a single day and returns a +single object with the relevant data. + +Builds on the heritage of mxd04.py + +This software is hereby placed in the public domain. +patricia.castellanos@nasa.gov +""" + +import os +import sys +import numpy as np +from datetime import date, datetime, timedelta +from glob import glob +from dateutil.parser import isoparse +from pyobs.npz import NPZ +from binObs_ import binobs2d, binobs3d + +try: + from pyods import ODS # must be inported before pyhdf +except: + pass + +from netCDF4 import Dataset +from bits import BITS + +#--- + +DATE_START = datetime(1993,1,1,0,0,0) + +SDS = dict ( + DT_META = ('longitude', 'latitude', + 'sensor_zenith_angle', 'sensor_azimuth_angle', + 'solar_zenith_angle', 'solar_azimuth_angle', + 'Scattering_Angle', 'Glint_Angle'), + DT_LAND = ( 'Mean_Reflectance_Land', + 'Corrected_Optical_Depth_Land', + 'Surface_Reflectance_Land', + 'Aerosol_Cloud_Fraction_Land', + 'Land_Ocean_Quality_Flag', + 'Number_Pixels_Used_Land', + 'STD_Reflectance_Land'), + DT_OCEAN = ( 'Effective_Optical_Depth_Average_Ocean', + 'Optical_Depth_Small_Average_Ocean', + 'Effective_Radius_Ocean', + 'Asymmetry_Factor_Average_Ocean', + 'Angstrom_Exponent_1_Ocean', + 'Angstrom_Exponent_2_Ocean', + 'Aerosol_Cloud_Fraction_Ocean', + 'Mean_Reflectance_Ocean', + 'Land_Ocean_Quality_Flag' ), + DB_META = ('Longitude', 'Latitude', 'Scan_Start_Time', + 'Viewing_Zenith_Angle', 'Relative_Azimuth_Angle', + 'Solar_Zenith_Angle', + 'Scattering_Angle' ), + DB_LAND = ('Aerosol_Optical_Thickness_550_Land', + 'Spectral_Aerosol_Optical_Thickness_Land', + 'Spectral_Single_Scattering_Albedo_Land', + 'Spectral_Surface_Reflectance', + 'Spectral_TOA_Reflectance_Land', + 'Aerosol_Optical_Thickness_550_Land_Best_Estimate', + 'Aerosol_Optical_Thickness_550_STDV_Land', + 'Aerosol_Optical_Thickness_QA_Flag_Land', + 'Algorithm_Flag_Land', + 'Aerosol_Type_Land', + 'Number_Of_Pixels_Used_Land', + 'Number_Valid_Pixels'), + DB_OCEAN = ('Aerosol_Optical_Thickness_550_Ocean', + 'Spectral_Aerosol_Optical_Thickness_Ocean', + 'Spectral_TOA_Reflectance_Ocean', + 'Aerosol_Optical_Thickness_550_Ocean_Best_Estimate', + 'Aerosol_Optical_Thickness_550_STDV_Ocean', + 'Aerosol_Optical_Thickness_QA_Flag_Ocean', + 'Algorithm_Flag_Ocean', + 'Aerosol_Type_Ocean', + 'Number_Of_Pixels_Used_Ocean', + 'Number_Valid_Pixels', + 'Wind_Speed', + 'Wind_Direction') + ) +# NOTE: DEEP BLUE does not have cloud information in their files. + + + +# AOD Shannels +CHANNELS = dict ( + DT_LAND = ( 480., 550., 670., 2250.), + DT_OCEAN = ( 480., 550., 670., 860., 1240., 1600., 2250. ), + DB_LAND = ( 412, 488, 550, 670 ), # file has 550 separate + DB_OCEAN = (488., 550., 670., 865., 1240., 1640., 2250.), + SREF = ( 480., 670., 2250. ) + ) + + + +ALIAS = dict ( Longitude = 'lon', + Latitude = 'lat', + longitude = 'lon', + latitude = 'lat', + Viewing_Zenith_Angle = 'SensorZenith', + sensor_zenith_angle = 'SensorZenith', + sensor_azimuth_angle = 'SensorAzimuth', + Relative_Azimuth_Angle = 'RelativeAzimuth', + Solar_Zenith_Angle = 'SolarZenith', + solar_zenith_angle = 'SolarZenith', + solar_azimuth_angle = 'SolarAzimuth', + Scattering_Angle = 'ScatteringAngle', + Glint_Angle = 'GlintAngle', + Mean_Reflectance_Land = 'reflectance', + Surface_Reflectance_Land = 'sfc_reflectance', + Corrected_Optical_Depth_Land = 'aod', + Aerosol_Cloud_Fraction_Land = 'cloud', + Effective_Optical_Depth_Average_Ocean = 'aod', + Optical_Depth_Small_Average_Ocean = 'aod_fine', + Aerosol_Cloud_Fraction_Ocean = 'cloud', + Mean_Reflectance_Ocean = 'reflectance', + Spectral_Aerosol_Optical_Thickness_Land = 'aod3ch', + Aerosol_Optical_Thickness_550_Land = 'aod550', + Spectral_Surface_Reflectance = 'sfc_reflectance', + Spectral_TOA_Reflectance_Land = 'reflectance', + Spectral_Single_Scattering_Albedo_Land = 'ssa', + Algorithm_Flag_Land = 'atype', + Angstrom_Exponent_Land = 'angstrom', + Spectral_Aerosol_Optical_Thickness_Ocean = 'aod', + Aerosol_Optical_Thickness_550_Ocean = 'aod550', + Spectral_TOA_Reflectance_Ocean = 'reflectance', + Algorithm_Flag_Ocean = 'atype', + Angstrom_Exponent_Ocean = 'angstrom', + Number_Of_Pixels_Used_Ocean = 'Number_Of_Pixels_Used', + Number_Of_Pixels_Used_Land = 'Number_Of_Pixels_Used', + Aerosol_Optical_Thickness_QA_Flag_Land = 'qa_flag', + Aerosol_Optical_Thickness_QA_Flag_Ocean = 'qa_flag', + Land_Ocean_Quality_Flag = 'qa_flag', + Scan_Start_Time = 'Time', + ) + +BAD, MARGINAL, GOOD, BEST = ( 0, 1, 2, 3 ) # DT QA marks +# for DB 0 = no retrieval, 1 = poor, 2 = moderate, 3 = good + +translate_sat = {'Suomi-NPP': 'SNPP'} + + +KX = dict ( SNPP_DT_OCEAN = 301, + SNPP_DT_LAND = 302, + SNPP_DB_OCEAN = 310, + SNPP_DB_LAND = 311, + ) + +KT = dict ( AOD = 45, ) + +IDENT = dict ( SNPP_DT_OCEAN = 'vsnppdto', + SNPP_DT_LAND = 'vsnppdtl', + SNPP_DB_OCEAN = 'vsnppdbo', + SNPP_DB_LAND = 'vsnppdbl', + ) + +MISSING = 999.999 + +#........................................................................... + +class Vx04_L2(object): + """ + This class implements the VIIRS Level 2 AEROSOL products. + VIIRS currently flies on the Suomi NPP (SNPP) and NOAA-20 (FKA JPSS-1) satellites + VIIRS will be flown on JPSS-2,3, and 4 + The Level2 AEROSOL products are implementations of the MODIS Deep Blue and Dark Target algoritms, + referred to as MOD04 (TERRA satellite) and MYD04 (AQUA satellite). + Therefore, we will refer to these products at VSNPP04 and VN2004, respectively + """ + + def __init__ (self,Path,Algo,Surface,syn_time=None,nsyn=8,Verb=0, + only_good=True,SDS=SDS,alias=None): + """ + Reads individual granules or a full day of Level 2 Vx04 files + present on a given *Path* and returns a single object with + all data concatenated for a given algorithm. On input, + + Required parameters: + Path -- can be a single file, a single directory, of a list + of files and directories. Directories are + transversed recursively. If a non Vx04 Level 2 + file is encountered, it is simply ignored. + Algo -- Algorithm: DT or DB + Surface -- OCEAN or LAND + + Optional parameters: + syn_type --- synoptic time + nsyn --- number of synoptic times per day + only_good --- keep only *good* observations + Verb -- Verbose level: + 0 - really quiet (default) + 1 - Warns if invalid file is found + 2 - Prints out non-zero number of aerosols in each file. + SDS --- Variables to be read from L2 Aerosol files. Must + be a dictionary with keys '{Algo}_META' and '{Algo}_{Surface}' + ALIAS --- dictionary of alises for SDSs + + """ + + algo = '{}_{}'.format(Algo,Surface) + if algo not in ('DT_LAND', 'DT_OCEAN', 'DB_LAND', 'DB_OCEAN'): + raise ValueError, "invalid algorithm "+algo+" --- must be DT_LAND, DT_OCEAN, DB_LAND, DB_OCEAN" + +# Initially are lists of numpy arrays for each granule +# ------------------------------------------------ + self.verb = Verb + self.sat = None # Satellite name + self.col = None # collection, e.g., 011 + self.algo = algo + self.SDS = SDS['{}_META'.format(Algo)] + SDS[algo] + self.SDS_META = SDS['{}_META'.format(Algo)] + + # Add/Substitute some aliases if given + # ------------------------------------ + self.ALIAS = ALIAS.copy() + if alias is not None: + for a in alias: self.ALIAS[a] = alias[a] + + + # Create empty lists for SDS to be read from file + # ----------------------------------------------- + for name in self.SDS: + self.__dict__[name] = [] + + # Read each granule, appending them to the list + # --------------------------------------------- + if type(Path) is list: + if len(Path) == 0: + self.nobs = 0 + print "WARNING: Empty Vx04_L2 object created" + return + else: + Path = [Path, ] + self._readList(Path) + + #Protect against empty MXD04 files + # -------------------------------- + if len(self.Scattering_Angle) == 0: + self.nobs = 0 + print "WARNING: Empty MxD04_L2 object created" + return + + # Make each attribute a single numpy array + # ---------------------------------------- + if 'DT' in self.algo: + self.SDS += ('Scan_Start_Time',) + for sds in self.SDS: + try: + self.__dict__[sds] = np.ma.concatenate(self.__dict__[sds]) + except: + print "Failed concatenating "+sds + + + # Determine index of "good" observations + # -------------------------------------- + if self.algo == 'DT_LAND': + self.iGood = self.Land_Ocean_Quality_Flag == BEST + elif self.algo == 'DT_OCEAN': + self.iGood = self.Land_Ocean_Quality_Flag > BAD + elif self.algo == 'DB_LAND': + self.iGood = self.Aerosol_Optical_Thickness_QA_Flag_Land > BAD # for now + elif self.algo == 'DB_OCEAN': + self.iGood = self.Aerosol_Optical_Thickness_QA_Flag_Ocean > BAD + else: + raise ValueError, 'invalid algorithm (very strange)' + + + # Keep only "good" observations + # ----------------------------- + if only_good: + m = self.iGood + for sds in self.SDS: + rank = len(self.__dict__[sds].shape) + if rank == 1: + self.__dict__[sds] = self.__dict__[sds][m] + elif rank == 2: + self.__dict__[sds] = self.__dict__[sds][m,:] + else: + raise IndexError, 'invalid rank=%d'%rank + self.iGood = self.iGood[m] + + + # Make aliases for compatibility with older code + # ---------------------------------------------- + Alias = self.ALIAS.keys() + for sds in self.SDS: + if sds in Alias: + self.__dict__[self.ALIAS[sds]] = self.__dict__[sds] + + # Calculate Glint angle for Deep Blue + # ---------------------- + if 'DB' in self.algo: + sza = np.radians(self.SolarZenith) + vza = np.radians(self.SensorZenith) + raa = np.radians(self.RelativeAzimuth) + cglint = np.cos(sza)*np.cos(vza) + np.sin(sza)*np.sin(vza)*np.cos(np.pi-raa) + self.GlintAngle = np.arccos(cglint) + + + # Create corresponding python time + # -------------------------------- + if 'DB' in self.algo: + self.Time = np.array([DATE_START+timedelta(seconds=s) for s in self.Scan_Start_Time]) + + # ODS friendly attributes + # ----------------------- + self.nobs = self.Scattering_Angle.shape[0] + self.kx = KX[self.sat+'_'+self.algo] + self.ident = IDENT[self.sat+'_'+self.algo] + self.channels = CHANNELS[self.algo] + if Algo == 'DT': + self.sChannels = CHANNELS["SREF"] # LAND surface reflectivity (not the same as algo) + else: + self.sChannels = CHANNELS[self.algo] # deep blue surface bands are the same as algorithm + + if 'DB' in self.algo: + self.rChannels = self.Reflectance_Bands + elif self.algo == 'DT_LAND': + self.rChannels = np.array([480.,670.,2250.]) + elif self.algo == 'DT_OCEAN': + self.rChannels = np.array([480.,550.,670.,860.,1240.,1600.,2250.]) + + + if syn_time == None: + self.syn_time = None + self.time = None + self.nymd = None + self.nhms = None + self.nsyn = None + else: + Dt = [ t-syn_time for t in self.Time ] + self.time = np.array([ (86400.*dt.days+dt.seconds)/60. for dt in Dt ]) # in minutes + self.syn_time = syn_time + self.nymd = 10000 * syn_time.year + 100*syn_time.month + syn_time.day + self.nhms = 10000 * syn_time.hour + 100*syn_time.minute + syn_time.second + self.nsyn = nsyn # number of synoptic times per day + + # Concatenate AOD channels for Deep Blue + # -------------------------------------- + if self.algo == 'DB_LAND': + try: + self.aod = np.ones((self.nobs,4)) + self.aod[:,0] = self.aod3ch[:,0] + self.aod[:,1] = self.aod3ch[:,1] + self.aod[:,2] = self.aod550[:] + self.aod[:,3] = self.aod3ch[:,2] + except: + pass # don't fuss, aod3ch may not have been read + + # Create a pseudo cloud fraction for Deep Blue + if Algo == 'DB': + self.cloud = self.Number_Of_Pixels_Used.astype(float)/self.Number_Valid_Pixels.astype(float) + +#--- + def _readList(self,List): + """ + Recursively, look for files in list; list items can + be files or directories. + """ + for item in List: + if os.path.isdir(item): self._readDir(item) + elif os.path.isfile(item): + if 'DB' in self.algo: + self._readGranuleDB(item) + else: + self._readGranuleDT(item) + else: + print "%s is not a valid file or directory, ignoring it"%item +#--- + def _readDir(self,dir): + """Recursively, look for files in directory.""" + for item in os.listdir(dir): + path = dir + os.sep + item + if os.path.isdir(path): self._readDir(path) + elif os.path.isfile(path): self._readGranule(path) + else: + print "%s is not a valid file or directory, ignoring it"%item + +#--- + def _readGranuleDB(self,filename): + """Reads one Vx04 Deep Blue granule with Level 2 aerosol data.""" + + # Don't fuss if the file cannot be opened + # --------------------------------------- + try: + if self.verb: + print "[] Working on "+filename + nc = Dataset(filename) + except: + if self.verb > 2: + print "- %s: not recognized as an netCDF file"%filename + return + + # Read select variables (reshape to allow concatenation later) + # ------------------------------------------------------------ + for sds in self.SDS: + v = nc.variables[sds][:] + a = nc.variables[sds].ncattrs() + if 'scale_factor' in a: + scale = nc.variables[sds].getncattr('scale_factor') + v = scale*v + if 'add_offset' in a: + add = nc.variables[sds].getncattr('add_offset') + v = v + add + + if len(v.shape) == 3: + if "TOA_Reflectance" in sds: + i, j, k = v.shape + v = v.reshape((i*j,k)) + else: + i, j, k = v.shape + v = v.reshape((i,j*k)).T + elif len(v.shape) == 2: + v = v.ravel() + else: + raise IndexError, "invalid shape for SDS <%s>"%sds + self.__dict__[sds].append(v) + + +# Satellite name +# -------------- + if self.sat is None: + self.sat = translate_sat[nc.platform] + +# Collection +# ---------- + if self.col is None: + self.col = nc.product_name.split('.')[-3] + +# Reflectance Bands +# ------------------ + self.Reflectance_Bands = nc.variables['Reflectance_Bands'][:] + +#--- + def _readGranuleDT(self,filename): + """Reads one Vx04 Dark Target granule with Level 2 aerosol data.""" + + # Don't fuss if the file cannot be opened + # --------------------------------------- + try: + if self.verb: + print "[] Working on "+filename + nc = Dataset(filename) + data = nc.groups['geophysical_data'] + loc = nc.groups['geolocation_data'] + except: + if self.verb > 2: + print "- %s: not recognized as an netCDF file"%filename + return + + # Read select variables (reshape to allow concatenation later) + # ------------------------------------------------------------ + for sds in self.SDS: + if sds in self.SDS_META: + vobj = loc.variables[sds] + else: + vobj = data.variables[sds] + + a = vobj.ncattrs() + v = vobj[:] + + # commenting this out for now + # all the scale factors in the DT files seem to be wrong + # maybe it will be fixed in future versions, so keeping + # code for now +# if 'scale_factor' in a: +# scale = vobj.getncattr('scale_factor') +# v = scale*v +# if 'add_offset' in a: +# add = vobj.getncattr('add_offset') +# v = v + add + + if len(v.shape) == 3: + i, j, k = v.shape + v = v.reshape((i*j,k)) + elif len(v.shape) == 2: + v = v.ravel() + else: + raise IndexError, "invalid shape for SDS <%s>"%sds + self.__dict__[sds].append(v) + + +# Satellite name +# -------------- + if self.sat is None: + self.sat = translate_sat[nc.platform] + +# Collection +# ---------- + if self.col is None: + self.col = nc.product_name.split('.')[-3] + +# Time + if 'Scan_Start_Time' not in self.__dict__: + self.Scan_Start_Time = [np.repeat(isoparse(nc.time_coverage_start[:-1]),i*j)] + else: + self.Scan_Start_Time.append(np.repeat(isoparse(nc.time_coverage_start[:-1]),i*j)) + +#--- + + def reduce(self,I): + """ + Reduce observations according to index I. + """ + Nicknames = self.ALIAS.values() + for name in self.__dict__: + if name in Nicknames: + continue # alias do not get reduced + q = self.__dict__[name] + if type(q) is type(self.lon): + if len(q) == self.nobs: + # print "{} Reducing "+name + self.__dict__[name] = q[I] + + Alias = self.ALIAS.keys() + for sds in self.SDS: + if sds in Alias: + self.__dict__[self.ALIAS[sds]] = self.__dict__[sds] # redefine aliases + + self.nobs = len(self.lon) + +#--- + def write(self,filename=None,dir='.',expid=None,Verb=1): + """ + Writes the un-gridded data to a numpy npz file. + """ + + # Stop here is no good obs available + # ---------------------------------- + if self.nobs == 0: + return # no data to work with + if any(self.iGood) == False: + return # no good data to work with + + if expid == None: + expid = self.algo + + if filename is None: + filename = '%s/%s.viirs.%d_%02dz.npz'%(dir,expid,self.nymd,self.nhms/10000) + + version = 1 # File format version + meta = [self.nymd,self.nhms,self.nobs,self.nch,self.kx,version] + savez(filename, + meta = meta, + lon = self.lon, + lat = self.lat, + ks = self.ks, + channels = self.channels, + qa_flag = self.qa_flag, + SolarZenith = self.SolarZenith, + SensorZenith = self.SensorZenith, + RelativeAzimuth = self.RelativeAzimuth, + ScatteringAngle = self.ScatteringAngle, + cloud = self.cloud, + aod = self.aod, + reflectance = self.reflectance) + + if Verb >=1: + print "[w] Wrote file "+filename +#--- + + def writeODS(self,filename=None,dir='.',expid=None,channels=None, + revised=False,nsyn=8,Verb=1): + """ + Writes the un-gridded data to an ODS file. If *revised* + is True, the revised *aod_* parameter is written to file. + """ + + if self.syn_time == None: + raise ValuError, "synoptic time missing, cannot write ODS" + + # Stop here if no good obs available + # ---------------------------------- + if self.nobs == 0: + return # no data to work with + if any(self.iGood) == False: + return # no good data to work with + + if expid == None: + expid = self.algo + + if filename is None: + filename = '%s/%s.obs.%d_%02dz.ods'%(dir,expid,self.nymd,self.nhms/10000) + + if channels is None: + channels = self.channels + + # Create and populated ODS object + # ------------------------------- + ns = self.nobs + nobs = len(channels) * ns + ods = ODS(nobs=nobs, kx=self.kx, kt=KT['AOD']) + i = 0 + ks = np.arange(ns) + 1 + for ch in channels: + I = range(i,i+ns) + j = channels.index(ch) + ods.ks[I] = ks + ods.lat[I] = self.lat[:] + ods.lon[I] = self.lon[:] + ods.time[I] = self.time[:].astype('int') + ods.lev[I] = ch + ods.qch[I] = self.qa_flag[:].astype('int') + if revised: + ods.obs[I] = self.aod_[:,j].astype('float32') + ods.xvec[I] = self.aod[:,j].astype('float32') + else: + ods.obs[I] = self.aod[:,j].astype('float32') + i += ns + + # Handle corrupted coordinates + # ---------------------------- + iBad = (ods.lon<-180) | (ods.lon>180.) | \ + (ods.lat<-90) | (ods.lat>90.) | \ + (abs(ods.time)>1440./self.nsyn) + ods.lon[iBad] = 0. + ods.lat[iBad] = -90. + ods.time[iBad] = 0. + ods.qcx[iBad] = 2 + + # Exclusion flag + # -------------- + iGood = (ods.qch>0) & (ods.obs<10.) & (ods.qcx==0) + ods.qcx[:] = 1 # All bad... + ods.qcx[iGood] = 0 # ... unless good + + ods_ = ods.select(qcx=0) + if Verb >=1: + print "[w] Writing file <"+filename+"> with %d observations"%ods_.nobs + + ods_.write(filename,self.nymd,self.nhms,nsyn=8,ftype='pre_anal') + +#--- + def writeg(self,filename=None,dir='.',expid=None,refine=8,res=None, + channels=None,Verb=1): + """ + Writes gridded MODIS measurements to file. + + refine -- refinement level for a base 4x5 GEOS-5 grid + refine=1 produces a 4 x 5 grid + refine=2 produces a 2 x2.50 grid + refine=4 produces a 1 x1,25 grid + refine=8 produces a 0.50x0.625 grid + refine=16 produces a 0.25x0.3125 grid + Alternatively, one can specify the grid resolution with a + single letter: + + res -- single letter denoting GEOS-5 resolution, + res='a' produces a 4 x 5 grid + res='b' produces a 2 x2.50 grid + res='c' produces a 1 x1,25 grid + res='d' produces a 0.50x0.625 grid + res='e' produces a 0.25x0.3125 grid + + NOTE: *res*, if specified, supersedes *refine*. + + Verb -- Verbose level: + 0 - really quiet (default) + 1 - Warns if invalid file is found + 2 - Prints out non-zero number of fires in each file. + + + """ + from gfio import GFIO + + # Stop here is no good obs available + # ---------------------------------- + if self.nobs == 0: + return # no data to work with + if any(self.iGood) == False: + return # no good data to work with + + if expid == None: + expid = self.algo + +# Output grid resolution +# ---------------------- + if res is not None: + if res=='a': refine = 1 + if res=='b': refine = 2 + if res=='c': refine = 4 + if res=='d': refine = 8 + if res=='e': refine = 16 + +# Lat lon grid +# ------------ + dx = 5. / refine + dy = 4. / refine + im = int(360. / dx) + jm = int(180. / dy + 1) + + glon = np.linspace(-180.,180.,im,endpoint=False) + glat = np.linspace(-90.,90.,jm) + + if channels is None: + channels = self.channels + + levs = np.array(channels) + + nch = len(channels) + nymd = self.nymd + nhms = self.nhms + + vtitle = [ 'Aerosol Optical Depth', + 'Aerosol Optical Depth (Revised)', + 'Aerosol Optical Depth (Fine Mode)', + 'Cloud Fraction' ] + + vname = ['tau', 'tau_', 'tau_fine', 'cloud' ] + vunits = [ '1', '1', '1', '1', ] + kmvar = [ nch, nch, nch, 0 ] + + title = 'Gridded MODIS Aerosol Retrievals' + source = 'NASA/GSFC/GMAO GEOS-5 Aerosol Group' + contact = 'arlindo.dasilva@nasa.gov' + + if filename is None: + filename = '%s/%s.sfc.%d_%02dz.nc4'%(dir,expid,self.nymd,self.nhms/10000) + + # Create the file + # --------------- + f = GFIO() + f.create(filename, vname, nymd, nhms, + lon=glon, lat=glat, levs=levs, levunits='nm', + vtitle=vtitle, vunits=vunits,kmvar=kmvar,amiss=MISSING, + title=title, source=source, contact=contact) + + # Subset AOD at specified channels + # -------------------------------- + I = [] + for ch in channels: + i = list(self.channels).index(ch) + I = I + [i,] + aod = self.aod[:,I] + + # Fine mode + # --------- + try: + aod_fine = self.aod_fine[:,I] + except: + aod_fine = MISSING * np.ones(aod.shape) # will compress like a charm + + # The Revised AOD may not exist + # ------------------------------- + try: + aod_ = self.aod_[:,I] + except: + aod_ = MISSING * np.ones(aod.shape) # will compress like a charm + + # Grid variable and write to file + # ------------------------------- + f.write('tau', nymd, nhms, + binobs3d(self.lon,self.lat,aod,im,jm,MISSING) ) + f.write('tau_', nymd, nhms, + binobs3d(self.lon,self.lat,aod_,im,jm,MISSING) ) + f.write('tau_fine', nymd, nhms, + binobs3d(self.lon,self.lat,aod_fine,im,jm,MISSING) ) + f.write('cloud', nymd, nhms, + binobs2d(self.lon,self.lat,self.cloud,im,jm,MISSING) ) + +# try: +# f.close() +# except: +# pass + + if Verb >=1: + print "[w] Wrote file "+filename + +#--- + def addVar(self,ga,expr='mag(u10m,v10m)',vname='wind',clmYear=None,tight=True): + """ + Given a grads object *ga* having the correct MERRA file as default, + interpolates *var* to obs location and saves it as an attribute + named *vname*. + + If *tight* is True, domain will be restricted conserve memory. This feature + has proven somewhat unstable for reasons yet TBD. + """ + + U = MISSING * np.ones(self.nobs) + if vname == None: + vname = expr + + # nearest time + # ------------ + t = _gatime(self.nymd,self.nhms) + if clmYear != None: + t = t[:-4] + str(clmYear) # replace year + ga('set time '+t,Quiet=True) + + # To conserve memory, restrict domain with 1 gridpoint halo + # --------------------------------------------------------- + if tight: + fh = ga.query("file") + x1, x2 = self.lon.min(),self.lon.max() + y1, y2 = self.lat.min(),self.lat.max() + ga('set lon %f %f'%(x1,x2),Quiet=True) + ga('set lat %f %f'%(y1,y2),Quiet=True) + qh = ga.query("dims") + x1, x2 = (qh.xi[0]-1,qh.xi[1]+1) + y1, y2 = (max(1,qh.yi[0]-1),min(fh.ny,qh.yi[1]+1)) # in [1,ny] + ga('set x %d %d'%(x1,x2),Quiet=True) # make sure x range is int + ga('set y %d %d'%(y1,y2),Quiet=True) # make sure y range is int + expr_ = ga.exp(expr) + else: + expr_ = ga.expr(expr) + u, levs = ga.interp(expr_, self.lon, self.lat ) + U = u.data + if len(np.shape(U)) == 0: + U = U * np.ones(1) # so that we can slice it later + + self.__dict__[vname] = U + +#--- + def getCoxMunk(self,filename='/nobackup/NNR/Misc/coxmunk_lut.npz',channel=550): + """ + Returns ocean albedo. + """ + + # Get precomputed albedo LUT + # -------------------------- + lut = NPZ(filename) + + # Trimmed wind speed + # ------------------ + w10m = self.wind.copy() + w10m[w10m<0] = 0 + w10m[w10m>50.] = 50. + + j = list(lut.channels).index(channel) + + # Interpolate albedo + # ------------------ + albedo = np.zeros(len(w10m)) + albedo[:] = np.interp(w10m,lut.speed,lut.albedo[:,j]) + + self.albedo = albedo + +#............................................................................ + +def granules ( path, algo, inst, syn_time, coll='011', nsyn=8, verbose=False ): + """ + Returns a list of Vx04 granules for a given algorithm at given synoptic time. + On input, + + path --- mounting point for the MxD04 Level 2 files + algo --- either DT or DB + inst --- SNPP + syn_time --- synoptic time (timedate format) + + coll --- collection: 011 (optional) + nsyn --- number of synoptic times per day (optional) + + """ + + # Get product name + # ----------------- + prod = 'AER{}_{}'.format(algo,inst) + + # Determine synoptic time range + # ----------------------------- + dt = timedelta(seconds = 12. * 60. * 60. / nsyn) + t1, t2 = (syn_time-dt,syn_time+dt) + + # Find VIIRS granules in synoptic time range + # ------------------------------------------ + dt = timedelta(minutes=6) + t = datetime(t1.year,t1.month,t1.day,t1.hour,0,0) + Granules = [] + while t < t2: + if t >= t1: + doy = t.timetuple()[7] + basen = "%s/%s/Level2/%s/%04d/%03d/AER%s_L2_VIIRS_%s.A%04d%03d.%02d%02d.%s.*.nc"\ + %(path,coll,prod,t.year,doy,algo,inst,t.year,doy,t.hour,t.minute,coll) + try: + filen = glob(basen)[0] + Granules += [filen,] + if verbose: + print " [x] Found "+filen + except: + pass + t += dt + + if len(Granules) == 0: + print "WARNING: no %s collection %s granules found for"%(prod,coll), syn_time + + return Granules + +#-- + +def print_stats(name,x=None): + "Prints simple stats" + from pylab import prctile + if type(name) is not str: + x = name + name = 'mean,stdv,rms,min,25%,median,75%,max: ' + if name == '__header__': + print '' + n = (80 - len(x))/2 + print n * ' ' + x + print n * ' ' + len(x) * '-' + print '' + print ' Name mean stdv rms min 25% median 75% max' + print ' --------- ------- ------- ------- ------- ------- ------- ------- -------' + elif name == '__sep__': + print ' --------- ------- ------- ------- ------- ------- ------- ------- -------' + elif name == '__footer__': + print ' --------- ------- ------- ------- ------- ------- ------- ------- -------' + print '' + else: + ave = x.mean() + std = x.std() + rms = np.sqrt(ave*ave+std*std) + prc = prctile(x) + print '%10s %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f '%\ + (name,ave,std,rms,prc[0],prc[1],prc[2],prc[3],prc[4]) + +#-- + +def _gatime(nymd,nhms): + Months = ('jan','feb','mar','apr','may','jun','jul','aug','sep','oct','nov','dec') + cymd = "%08d"%int(nymd) + chms = "%06d"%int(nhms) + t = chms[0:2]+":"+chms[2:4]+"Z"+\ + cymd[6:8]+Months[int(cymd[4:6])-1]+cymd[0:4] + return t + +#............................................................................ + +if __name__ == "__main__": + + syn_time = datetime(2013,10,26,10,0,0) + Files = granules('/nobackup/VIIRS/','DB','SNPP',syn_time,coll='011') + + db_ocean = Vx04_L2(Files,'DB','OCEAN',syn_time,Verb=1,only_good=True) + From 700e68737c062cd34de66ee68e234bf204ad5ded Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 29 Sep 2022 10:00:19 -0400 Subject: [PATCH 02/67] PC added a function to binObs to return gridded obs counts --- src/Shared/GMAO_Shared/GMAO_pyobs/binObs_py.F | 63 +++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/binObs_py.F b/src/Shared/GMAO_Shared/GMAO_pyobs/binObs_py.F index 0c04a542..4c573398 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/binObs_py.F +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/binObs_py.F @@ -201,6 +201,69 @@ subroutine binObs3D(lon,lat,obs,nobs,gObs,im,jm,km,missing) end do + end +!----------- + + subroutine binObsCnt3D(lon,lat,obs,nobs,nSample,im,jm,km,missing) +!vb use MAPL_BaseMod, only: MAPL_GetHorzijIndex + implicit NONE + integer, intent(in) :: im, jm, km,nobs + real, intent(in) :: lon(nobs) + real, intent(in) :: lat(nobs) + real, intent(in) :: obs(nobs,km) + real, intent(in) :: missing + real, intent(out) :: nSample(im,jm,km) +! +! Count the Bined obs. It assumes a global 3D GEOS-5 A-Grid: +! P. Castellanos Sep 2022 +! Longitudes in [-180,10] +! Latitudes in [-90,90] +! +! + + integer :: i, j, k, n + real :: dLon, dLat, xLon + integer :: ics(nobs),jcs(nobs) + real, parameter :: radToDeg = 57.2957795 + + dLon = 360. / im + dLat = 180. / ( jm - 1.) + + nSample = 0.0 + do k = 1, km + + if (jm == 6*im)then + print*, "warning, function under dev" +!vb call MAPL_GetHorzijIndex(im,jm,nobs,ics,jcs, +!vb & lon=lon/radToDeg,lat=lat/radToDeg) + +!vb if (abs(obs(n,k)-missing) > 0.01*abs(missing)) then +!vb nSample(ics(n),jcs(n),k) = nSample(ics(n),jcs(n),k) + 1.0 + else + + do n = 1,nobs + + if ( (abs(lon(n))>180.) .OR. (abs(lat(n))>90.) ) cycle + + xLon = lon(n) + if ( xLon >= 180. ) xLon = xLon - 360. + + i = 1 + nint((xlon + 180. ) / dLon) + j = 1 + nint((lat(n) + 90. ) / dLat) + + if ( i>im ) i = i - im + if ( i<1 ) i = i + im + + if (abs(obs(n,k)-missing) > 0.01*abs(missing)) then + nSample(i,j,k) = nSample(i,j,k) + 1.0 + end if + + end do + + end if + + end do + end !------------- subroutine binRms2D(lon,lat,obs,nobs,gObs,im,jm,missing) From d68e7dcaab7f117795ffa747493fea5736c5d1fe Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Tue, 4 Oct 2022 12:45:48 -0400 Subject: [PATCH 03/67] PC - new vxd04_nnr.py script --- src/Applications/GAAS_App/CMakeLists.txt | 1 + src/Applications/GAAS_App/vx04_nnr.py | 518 +++++++++++++++++++++++ 2 files changed, 519 insertions(+) create mode 100644 src/Applications/GAAS_App/vx04_nnr.py diff --git a/src/Applications/GAAS_App/CMakeLists.txt b/src/Applications/GAAS_App/CMakeLists.txt index c4787a11..3932fe11 100644 --- a/src/Applications/GAAS_App/CMakeLists.txt +++ b/src/Applications/GAAS_App/CMakeLists.txt @@ -61,6 +61,7 @@ set (PYFILES avhrr_nnr.py geo04_nnr.py mxd04_nnr.py + vx04_nnr.py ) set (RCFILES diff --git a/src/Applications/GAAS_App/vx04_nnr.py b/src/Applications/GAAS_App/vx04_nnr.py new file mode 100644 index 00000000..e980916a --- /dev/null +++ b/src/Applications/GAAS_App/vx04_nnr.py @@ -0,0 +1,518 @@ +""" +This module implements the VIIRS NNR AOD retrievals. + +This version works from VIIRS DT & DT Level 2 files. + +""" +import os, sys +import warnings +from pyobs.vx04 import Vx04_L2, MISSING, granules, BEST +from ffnet import loadnet +from numpy import c_ as cat +from numpy import copy, ones, sin, cos, exp, arccos, pi, any, log +import numpy as np +from pyobs.bits import BITS + +# SDS to be read in +# ------------ +SDS = dict ( + DT_META = ('longitude', 'latitude', + 'sensor_zenith_angle', 'sensor_azimuth_angle', + 'solar_zenith_angle', 'solar_azimuth_angle', + 'Scattering_Angle', 'Glint_Angle'), + DT_LAND = ( 'Mean_Reflectance_Land', + 'Corrected_Optical_Depth_Land', + 'Surface_Reflectance_Land', + 'Aerosol_Cloud_Fraction_Land', + 'Land_Ocean_Quality_Flag'), + DT_OCEAN = ( 'Effective_Optical_Depth_Average_Ocean', + 'Aerosol_Cloud_Fraction_Ocean', + 'Mean_Reflectance_Ocean', + 'Land_Ocean_Quality_Flag' ), + DB_META = ('Longitude', 'Latitude', 'Scan_Start_Time', + 'Viewing_Zenith_Angle', 'Relative_Azimuth_Angle', + 'Solar_Zenith_Angle', + 'Scattering_Angle' ), + DB_LAND = ('Aerosol_Optical_Thickness_550_Land', + 'Spectral_Aerosol_Optical_Thickness_Land', + 'Spectral_Surface_Reflectance', + 'Spectral_TOA_Reflectance_Land', + 'Aerosol_Optical_Thickness_QA_Flag_Land', + 'Algorithm_Flag_Land', + 'Number_Of_Pixels_Used_Land', + 'Number_Valid_Pixels'), + DB_OCEAN = ('Aerosol_Optical_Thickness_550_Ocean', + 'Spectral_Aerosol_Optical_Thickness_Ocean', + 'Spectral_TOA_Reflectance_Ocean', + 'Aerosol_Optical_Thickness_QA_Flag_Ocean', + 'Algorithm_Flag_Ocean', + 'Number_Of_Pixels_Used_Ocean', + 'Number_Valid_Pixels') + ) + + +# Channels for TOA reflectance +# ----------------------------- +CHANNELS = dict ( + DT_LAND = ( 480, 670, 2250), + DT_OCEAN = ( 480, 550, 670, 860, 1240, 1600, 2250 ), + DB_LAND = ( 412, 488, 550, 670, 865, 1240, 1640, 2250 ), + DB_OCEAN = ( 412, 488, 550, 670, 865, 1240, 1640, 2250 ), + ) + +SCHANNELS = dict ( + DT_LAND = ( 480, 670, 2250 ), + DB_LAND = ( 412, 488, 670 ), + ) + + +# Translate Inputs between NNR and MODIS classes +# ----------------------------------------------- +TranslateInput = dict ( mRef412 = ('reflectance',412), + mRef480 = ('reflectance',480), + mRef550 = ('reflectance',550), + mRef670 = ('reflectance',670), + mRef860 = ('reflectance',860), + mRef1240 = ('reflectance',1240), + mRef1600 = ('reflectance',1600), + mRef2250 = ('reflectance',2250), + mSre412 = ('sfc_reflectance',412), + mSre480 = ('sfc_reflectance',480), + mSre670 = ('sfc_reflectance',670), + mSre2250 = ('sfc_reflectance',2250), + ) + +for var in ( 'ScatteringAngle','GlintAngle', + 'SolarAzimuth', 'SolarZenith', + 'SensorAzimuth','SensorZenith', + 'cloud','qa_flag' ): + TranslateInput[var] = (var,) + +# Translate Targets between ANET and MODIS classes +# ------------------------------------------------ +TranslateTarget = dict ( aTau440 = ( 'aod_', 440 ), + aTau470 = ( 'aod_', 470 ), + aTau500 = ( 'aod_', 500 ), + aTau550 = ( 'aod_', 550 ), + aTau660 = ( 'aod_', 660 ), + aTau870 = ( 'aod_', 870 ), + ) + +class Vx04_NNR(Vx04_L2): + """ + This class extends VIIRS by adding methods for producing + NNR AOD retrievals based on the Neural Net model developed + with class *abc_c6*. + """ + + def __init__(self,l2_path,sat,algo,syn_time,aer_x, + cloud_thresh=0.70, + glint_thresh=40.0, + scat_thresh=170.0, + cloudFree=None, + aodmax=1.0, + coll='011',verbose=0): + """ + Contructs a VX04 object from VIIRS Aerosol Level 2 + granules. On input, + + l2_path --- top directory for the VIIRS Level 2 files; + it must have subdirectories AERDB_inst and AERDT_inst. + sat --- either *SNPP* (Suomi-NPP) or *NOAA20* (NOAA-20) + algo --- aerosol algorithm: DT_LAND, DT_OCEAN, DB_LAND, or DB_OCEAN + syn_time --- synoptic time + + 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 + + The following attributes are also defined: + fractions dust, sea salt, BC+OC, sulfate + + It also performs Q/C, setting attribute iGood. On, + input, *cloud_thresh* is the cloud fraction limit. + + (For right now this is not implemented) + When DEEP BLUE algorithm is requested, filters for + retrievals where DARK TARGET obs are unavailable. + """ + + self.verbose = verbose + self.algo = algo + self.cloudFree = cloudFree + self.aodmax = aodmax + + # Initialize superclass + # --------------------- + Files = granules(l2_path,algo,sat,syn_time,coll=coll) + Vx04_L2.__init__(self,Files,algo,syn_time, + only_good=True, + SDS=SDS, + Verb=verbose) + + if self.nobs < 1: + return # no obs, nothing to do + + # Channels + # ----------------------------- + self.rChannels = CHANNELS[algo] + if algo in SCHANNELS: + self.sChannels = SCHANNELS[algo] + + + # DB Algorithm only used when Dark Target data is unavailable + # (Not currently implemented) + # -------------------------------------------------------------- + +# if algo == "DEEP": +# # Get DARK TARGET qa_flag +# self.qa_flag_lnd = BITS(self.Quality_Assurance_Land[:,0])[1:4] +# lndGood = self.qa_flag_lnd == BEST +# lndGood = lndGood & (self.cloud_lnd < cloud_thresh) +# rChannels = CHANNELS["LAND"] +# sChannels = SCHANNELS["LAND"] +# for i,c in enumerate(rChannels): +# lndGood = lndGood & (self.reflectance_lnd[:,i]>0) + +# for i,c in enumerate(sChannels): +# lndGood = lndGood & (self.sfc_reflectance_lnd[:,i]>0) + +# self.iGood = (self.qa_flag == BEST) & ~lndGood + +# # Keep only "good" observations +# # ----------------------------- +# m = self.iGood +# for sds in self.SDS: +# rank = len(self.__dict__[sds].shape) +# if rank == 1: +# self.__dict__[sds] = self.__dict__[sds][m] +# elif rank == 2: +# self.__dict__[sds] = self.__dict__[sds][m,:] +# else: +# raise IndexError, 'invalid rank=%d'%rank + +# # Reset aliases +# for sds in self.SDS: +# if sds in self.ALIAS: +# self.__dict__[self.ALIAS[sds]] = self.__dict__[sds] + + +# self.qa_flag = self.qa_flag[m] +# self.aod = self.aod[m,:] +# self.time = self.time[m] +# self.Time = self.Time[m] +# self.iGood = self.iGood[m] +# self.nobs = self.Longitude.shape[0] + +# if self.nobs < 1: +# return # no obs, nothing to do + + + # Q/C + # --- + self.iGood = self.cloud<cloud_thresh +# if algo == "LAND": +# self.iGood = self.iGood & (self.cloud_deep<cloud_thresh) +# elif algo == "DEEP": +# self.iGood = self.iGood & (self.cloud_lnd<cloud_thresh) + + for i,c in enumerate(self.rChannels): + self.iGood = self.iGood & (self.reflectance[:,i]>0) + + if algo in SCHANNELS: + for i,c in enumerate(self.sChannels): + self.iGood = self.iGood & (self.sfc_reflectance[:,i]>0) + + if "OCEAN" in algo: + self.iGood = self.iGood & (self.GlintAngle > glint_thresh) + + if "LAND" in algo: + self.iGood = self.iGood & (self.ScatteringAngle < scat_thresh) + + if any(self.iGood) == False: + print "WARNING: Strange, no good obs left to work with" + return + + # Create attribute for holding NNR predicted AOD + # ---------------------------------------------- + self.aod_ = MISSING * ones((self.nobs,len(self.channels))) + + # Make sure same good AOD is kept for gridding + # -------------------------------------------- + if len(self.aod.shape) == 1: + self.aod.shape = self.aod.shape + (1,) + self.aod[self.iGood==False,:] = MISSING + + + # Angle transforms: for NN calculations we work with cosine of angles + # ------------------------------------------------------------------- + self.ScatteringAngle = cos(self.ScatteringAngle*pi/180.0) + self.SensorAzimuth = cos(self.SensorAzimuth*pi/180.0) + self.SensorZenith = cos(self.SensorZenith*pi/180.0) + self.SolarAzimuth = cos(self.SolarAzimuth*pi/180.0) + self.SolarZenith = cos(self.SolarZenith*pi/180.0) + self.GlintAngle = cos(self.GlintAngle*pi/180.0) + + # Get fractional composition + # ------------------------------ + self.speciate(aer_x,Verbose=verbose) + + + def speciate(self,aer_x,Verbose=False): + """ + Use GAAS to derive fractional composition. + """ + + self.sampleFile(aer_x,onlyVars=('TOTEXTTAU', + 'DUEXTTAU', + 'SSEXTTAU', + 'BCEXTTAU', + 'OCEXTTAU', + 'SUEXTTAU', + ),Verbose=Verbose) + + s = self.sample + I = (s.TOTEXTTAU<=0) + s.TOTEXTTAU[I] = 1.E30 + self.fdu = s.DUEXTTAU / s.TOTEXTTAU + self.fss = s.SSEXTTAU / s.TOTEXTTAU + self.fbc = s.BCEXTTAU / s.TOTEXTTAU + self.foc = s.OCEXTTAU / s.TOTEXTTAU + self.fcc = self.fbc + self.foc + self.fsu = s.SUEXTTAU / s.TOTEXTTAU + + # Special handle nitrate (treat it as it were sulfate) + # ---------------------------------------------------- + try: + self.sampleFile(aer_x,onlyVars=('NIEXTTAU',),Verbose=Verbose) + self.fsu += self.sample.NIEXTTAU / s.TOTEXTTAU + except: + pass # ignore it for systems without nitrates + + del self.sample + +#--- + def sampleFile(self, inFile, npzFile=None, onlyVars=None, Verbose=False): + """ + Interpolates all variables of inFile and optionally + save them to file *npzFile* + """ + from gfio import GFIO, GFIOctl, GFIOHandle + + # Instantiate grads and open file + # ------------------------------- + name, ext = os.path.splitext(inFile) + if ext in ( '.nc4', '.nc', '.hdf'): + fh = GFIO(inFile) # open single file + if fh.lm == 1: + timeInterp = False # no time interpolation in this case + else: + raise ValueError, "cannot handle files with more tha 1 time, use ctl instead" + else: + fh = GFIOctl(inFile) # open timeseries + timeInterp = True # perform time interpolation + tymes = np.array([self.syn_time]*self.nobs) + + self.sample = GFIOHandle(inFile) + if onlyVars is None: + onlyVars = fh.vname + + lons = self.lon + lats = self.lat + + + + # Loop over variables on file + # --------------------------- + for v in onlyVars: + if Verbose: + print "<> Sampling ", v + if timeInterp: + var = fh.sample(v,lons,lats,tymes,Verbose=Verbose) + else: + var = fh.interp(v,lons,lats) + if (var.size == 1) & (len(var.shape) == 0): + var.shape = (1,) #protect against when only one value is returned and shape=() + if len(var.shape) == 1: + self.sample.__dict__[v] = var + elif len(var.shape) == 2: + var = var.T # shape should be (nobs,nz) + self.sample.__dict__[v] = var + else: + raise IndexError, 'variable <%s> has rank = %d'%(v,len(var.shape)) + + if npzFile is not None: + savez(npzFile,**self.sample.__dict__) + + + def _loadNet(self,nnFile): + """ + Loads the Neural Net weights created with class ABC. + """ + self.net = loadnet(nnFile) + + def _getInputs(self): + """ + Get Inputs for Neural Net. + """ + + # Loop over inputs + # ---------------- + first = True + for inputName in self.net.InputNames: + try: + iName = TranslateInput[inputName] + except: + iName = inputName + + if self.verbose>0: + print 'Getting NN input ',iName + + # Retrieve input + # -------------- + if type(iName) is str: + input = self.__dict__[iName][:] + + elif len(iName) == 2: + name, ch = iName + if 'mSre' in inputName: # LAND, surface reflectivity + k = list(self.sChannels).index(ch) # index of channel + elif 'mRef' in inputName: # TOA reflectances + k = list(self.rChannels).index(ch) # index of channel + + input = self.__dict__[name][:,k] + + elif len(iName) == 1: + name = iName[0] + input = self.__dict__[name][:] + + else: + raise ValueError, "strange, len(iName)=%d"%len(iName) + + # Concatenate Inputs + # ------------------ + if first: + inputs = input + first = False + else: + inputs = cat[inputs,input] + + # Keep only good observations + # --------------------------- + return inputs[self.iGood,:] + +#-- + def apply(self,nnFile): + """ + Apply bias correction to AOD. + """ + + # Stop here is no good obs available + # ---------------------------------- + if self.nobs == 0: + return # no data to work with + if any(self.iGood) == False: + return # no good data to work with + + # Load the Neural Net + # ------------------- + self._loadNet(nnFile) + + # Evaluate NN on inputs + # --------------------- + targets = self.net(self._getInputs()) + + + # Targets do not have to be in standard retrieval + # ---------------------------------------------- + for i,targetName in enumerate(self.net.TargetNames): + name, ch = TranslateTarget[targetName] + try: + k = list(self.channels).index(ch) # index of channel + except: + # add new target channel to end + self.channels += (ch,) + self.aod = np.append(self.aod,MISSING*ones((self.nobs,1)),axis=1) + self.aod_ = np.append(self.aod_,MISSING*ones((self.nobs,1)),axis=1) + + # Replace targets with unbiased values + # ------------------------------------ + self.channels_ = [] # channels being revised + for i,targetName in enumerate(self.net.TargetNames): + name, ch = TranslateTarget[targetName] + if self.verbose>0: + if self.net.laod: + print "NN Retrieving log(AOD+0.01) at %dnm "%ch + else: + print "NN Retrieving AOD at %dnm "%ch + k = list(self.channels).index(ch) # index of channel + self.channels_ = self.channels_ + [ch,] + if self.net.laod: + result = exp(targets[:,i]) - 0.01 # inverse + else: + result = targets[:,i] + + self.__dict__[name][self.iGood,k] = result + + + # Do extra cloud filtering if required + if self.cloudFree is not None: + cloudy = (self.cloud>=self.cloudFree) +# 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) + + 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] ) + + 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 + + self.iGood[self.iGood][contaminated] = False + + + + + +#--- + + __call__= apply + + +#--- + +if __name__ == "__main__": + + from datetime import datetime + + l2_path = '/nobackup/VIIRS/' + algo = 'DB_LAND' + sat = 'SNPP' + coll = '011' + aer_x = '/nobackup/NNR/Misc/tavg1_2d_aer_Nx' + + syn_time = datetime(2016,12,19,15,0,0) + + if algo == 'OCEAN': + nn_file = '/nobackup/NNR/Net/nnr_003.mydo_Tau.net' + elif algo == 'LAND': + nn_file = '/nobackup/NNR/Net/nnr_003.mydl_Tau.net' + elif algo == 'DEEP': + nn_file = '/nobackup/NNR/Net/nnr_003.mydd_Tau.net' + + m = Vx04_NNR(l2_path,algo.upper(),sat,syn_time,aer_x, + coll=coll, + cloud_thresh=0.7, + verbose=True) + + m.apply(nn_file) + aod = m.aod_ From fe5c8c56abc412af2a90e9122ce8ef18fbcd4d9c Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Tue, 4 Oct 2022 12:53:43 -0400 Subject: [PATCH 04/67] PC - fix surface reflectance, and add counts to binobs --- .../GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 49 +++++++++++-------- 1 file changed, 28 insertions(+), 21 deletions(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index d9820f3c..424001e5 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -16,7 +16,7 @@ from glob import glob from dateutil.parser import isoparse from pyobs.npz import NPZ -from binObs_ import binobs2d, binobs3d +from binObs_ import binobs2d, binobs3d, binobscnt3d try: from pyods import ODS # must be inported before pyhdf @@ -90,7 +90,8 @@ DT_OCEAN = ( 480., 550., 670., 860., 1240., 1600., 2250. ), DB_LAND = ( 412, 488, 550, 670 ), # file has 550 separate DB_OCEAN = (488., 550., 670., 865., 1240., 1640., 2250.), - SREF = ( 480., 670., 2250. ) + DT_SREF = ( 480., 670., 2250. ), + DB_SREF = (412., 488., 670. ), ) @@ -170,7 +171,7 @@ class Vx04_L2(object): Therefore, we will refer to these products at VSNPP04 and VN2004, respectively """ - def __init__ (self,Path,Algo,Surface,syn_time=None,nsyn=8,Verb=0, + def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, only_good=True,SDS=SDS,alias=None): """ Reads individual granules or a full day of Level 2 Vx04 files @@ -182,8 +183,7 @@ def __init__ (self,Path,Algo,Surface,syn_time=None,nsyn=8,Verb=0, of files and directories. Directories are transversed recursively. If a non Vx04 Level 2 file is encountered, it is simply ignored. - Algo -- Algorithm: DT or DB - Surface -- OCEAN or LAND + algo -- Algorithm: DT_LAND, DT_OCEAN, DB_LAND or DB_OCEAN Optional parameters: syn_type --- synoptic time @@ -199,7 +199,6 @@ def __init__ (self,Path,Algo,Surface,syn_time=None,nsyn=8,Verb=0, """ - algo = '{}_{}'.format(Algo,Surface) if algo not in ('DT_LAND', 'DT_OCEAN', 'DB_LAND', 'DB_OCEAN'): raise ValueError, "invalid algorithm "+algo+" --- must be DT_LAND, DT_OCEAN, DB_LAND, DB_OCEAN" @@ -209,6 +208,7 @@ def __init__ (self,Path,Algo,Surface,syn_time=None,nsyn=8,Verb=0, self.sat = None # Satellite name self.col = None # collection, e.g., 011 self.algo = algo + Algo, Surface = algo.split('_') self.SDS = SDS['{}_META'.format(Algo)] + SDS[algo] self.SDS_META = SDS['{}_META'.format(Algo)] @@ -295,14 +295,16 @@ def __init__ (self,Path,Algo,Surface,syn_time=None,nsyn=8,Verb=0, sza = np.radians(self.SolarZenith) vza = np.radians(self.SensorZenith) raa = np.radians(self.RelativeAzimuth) - cglint = np.cos(sza)*np.cos(vza) + np.sin(sza)*np.sin(vza)*np.cos(np.pi-raa) - self.GlintAngle = np.arccos(cglint) + cglint = np.cos(sza)*np.cos(vza) + np.sin(sza)*np.sin(vza)*np.cos(raa) + self.GlintAngle = np.degrees(np.arccos(cglint)) # Create corresponding python time # -------------------------------- if 'DB' in self.algo: self.Time = np.array([DATE_START+timedelta(seconds=s) for s in self.Scan_Start_Time]) + else: + self.Time = np.array(self.Time) # masked datetime arrays aren't friendly # ODS friendly attributes # ----------------------- @@ -310,10 +312,8 @@ def __init__ (self,Path,Algo,Surface,syn_time=None,nsyn=8,Verb=0, self.kx = KX[self.sat+'_'+self.algo] self.ident = IDENT[self.sat+'_'+self.algo] self.channels = CHANNELS[self.algo] - if Algo == 'DT': - self.sChannels = CHANNELS["SREF"] # LAND surface reflectivity (not the same as algo) - else: - self.sChannels = CHANNELS[self.algo] # deep blue surface bands are the same as algorithm + if Surface == 'LAND': + self.sChannels = CHANNELS["{}_SREF".format(Algo)] # LAND surface reflectivity (not the same as algo) if 'DB' in self.algo: self.rChannels = self.Reflectance_Bands @@ -706,11 +706,13 @@ def writeg(self,filename=None,dir='.',expid=None,refine=8,res=None, vtitle = [ 'Aerosol Optical Depth', 'Aerosol Optical Depth (Revised)', 'Aerosol Optical Depth (Fine Mode)', + 'Aerosol Optical Depth Obs Count', + 'Aerosol Optical Depth (Revised) Obs Count', 'Cloud Fraction' ] - vname = ['tau', 'tau_', 'tau_fine', 'cloud' ] - vunits = [ '1', '1', '1', '1', ] - kmvar = [ nch, nch, nch, 0 ] + vname = ['tau', 'tau_', 'tau_fine', 'count_tau', 'count_tau_', 'cloud' ] + vunits = [ '1', '1', '1', '1', '1', '1', ] + kmvar = [ nch, nch, nch, nch, nch, 0 ] title = 'Gridded MODIS Aerosol Retrievals' source = 'NASA/GSFC/GMAO GEOS-5 Aerosol Group' @@ -757,6 +759,10 @@ def writeg(self,filename=None,dir='.',expid=None,refine=8,res=None, binobs3d(self.lon,self.lat,aod_,im,jm,MISSING) ) f.write('tau_fine', nymd, nhms, binobs3d(self.lon,self.lat,aod_fine,im,jm,MISSING) ) + f.write('count_tau', nymd, nhms, + binobscnt3d(self.lon,self.lat,aod,im,jm,MISSING) ) + f.write('count_tau_', nymd, nhms, + binsobscnt3d(self.lon,self.lat,aod_,im,jm,MISSING) ) f.write('cloud', nymd, nhms, binobs2d(self.lon,self.lat,self.cloud,im,jm,MISSING) ) @@ -840,14 +846,14 @@ def getCoxMunk(self,filename='/nobackup/NNR/Misc/coxmunk_lut.npz',channel=550): #............................................................................ -def granules ( path, algo, inst, syn_time, coll='011', nsyn=8, verbose=False ): +def granules ( path, algo, sat, syn_time, coll='011', nsyn=8, verbose=False ): """ Returns a list of Vx04 granules for a given algorithm at given synoptic time. On input, path --- mounting point for the MxD04 Level 2 files - algo --- either DT or DB - inst --- SNPP + algo --- either DT_LAND, DT_OCEAN, DB_LAND or DB_OCEAN + sat --- SNPP syn_time --- synoptic time (timedate format) coll --- collection: 011 (optional) @@ -857,7 +863,8 @@ def granules ( path, algo, inst, syn_time, coll='011', nsyn=8, verbose=False ): # Get product name # ----------------- - prod = 'AER{}_{}'.format(algo,inst) + Algo = algo.split('_')[0] + prod = 'AER{}_{}'.format(Algo,sat) # Determine synoptic time range # ----------------------------- @@ -873,7 +880,7 @@ def granules ( path, algo, inst, syn_time, coll='011', nsyn=8, verbose=False ): if t >= t1: doy = t.timetuple()[7] basen = "%s/%s/Level2/%s/%04d/%03d/AER%s_L2_VIIRS_%s.A%04d%03d.%02d%02d.%s.*.nc"\ - %(path,coll,prod,t.year,doy,algo,inst,t.year,doy,t.hour,t.minute,coll) + %(path,coll,prod,t.year,doy,Algo,sat,t.year,doy,t.hour,t.minute,coll) try: filen = glob(basen)[0] Granules += [filen,] @@ -884,7 +891,7 @@ def granules ( path, algo, inst, syn_time, coll='011', nsyn=8, verbose=False ): t += dt if len(Granules) == 0: - print "WARNING: no %s collection %s granules found for"%(prod,coll), syn_time + print "WARNING: no %s collection %s granules found for"%(algo,coll), syn_time return Granules From 4cb763d99e281debf0127e609b2793ac786f81d1 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Wed, 5 Oct 2022 17:47:31 -0400 Subject: [PATCH 05/67] PC - trying to get a DT relative azimuth angle. not provided in files --- src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index 424001e5..866c69e0 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -297,7 +297,16 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, raa = np.radians(self.RelativeAzimuth) cglint = np.cos(sza)*np.cos(vza) + np.sin(sza)*np.sin(vza)*np.cos(raa) self.GlintAngle = np.degrees(np.arccos(cglint)) - + elif 'DT' in self.algo: + # this kind of seems to match the DB RAA + # the DT sensor and solar azimuth angles seem off + # I can't find a DB definition for RAA so this is close enough + raa = self.SolarZenith - self.SensoZenith + ii = raa <0 + raa[ii] = raa[ii] + 180. + ii = raa < 0 + raa[ii] = raa[ii]*-1. + self.RelativeAzimuth = raa # Create corresponding python time # -------------------------------- From 23b00b0eaa1b82f2f501989cd091a620c328f358 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 10 Oct 2022 18:17:18 -0400 Subject: [PATCH 06/67] PC - add new viirs scripts to cmakelists --- src/Components/misc/obs_aod/ABC/CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Components/misc/obs_aod/ABC/CMakeLists.txt b/src/Components/misc/obs_aod/ABC/CMakeLists.txt index cb7fc02a..269916d7 100644 --- a/src/Components/misc/obs_aod/ABC/CMakeLists.txt +++ b/src/Components/misc/obs_aod/ABC/CMakeLists.txt @@ -23,6 +23,7 @@ set (PYFILES abc_avhrr.py abc_c6_aux.py abc_c6.py + abc_viirs.py abc_misr.py abc_modis.py albedo.py @@ -35,6 +36,7 @@ set (PYFILES evaluate_avhrr.py evaluate_cdr.py giant.py + giant_viirs.py man.py mcd43c.py nn.py From dcd6c7f099376b993008dc8f666e95c952b30528 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 10 Oct 2022 18:18:23 -0400 Subject: [PATCH 07/67] PC - fix time format and add derived angles to sds list --- src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index 866c69e0..71649253 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -297,16 +297,18 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, raa = np.radians(self.RelativeAzimuth) cglint = np.cos(sza)*np.cos(vza) + np.sin(sza)*np.sin(vza)*np.cos(raa) self.GlintAngle = np.degrees(np.arccos(cglint)) + self.SDS += ('GlintAngle',) elif 'DT' in self.algo: # this kind of seems to match the DB RAA # the DT sensor and solar azimuth angles seem off # I can't find a DB definition for RAA so this is close enough - raa = self.SolarZenith - self.SensoZenith + raa = self.SolarZenith - self.SensorZenith ii = raa <0 raa[ii] = raa[ii] + 180. ii = raa < 0 raa[ii] = raa[ii]*-1. self.RelativeAzimuth = raa + self.SDS += ('RelativeAzimuth',) # Create corresponding python time # -------------------------------- @@ -325,7 +327,7 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, self.sChannels = CHANNELS["{}_SREF".format(Algo)] # LAND surface reflectivity (not the same as algo) if 'DB' in self.algo: - self.rChannels = self.Reflectance_Bands + self.rChannels = self.Reflectance_Bands # [ 412., 488., 550., 670., 865., 1240., 1640., 2250.] elif self.algo == 'DT_LAND': self.rChannels = np.array([480.,670.,2250.]) elif self.algo == 'DT_OCEAN': @@ -528,7 +530,13 @@ def reduce(self,I): if sds in Alias: self.__dict__[self.ALIAS[sds]] = self.__dict__[sds] # redefine aliases - self.nobs = len(self.lon) + self.nobs = len(self.lon) + # Create corresponding python time + # -------------------------------- + if 'DB' in self.algo: + self.Time = np.array([DATE_START+timedelta(seconds=s) for s in self.Scan_Start_Time]) + else: + self.Time = np.array(self.Time) # masked datetime arrays aren't friendly #--- def write(self,filename=None,dir='.',expid=None,Verb=1): From 9b2e08fe9c6585158fe3f7fbb1767512c0c10e87 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 10 Oct 2022 18:18:55 -0400 Subject: [PATCH 08/67] PC - new abs and giant script for viirs --- src/Components/misc/obs_aod/ABC/abc_viirs.py | 873 ++++++++++++++++++ .../misc/obs_aod/ABC/giant_viirs.py | 608 ++++++++++++ 2 files changed, 1481 insertions(+) create mode 100644 src/Components/misc/obs_aod/ABC/abc_viirs.py create mode 100644 src/Components/misc/obs_aod/ABC/giant_viirs.py diff --git a/src/Components/misc/obs_aod/ABC/abc_viirs.py b/src/Components/misc/obs_aod/ABC/abc_viirs.py new file mode 100644 index 00000000..a0aae9bf --- /dev/null +++ b/src/Components/misc/obs_aod/ABC/abc_viirs.py @@ -0,0 +1,873 @@ +""" + This module implements a Neural Net based MODIS Collection 6 Neural Net Retrieval. + + Arlindo da Silva, June 2015. + +""" + +import os, sys +from matplotlib.pyplot import cm, imshow, plot, figure +from matplotlib.pyplot import xlabel, ylabel, title, grid, savefig, legend +import matplotlib.pyplot as plt +from matplotlib.ticker import MultipleLocator +import matplotlib.patches as mpatches +from numpy import c_ as cat +from numpy import random, sort, pi, load, cos, log, std, exp +from numpy import reshape, arange, ones, zeros, interp, sqrt +from numpy import meshgrid, concatenate, squeeze +import numpy as np +from giant_viirs import MISSING, DT_LAND, DT_OCEAN, DB_LAND, DB_OCEAN +from nn import NN, _plotKDE +import itertools +from sklearn.linear_model import LinearRegression +from multiprocessing import cpu_count +from abc_c6_aux import SummarizeCombinations, get_Iquartiles, get_Ispecies, get_ImRef +from abc_c6_aux import make_plots, make_plots_angstrom, TestStats, SummaryPDFs +from brdf import rtlsReflectance +from mcd43c import BRDF + +# ------ + +VARNAMES = {'cloud': 'MOD04 Cloud Fraction', + 'ScatteringAngle': 'Scattering Angle', + 'GlintAngle': 'Glint Angle', + 'AMF': 'Air Mass Factor', + 'SolarZenith': 'Solar Zenith Angle', + 'CoxMunkLUT': 'Cox-Munk White Sky Albedo', + 'COxMunkBRF': 'Cox-Munk Bidirectional Surface Reflectance', + 'MOD43BClimAlbedo': 'MOD43B Albedo Climatology', + 'fdu': 'MERRA2 Fraction Dust Aerosol', + 'fcc': 'MERRA2 Fraction Carbonaceous Aerosol', + 'fsu': 'MERRA2 Fraction Sulfate Aerosol', + 'year': 'Year'} +#-------------------------------------------------------------------------------------- +class SETUP(object): + def setupNN(self,retrieval,expid, + nHidden=None, + nHLayers=1, + combinations=False, + Input_nnr = ['mRef470','mRef550','mRef660', 'mRef870', + 'mRef1200','mRef1600','mRef2100', + 'ScatteringAngle', 'GlintAngle', + 'AMF', 'SolarZenith', + 'cloud', 'albedo','fdu','fcc','fsu' ], + Input_const = None, + Target = ['aTau550',], + K=None): + + + self.retrieval = retrieval + # Create outdir if it doesn't exist + # --------------------------------- + self.outdir = "./{}/".format(expid) + if not os.path.exists(self.outdir): + os.makedirs(self.outdir) + + self.plotdir = self.outdir + + # save some inputs + # ----------------- + self.expid = expid + self.Target = Target + self.nTarget = len(Target) + self.K = K + self.nHidden = nHidden + + # figure out if you need to calculate angstrom exponent + angstrom = False + for tname in Target: + if not angstrom: + if 'AE' in tname: + angstrom = True + self.angstrom = angstrom + + # if angstrom is being trained + # find the base wavelength + # calculate angstrom with respect to the base wavelength + # ------------------------------------------------------- + if angstrom: + # find base wavelength + for i,tname in enumerate(Target): + if 'Tau' in tname: + base_name = tname + base_wavs = tname.split('Tau')[-1] + base_wav = float(base_wavs) + base_tau = self.__dict__[tname] + base_wav_i = i + + self.AE_base_wav = base_wav + self.AE_base_wav_i = base_wav_i + + # Calculate the angstrom exponent + # with respect to the base wavelength + for tname in Target: + if 'Tau' not in tname: + wavs = tname.split('AE')[-1] + wav = float(wavs) + tau = self.__dict__['aTau'+wavs] + AE = -1.*np.log(tau/base_tau)/np.log(wav/base_wav) + self.__dict__['aAE'+wavs] = AE + + # Balance the dataset before splitting + # No aerosol type should make up more that 35% + # of the total number of obs + # -------------------------------------- + self.iValid = self.balance(int(self.nobs*0.35)) + + # Flatten Input_nnr into one list + # ------------------------------- + input_nnr = flatten_list(Input_nnr) + + # Create list of combinations + # --------------------------- + if combinations: + self.comblist, self.combgroups = get_combinations(Input_nnr,Input_const) + else: + self.comblist = [input_nnr] + + # Initialize arrays to hold stats + # ------------------------------ + self.nnr = STATS(K,self.comblist,self.nTarget) + self.orig = STATS(K,self.comblist,self.nTarget) + + # Initialize K-folding + # -------------------- + if K is None: + self.iTest = ones([self.nobs]).astype(bool) + self.iTrain = self.iValid + else: + self.kfold(K=K) + + # Create list of topologies + # ------------------------- + self.topology = [] + if not combinations: + if self.nHidden is None: + self.nHidden = len(input_nnr) + else: + self.nHidden = nHidden + + self.topology.append((len(input_nnr),) + (self.nHidden,)*nHLayers + (len(Target),)) + + else: + for c,Input in enumerate(self.comblist): + if nHidden is None: + self.nHidden = len(Input) + else: + self.nHidden = nHidden + + self.topology.append((len(Input),) + (self.nHidden,)*nHLayers + (len(Target),)) + + self.combinations = combinations + +#---------------------------------------------------------------------------- +class ABC(object): + + """ + Common Subroutines to all the ABC Classes + """ + + def __init__(self,fname,Albedo,coxmunk_lut=None,NDVI=False): + + # Get Auxiliary Data + self.fnameRoot = fname[:-3] + self.setWind() + self.setAlbedo(Albedo,coxmunk_lut=coxmunk_lut) + self.setSpec() + if NDVI: + self.setNDVI() + + def setWind(self): + # Read in wind + # ------------------------ + self.wind = load(self.fnameRoot + "_MERRA2.npz")['wind'] + self.giantList.append('wind') + self.Wind = '' #need this for backwards compatibility + + def setAlbedo(self,Albedo,coxmunk_lut=None): + # Define wind speed dependent ocean albedo + # ---------------------------------------- + if Albedo is not None: + for albedo in Albedo: + if albedo == 'CoxMunkLUT': + self.getCoxMunk(coxmunk_lut) + self.giantList.append(albedo) + elif albedo == 'MCD43C1': + self.setBRDF() + elif 'CxAlbedo' in albedo : + self.setCoxMunkBRF(albedo) + else: + self.__dict__[albedo] = squeeze(load(self.fnameRoot+'_'+albedo+'.npz')["albedo"]) + self.giantList.append(albedo) + + def setSpec(self): + # Read in Aerosol Fractional Composition + # -------------------------------------- + names = ('fdu','fss','fcc','fsu') + for name in names: + self.__dict__[name] = load(self.fnameRoot + "_MERRA2.npz")[name] + self.giantList.append(name) + + def setCoxMunkBRF(self,albedo): + # Read in Cox Munk Bidirectional surface reflectance + # -------------------------------------------------- + names = ['470','550','660','870','1200','1600','2100'] + for ch in names: + name = 'CxAlbedo' + ch + self.__dict__[name] = squeeze(load(self.fnameRoot+'_CxAlbedo.npz')[name]) + self.giantList.append(name) + + def setBRDF(self): + # Read in MCD43C1 BRDF + # Calculate bidirectional surface reflectance + # --------------------------------------------- + brdf = BRDF(self.nobs) + names = ('BRDFvis','BRDFnir','BRDF470','BRDF550', + 'BRDF650','BRDF850','BRDF1200','BRDF1600', + 'BRDF2100') + for name in brdf.__dict__: + brdf.__dict__[name] = load(self.fnameRoot + "_MCD43C1.npz")[name] + + for name in names: + ch = name[4:] + Kiso = brdf.__dict__['Riso'+ch] + Kvol = brdf.__dict__['Rvol'+ch] + Kgeo = brdf.__dict__['Rgeo'+ch] + self.__dict__[name] = rtlsReflectance(Kiso,Kgeo,Kvol, + self.SolarZenith,self.SensorZenith, + self.SolarAzimuth,self.SensorAzimuth) + self.giantList.append(name) + + + def setNDVI(self): + # Read in NDVI + # ------------- + names = ('NDVI','EVI','NIRref') + for name in names: + self.__dict__[name] = load(self.fnameRoot + "_NDVI.npz")[name] + self.giantList.append(name) + + def outlierRemoval(self,outliers): + + # # Outlier removal based on log-transformed AOD + # # -------------------------------------------- + if outliers > 0.: + d = log(self.mTau550[self.iValid]+0.01) - log(self.aTau550[self.iValid]+0.01) + if self.verbose>0: + print "Outlier removal: %d sig_d = %f nGood=%d "%(-1,std(d),d.size) + for iter in range(3): + iValid = (abs(d)<outliers*std(d)) + self.iValid[self.iValid] = iValid + d = log(self.mTau550[self.iValid]+0.01) - log(self.aTau550[self.iValid]+0.01) + if self.verbose>0: + print "Outlier removal: %d sig_d = %f nGood=%d "%(iter,std(d),d.size) + + def angleTranform(self): + # Angle transforms: for NN work we work with cosine of angles + # ----------------------------------------------------------- + self.ScatteringAngle = cos(self.ScatteringAngle*pi/180.0) + self.SensorZenith = cos(self.SensorZenith*pi/180.0) + self.SolarZenith = cos(self.SolarZenith*pi/180.0) + self.GlintAngle = cos(self.GlintAngle*pi/180.0) + self.AMF = (1/self.SolarZenith) + (1/self.SensorZenith) + self.giantList.append('AMF') + + def setYear(self): + # Year + #------- + self.year = np.array([t.year for t in self.tyme]) + self.giantList.append('year') + + def addFilter(self,aFilter): + if aFilter is not None: + filters = [] + for f in aFilter: + filters.append(self.__dict__[f]>0) + + oiValid = reduce(lambda x,y: x&y,filters) + self.iValid = self.iValid & oiValid + +#---------------------------------------------------------------------------- +class ABC_DT_Ocean (DT_OCEAN,NN,SETUP,ABC): + + def __init__ (self,fname, + coxmunk_lut='/nobackup/NNR/Misc/coxmunk_lut.npz', + outliers=3., + laod=True, + verbose=0, + cloud_thresh=0.70, + glint_thresh=40.0, + Albedo=None, + aFilter=None, + tymemax=None): + """ + Initializes the AOD Bias Correction (ABC) for the VIIRS DT Ocean algorithm. + + On Input, + + fname --- file name for the CSV file with the co-located MODIS/AERONET + data (see class OCEAN) + outliers -- number of standard deviations for outlinear removal. + laod --- if True, targets are log-transformed AOD, log(Tau+0.01) + tymemax --- truncate the data record in the giant file at tymemax. + set to None to read entire data record. + + + Reads in two Albedo variables + albedo - cox munk lut that parameterizes albedo with wind speed + Requires coxmunk_lut npz file. + BRF - cox munk bidirection reflectance computed with VLIDORT + and stored in npz + Requires a NPZ file with the data. + Both require a wind speed npz file. + """ + + self.verbose = verbose + self.laod = laod + + DT_OCEAN.__init__(self,fname,tymemax=tymemax) # initialize superclass + + # Get Auxiliary Data + ABC.__init__(self,fname,Albedo,coxmunk_lut=coxmunk_lut) + + # Q/C + # --- + self.iValid = (self.qa>0) & \ + (self.aTau470 > -0.01) &\ + (self.aTau550 > -0.01) &\ + (self.aTau660 > -0.01) &\ + (self.aTau870 > -0.01) &\ + (self.mTau480 > -0.01) &\ + (self.mTau550 > -0.01) &\ + (self.mTau670 > -0.01) &\ + (self.mTau870 > -0.01) &\ + (self.mRef480 > 0.0) &\ + (self.mRef550 > 0.0) &\ + (self.mRef670 > 0.0) &\ + (self.mRef860 > 0.0) &\ + (self.mRef1240 > 0.0) &\ + (self.mRef1600 > 0.0) &\ + (self.mRef2250 > 0.0) &\ + (self.cloud <cloud_thresh) &\ + (self.cloud >= 0) &\ + (self.GlintAngle != MISSING ) &\ + (self.GlintAngle > glint_thresh) + + # Filter by additional variables + # ------------------------------ + self.addFilter(aFilter) + + # glint_thresh > 40 is a bit redundant b/c MOD04 should already give these a qa==0 or + # does not retrieve. However, there are a few cases (~200) where this does not happen. + # the GlingAngle is very close to 40, greater than 38. Not sure why these get through. + + # Outlier removal based on log-transformed AOD + # -------------------------------------------- + self.outlierRemoval(outliers) + + # Reduce the Dataset + # -------------------- + self.reduce(self.iValid) + self.iValid = ones(self.lon.shape).astype(bool) + + # Angle transforms: for NN work we work with cosine of angles + # ----------------------------------------------------------- + self.angleTranform() + + # Year + #------- + self.setYear() + +#---------------------------------------------------------------------------- +class ABC_DB_Ocean (DB_OCEAN,NN,SETUP,ABC): + + def __init__ (self,fname, + coxmunk_lut='/nobackup/NNR/Misc/coxmunk_lut.npz', + outliers=3., + laod=True, + verbose=0, + cloud_thresh=0.70, + glint_thresh=40.0, + Albedo=None, + aFilter=None, + tymemax=None): + """ + Initializes the AOD Bias Correction (ABC) for the VIIRS DB Ocean algorithm. + + On Input, + + fname --- file name for the CSV file with the co-located MODIS/AERONET + data (see class OCEAN) + outliers -- number of standard deviations for outlinear removal. + laod --- if True, targets are log-transformed AOD, log(Tau+0.01) + tymemax --- truncate the data record in the giant file at tymemax. + set to None to read entire data record. + + + Reads in two Albedo variables + albedo - cox munk lut that parameterizes albedo with wind speed + Requires coxmunk_lut npz file. + BRF - cox munk bidirection reflectance computed with VLIDORT + and stored in npz + Requires a NPZ file with the data. + Both require a wind speed npz file. + """ + + self.verbose = verbose + self.laod = laod + + DB_OCEAN.__init__(self,fname,tymemax=tymemax) # initialize superclass + + # Get Auxiliary Data + ABC.__init__(self,fname,Albedo,coxmunk_lut=coxmunk_lut) + + # Q/C + # --- + self.iValid = (self.qa>0) & \ + (self.aTau470 > -0.01) &\ + (self.aTau550 > -0.01) &\ + (self.aTau660 > -0.01) &\ + (self.aTau870 > -0.01) &\ + (self.mTau480 > -0.01) &\ + (self.mTau550 > -0.01) &\ + (self.mTau670 > -0.01) &\ + (self.mTau865 > -0.01) &\ + (self.mTau1240 > -0.01) &\ + (self.mTau1640 > -0.01) &\ + (self.mTau2250 > -0.01) &\ + (self.mRef412 > 0.0) &\ + (self.mRef488 > 0.0) &\ + (self.mRef550 > 0.0) &\ + (self.mRef670 > 0.0) &\ + (self.mRef865 > 0.0) &\ + (self.mRef1240 > 0.0) &\ + (self.mRef1640 > 0.0) &\ + (self.mRef2250 > 0.0) &\ + (self.cloud <cloud_thresh) &\ + (self.cloud >= 0) &\ + (self.GlintAngle != MISSING ) &\ + (self.GlintAngle > glint_thresh) + + # Filter by additional variables + # ------------------------------ + self.addFilter(aFilter) + + # glint_thresh > 40 is a bit redundant b/c MOD04 should already give these a qa==0 or + # does not retrieve. However, there are a few cases (~200) where this does not happen. + # the GlingAngle is very close to 40, greater than 38. Not sure why these get through. + + # Outlier removal based on log-transformed AOD + # -------------------------------------------- + self.outlierRemoval(outliers) + + # Reduce the Dataset + # -------------------- + self.reduce(self.iValid) + self.iValid = ones(self.lon.shape).astype(bool) + + # Angle transforms: for NN work we work with cosine of angles + # ----------------------------------------------------------- + self.angleTranform() + + # Year + #------- + self.setYear() +#---------------------------------------------------------------------------- + +class ABC_DT_Land (DT_LAND,NN,SETUP,ABC): + + def __init__ (self, fname, + Albedo=None, + alb_max = 0.25, + outliers=3., + laod=True, + verbose=0, + cloud_thresh=0.70, + aFilter=None, + NDVI=False, + tymemax=None): + """ + Initializes the AOD Bias Correction (ABC) for the VIIRS DT Land algorithm. + + On Input, + + fname --- file name for the CSV file with the co-located VIIRS/AERONET + data (see class OCEAN) + + Albedo --- albedo file name identifier; albedo file will be created + from this identifier (See below). + outliers -- number of standard deviations for outlinear removal. + laod --- if True, targets are log-transformed AOD, log(Tau+0.01) + tymemax --- truncate the data record in the giant file at tymemax. + set to None to read entire data record. + """ + + self.verbose = verbose + self.laod = laod + + DT_LAND.__init__(self,fname,tymemax=tymemax) # initialize superclass + + # Get Auxiliary Data + ABC.__init__(self,fname,Albedo,NDVI=NDVI) + + # Q/C: enforce QA=3 and albedo in (0,0.25), scattering angle<170 + # -------------------------------------------------------------- + self.iValid = (self.qa==3) & \ + (self.aTau470 > -0.01) & \ + (self.aTau550 > -0.01) & \ + (self.aTau660 > -0.01) & \ + (self.mTau480 > -0.01) & \ + (self.mTau550 > -0.01) & \ + (self.mTau670 > -0.01) & \ + (self.mTau2250> -0.01) & \ + (self.cloud<cloud_thresh) & \ + (self.cloud >= 0) & \ + (self.ScatteringAngle<170.) & \ + (self.mRef480 > 0) & \ + (self.mRef670 > 0) & \ + (self.mRef2250 > 0) & \ + (self.mSre480 > 0.0) & \ + (self.mSre670 > 0.0) & \ + (self.mSre2250> 0.0) + +# (self.mRef412 > 0) & \ +# (self.mRef440 > 0) & \ +# (self.mRef550 > 0) & \ +# (self.mRef870 > 0) & \ +# (self.mRef1200 > 0) & \ +# (self.mRef1600 > 0) & \ + + # Filter by additional variables + # ------------------------------ + self.addFilter(aFilter) + + + # Outlier removal based on log-transformed AOD + # -------------------------------------------- + self.outlierRemoval(outliers) + + # Reduce the Dataset + # -------------------- + self.reduce(self.iValid) + self.iValid = ones(self.lon.shape).astype(bool) + + # Angle transforms: for NN work we work with cosine of angles + # ----------------------------------------------------------- + self.angleTranform() +#---------------------------------------------------------------------------- + +class ABC_DB_Land (DB_LAND,NN,SETUP,ABC): + + def __init__ (self, fname, + Albedo=None, + outliers=3., + laod=True, + verbose=0, + cloud_thresh=0.70, + aFilter=None, + NDVI=False, + tymemax=None): + """ + Initializes the AOD Bias Correction (ABC) for the MODIS Land algorithm. + + On Input, + + fname --- file name for the CSV file with the co-located MODIS/AERONET + data (see class OCEAN) + + Albedo --- albedo file name identifier; albedo file will be created + from this identifier (See below). + outliers -- number of standard deviations for outlinear removal. + laod --- if True, targets are log-transformed AOD, log(Tau+0.01) + tymemax --- truncate the data record in the giant file at tymemax. + set to None to read entire data record. + """ + + self.verbose = verbose + self.laod = laod + + DB_LAND.__init__(self,fname,tymemax=tymemax) # initialize superclass + + # Get Auxiliary Data + ABC.__init__(self,fname,Albedo,NDVI=NDVI) + + # Q/C: enforce QA=3, scattering angle<170 + # -------------------------------------------------------------- + self.iValid = (self.qa==3) & \ + (self.aTau470 > -0.01) & \ + (self.aTau550 > -0.01) & \ + (self.aTau660 > -0.01) & \ + (self.mTau412 > -0.01) & \ + (self.mTau480 > -0.01) & \ + (self.mTau550 > -0.01) & \ + (self.mTau670 > -0.01) & \ + (self.mTau2250 > -0.01) & \ + (self.cloud<cloud_thresh) & \ + (self.cloud >= 0) & \ + (self.ScatteringAngle<170.) & \ + (self.mRef412 > 0) & \ + (self.mRef488 > 0) & \ + (self.mRef670 > 0) & \ + (self.mRef865 > 0) & \ + (self.mRef1240 > 0) & \ + (self.mRef1640 > 0) & \ + (self.mRef2250 > 0) & \ + (self.mSre412 > 0.0) & \ + (self.mSre488 > 0.0) & \ + (self.mSre670 > 0.0) + + + # Filter by additional variables + # ------------------------------ + self.addFilter(aFilter) + + # Outlier removal based on log-transformed AOD + # -------------------------------------------- + self.outlierRemoval(outliers) + + # Reduce the Dataset + # -------------------- + self.reduce(self.iValid) + self.iValid = ones(self.lon.shape).astype(bool) + + + # Angle transforms: for NN work we work with cosine of angles + # ----------------------------------------------------------- + self.angleTranform() + +#---------------------------------------------------------------------------- +class STATS(object): + + def __init__ (self,K,comblist,nTarget): + c = max([len(comblist),1]) + if K is None: + k = 1 + else: + k = K + + self.slope = np.ones([k,c,nTarget])*np.nan + self.intercept = np.ones([k,c,nTarget])*np.nan + self.R = np.ones([k,c,nTarget])*np.nan + self.rmse = np.ones([k,c,nTarget])*np.nan + self.mae = np.ones([k,c,nTarget])*np.nan + self.me = np.ones([k,c,nTarget])*np.nan + +#--------------------------------------------------------------------- +def _train(mxd,expid,c): + + ident = mxd.ident + outdir = mxd.outdir + + nHidden = mxd.nHidden + topology = mxd.topology[c] + + Input = mxd.comblist[c] + Target = mxd.Target + + print "-"*80 + print "--> nHidden = ", nHidden + print "--> Inputs = ", Input + + n = cpu_count() + kwargs = {'nproc' : n} + if mxd.K is None: + mxd.train(Input=Input,Target=Target,nHidden=nHidden,topology=topology,**kwargs) + mxd.savenet(outdir+"/"+expid+'_Tau.net') + else: + k = 1 + for iTrain, iTest in mxd.kf: + I = arange(mxd.nobs) + iValid = I[mxd.iValid] + mxd.iTrain = iValid[iTrain] + mxd.iTest = iValid[iTest] + + mxd.train(Input=Input,Target=Target,nHidden=nHidden,topology=topology,**kwargs) + mxd.savenet(outdir+"/"+expid+'.k={}_Tau.net'.format(str(k))) + k = k + 1 +#-------------------------------------------------------------------------------------- + +def _test(mxd,expid,c,plotting=True): + ident = mxd.ident + outdir = mxd.outdir + if mxd.K is None: + if mxd.combinations: + inputs = expid.split('.') + found = False + for invars in itertools.permutations(inputs): + try: + netFile = outdir+"/"+".".join(invars)+'_Tau.net' + mxd.loadnet(netFile) + found = True + break + except: + pass + if found: break + + if not found: + print '{} not found. Need to train this combinatin of inputs'.format(netFile) + raise + else: + invars = mxd.comblist[0] + netFile = outdir+"/"+".".join(invars)+'_Tau.net' + + mxd.net = mxd.loadnet(netFile) + mxd.Input = mxd.comblist[c] + TestStats(mxd,mxd.K,c) + if plotting: + if mxd.angstrom: + make_plots_angstrom(mxd,expid,ident,I=mxd.iTest) + else: + make_plots(mxd,expid,ident,I=mxd.iTest) + else: + k = 1 + for iTrain, iTest in mxd.kf: + I = arange(mxd.nobs) + iValid = I[mxd.iValid] + mxd.iTrain = iValid[iTrain] + mxd.iTest = iValid[iTest] + + if mxd.combinations: + inputs = expid.split('.') + found = False + for invars in itertools.permutations(inputs): + try: + netFile = outdir+"/"+".".join(invars)+'.k={}_Tau.net'.format(str(k)) + mxd.loadnet(netFile) + found = True + print 'found file',netFile + break + except: + pass + + if not found: + print '{} not found. Need to train this combinatin of inputs'.format(netFile) + raise + else: + invars = mxd.comblist[0] + netFile = outdir+"/"+".".join(invars)+'.k={}_Tau.net'.format(str(k)) + + + mxd.net = mxd.loadnet(netFile) + mxd.Input = mxd.comblist[c] + TestStats(mxd,k-1,c) + if plotting: + if mxd.angstrom: + make_plots_angstrom(mxd,expid,'.k={}'.format(str(k)),I=mxd.iTest) + else: + make_plots(mxd,expid,'.k={}'.format(str(k)),I=mxd.iTest) + k = k + 1 + +#--------------------------------------------------------------------- +def _trainMODIS(mxdx): + + if not mxdx.combinations: + Input = mxdx.comblist[0] + _train(mxdx,'.'.join(Input),0) + else: + for c,Input in enumerate(mxdx.comblist): + _train(mxdx,'.'.join(Input),c) + +# ------------------------------------------------------------------- +def _testMODIS(mxdx): + + if not mxdx.combinations: + _test(mxdx,mxdx.expid,0,plotting=True) + else: + for c,Input in enumerate(mxdx.comblist): + _test(mxdx,'.'.join(Input),c,plotting=False) + + +#--------------------------------------------------------------------- +def get_combinations(Input_nnr,Input_const): + comblist = [] + combgroups = [] + for n in arange(len(Input_nnr)): + for invars in itertools.combinations(Input_nnr,n+1): + b = () + for c in invars: + if type(c) is list: + b = b + tuple(c) + else: + b = b + (c,) + #don't do both kinds of abledo together + if not (('BRF' in b) and ('albedo' in b)): + if Input_const is not None: + comblist.append(tuple(Input_const)+b) + combgroups.append((Input_const,)+invars) + else: + comblist.append(b) + combgroups.append(invars) + + if Input_const is not None: + comblist.insert(0,tuple(Input_const)) + combgroups.insert(0,(Input_const,)) + + return comblist,combgroups +#--------------------------------------------------------------------- +def flatten_list(Input_nnr): + Input = () + for i in Input_nnr: + if type(i) is list: + Input = Input + tuple(i) + else: + Input = Input + (i,) + + return list(Input) + +#------------------------------------------------------------------ + +if __name__ == "__main__": + + """ + Example Training/testing + """ + from abc_c6_aux import SummarizeCombinations, SummaryPDFs + + nHidden = None + nHLayers = 1 + combinations = True + Target = ['aTau550',] + Albedo = ['CoxMunkBRF'] + K = 2 + + filename = '/nobackup/6/NNR/Training/giant_C6_10km_Terra_20150921.nc' + retrieval = 'VSDT_OCEAN' + + + doTrain = True + doTest = True + expid = 'CoxMunkTest_wSZA' + Input_const = ['SolarZenith','ScatteringAngle', 'GlintAngle','mRef470','mRef550','mRef660', 'mRef870','mRef1200','mRef1600','mRef2100'] + Input_nnr = ['CoxMunkBRF',['fdu','fcc','fsu']] + aFilter = ['CoxMunkBRF'] + + expid = '{}_{}'.format(retrieval,expid) + + ocean = ABC_Ocean(filename,Albedo=Albedo,verbose=1,aFilter=aFilter) + # Initialize class for training/testing + # --------------------------------------------- + ocean.setupNN(retrieval, expid, + nHidden = nHidden, + nHLayers = nHLayers, + combinations = combinations, + Input_const = Input_const, + Input_nnr = Input_nnr, + Target = Target, + K = K) + + if Input_const is not None: + InputMaster = list((Input_const,) + tuple(Input_nnr)) + else: + InputMaster = Input_nnr + + # Do Training and Testing + # ------------------------ + if doTrain: + _trainMODIS(ocean) + + if doTest: + _testMODIS(ocean) + + if combinations: + SummarizeCombinations(ocean,InputMaster,yrange=None,sortname='slope') + SummaryPDFs(ocean,varnames=['mRef870','mRef660']) diff --git a/src/Components/misc/obs_aod/ABC/giant_viirs.py b/src/Components/misc/obs_aod/ABC/giant_viirs.py new file mode 100644 index 00000000..37ab9dec --- /dev/null +++ b/src/Components/misc/obs_aod/ABC/giant_viirs.py @@ -0,0 +1,608 @@ +#!/bin/env python + +import os +import sys + +from types import * +from netCDF4 import Dataset +from datetime import date, datetime, timedelta +from numpy import savez, meshgrid, array, concatenate, zeros, ones, \ + linspace, sqrt, load, shape, random, interp +import numpy as np +from dateutil.parser import parse as isoparse + +from pyobs.npz import NPZ + +META = ( "Date", + "Time", + "Latitude", + "Longitude", + "SolarZenith", + "SensorZenith", + "ScatteringAngle", + "GlintAngle") + +ANET = ( "mean_AOD0340", + "mean_AOD0380", + "mean_AOD0410", + "mean_AOD0440", + "mean_AOD0470intrp", + "mean_AOD0500", + "mean_AOD0550intrp", + "mean_AOD0660intrp", + "mean_AOD0670", + "mean_AOD0870", + "mean_AOD1020", + "mean_AOD1640", +# "mean_AOD2100intrp", + "mean_WaterVapor", + "nval_AOD0550intrp" + ) + + +xDT_LAND = ("mean_AOD0480corr-l", + "mean_AOD0550corr-l", + "mean_AOD0670corr-l", + "mean_AOD2250corr-l", +# "mean_mref0412-l", +# "mean_mref0443-l", + "mean_mref0480-l", +# "mean_mref0550-l", + "mean_mref0670-l", +# "mean_mref0745-l", +# "mean_mref0870-l", +# "mean_mref1200-l", +# "mean_mref1600-l", + "mean_mref2250-l", + "mean_surfre0480-l", + "mean_surfre0670-l", + "mean_surfre2250-l", + "mean_acfrac-l", + "mode_QA-l", +# "cval_cldpixdistavg", +# "nval_npu0550-l", +# "nval_AOD0550corr-l" + ) + +xDT_OCEAN = ( "mean_AOD0480ea-o", + "mean_AOD0550ea-o", + "mean_AOD0670ea-o", + "mean_AOD0860ea-o", + "mean_AOD1240ea-o", + "mean_AOD1600ea-o", + "mean_AOD2250ea-o", + "mean_mref0480-o", + "mean_mref0550-o", + "mean_mref0670-o", + "mean_mref0860-o", + "mean_mref1240-o", + "mean_mref1600-o", + "mean_mref2250-o", +# "mean_wspeed-o", + "mean_acfrac-o", + "mode_QAavg-o", +# "cval_cldpixdistavg", + "nval_AOD0550ea-o" + ) + +xDB_LAND = ( "mean_AOD0412dpbl-l", + "mean_AOD0488dpbl-l", +# "mean_AOD0550bestdpbl-l", + "mean_AOD0550dpbl-l", + "mean_AOD0670dpbl-l", + "mean_cfracdpbl-l", + "mean_mref0412dpbl-l", + "mean_mref0488dpbl-l", + "mean_mref0670dpbl-l", + "mean_mref0865dpbl-l", + "mean_mref1240dpbl-l", + "mean_mref1640dpbl-l", + "mean_mref2250dpbl-l", + "mean_surfre0412dpbl-l", + "mean_surfre0488dpbl-l", + "mean_surfre0670dpbl-l", + "mode_QAdpbl-l", +# "cval_cldpixdistavg", + "nval_AOD0550dpbl-l", + ) + +xDB_OCEAN = ("mean_AOD0488dpbl-o", +# "mean_AOD0550bestdpbl-l", + "mean_AOD0550dpbl-o", + "mean_AOD0670dpbl-o", + "mean_AOD0865dpbl-o", + "mean_AOD01240dpbl-o", + "mean_AOD01640dpbl-o", + "mean_AOD02250dpbl-o", + "mean_cfracdpbl-o", + "mean_mref0412dpbl-o", + "mean_mref0488dpbl-o", + "mean_mref0670dpbl-o", + "mean_mref0865dpbl-o", + "mean_mref1240dpbl-o", + "mean_mref1640dpbl-o", + "mean_mref2250dpbl-o", + "mode_QAdpbl-o", +# "cval_cldpixdistavg", + ) + +ALIAS = { + "Latitude" : 'lat', + "Longitude" : 'lon', + "mean_AOD0440" : 'aTau440', + "mean_AOD0470intrp" : 'aTau470', + "mean_AOD0500" : 'aTau500', + "mean_AOD0550intrp" : 'aTau550', + "mean_AOD0660intrp" : 'aTau660', + "mean_AOD0870" : 'aTau870', + "mean_AOD2100intrp" : 'aTau2100', + "mean_WaterVapor" : 'aWaterVapor', + "nval_AOD0550intrp" : 'aNcollo', + + "mean_AOD0470corr-l" : 'mTau470', + "mean_AOD0550corr-l" : 'mTau550', + "mean_AOD0660corr-l" : 'mTau660', + "mean_AOD2100corr-l" : 'mTau2100', + "mean_mref0412-l" : 'mRef412', + "mean_mref0443-l" : 'mRef440', + "mean_mref0470-l" : 'mRef470', + "mean_mref0550-l" : 'mRef550', + "mean_mref0660-l" : 'mRef660', + "mean_mref0745-l" : 'mRef745', + "mean_mref0870-l" : 'mRef870', + "mean_mref1200-l" : 'mRef1200', + "mean_mref1600-l" : 'mRef1600', + "mean_mref2100-l" : 'mRef2100', + "mean_surfre0470-l" : 'mSre470', + "mean_surfre0660-l" : 'mSre660', + "mean_surfre2100-l" : 'mSre2100', + "mean_acfrac-l" : 'cloud', + "mode_QA-l" : 'qa', + "nval_AOD0550corr-l" : 'mNcollo', + + "mean_AOD0470ea-o" : 'mTau470', + "mean_AOD0550ea-o" : 'mTau550', + "mean_AOD0660ea-o" : 'mTau660', + "mean_AOD0870ea-o" : 'mTau870', + "mean_AOD1200ea-o" : 'mTau1200', + "mean_AOD1600ea-o" : 'mTau1600', + "mean_AOD2100ea-o" : 'mTau2100', + "mean_mref0470-o" : 'mRef470', + "mean_mref0550-o" : 'mRef550', + "mean_mref0660-o" : 'mRef660', + "mean_mref0870-o" : 'mRef870', + "mean_mref1200-o" : 'mRef1200', + "mean_mref1600-o" : 'mRef1600', + "mean_mref2100-o" : 'mRef2100', + "mean_wspeed-o" : 'speed', + "mean_acfrac-o" : 'cloud', + "mode_QAavg-o" : 'qa', + "nval_AOD0550ea-o" : 'mNcollo', + + "mean_AOD0412dpbl-l" : 'mTau412', + "mean_AOD0470dpbl-l" : 'mTau470', + "mean_AOD0550dpbl-l" : 'mTau550', + "mean_AOD0660dpbl-l" : 'mTau660', + "mean_cfracdpbl-l" : 'cloud', + "mean_mref0412dpbl-l" : 'mRef412', + "mean_mref0470dpbl-l" : 'mRef470', + "mean_mref0660dpbl-l" : 'mRef660', + "mean_surfre0412dpbl-l" : 'mSre412', + "mean_surfre0470dpbl-l" : 'mSre470', + "mean_surfre0660dpbl-l" : 'mSre660', + "mode_QAdpbl-l" : 'qa', + "nval_AOD0550dpbl-l" : 'mNcollo', + + "mean_dtdbAOD0550" : 'mTau550comb', + "cval_cldpixdistavg" : 'clDist' + + + + } + +MISSING = 1.E20 + +# .................................................................................. +class CX_ALBEDO(object): + """ + Container for Cox-Munk albedo + """ + def __init__(self,s_channels): + for ch in s_channels: + self.__dict__['CxAlbedo' + ch] = [] + + +class GIANT(object): + + """ + Read Level-2 GIANT Co-located AERONET/MODIS/VIIRS aerosol files from + GSFC MODIS group. + """ + + + def __init__ (self,filename,xVars=(),only_good=True,tymemax=None): + """ + Creates an GIANT object defining the attributes corresponding + to the SDS of interest. + """ + + if 'SNPP' in filename: self.sat = 'SNPP' + else: self.sat = 'Unknown' + + self.only_good = only_good + Names = META+ANET+xVars + + # Simplify variable names + # ------------------------ + self.ALIAS = ALIAS.copy() + + # Read in variables + # ----------------- + print 'filename ',filename + nc = Dataset(filename) + Alias = self.ALIAS.keys() + self.giantList =[] + for name in Names: + data = nc.variables[name][:] + if name in Alias: + name = self.ALIAS[name] + + # old files use -9999.0 for fill value + # new files use masked arrays + # convert everythong to regular array filling with -9999.0 + # make sure _fill_value is -9999.0 + self.__dict__[name] = np.array(data) + self.giantList.append(name) + nc.close() + + # Form python tyme + # ---------------- + # new files have an ISO_DateTime variable + nc = Dataset(filename) + if 'ISO_DateTime' in nc.variables.keys(): + try: + iso = nc.variables['ISO_DateTime'][:] + self.tyme = array([isoparse(''.join(array(t))) for t in iso]) + except: + # old file only have Date and Time variables + D = self.Date[:,0:10] + T = self.Time[:,0:5] # they didn't save the seconds + # Bug in dataset, first field is blank + D[0] = D[1] + T[0] = T[1] + self.aTau550[0] = -9999.0 + self.tyme = array([ isoparse(''.join(D[i])+'T'+''.join(t)) for i, t in enumerate(T) ]) + nc.close() + + # Limit to the MERRA-2 time series + #--------------------------------- + if tymemax is not None: + tymemax = isoparse(tymemax) + I = self.tyme < tymemax + + for name in Names: + if name in Alias: + name = self.ALIAS[name] + self.__dict__[name] = self.__dict__[name][I] + + self.tyme = self.tyme[I] + + del self.Date, self.Time + self.giantList.remove('Date') + self.giantList.remove('Time') + self.giantList.append('tyme') + + # Record number of observations + # ----------------------------- + self.nobs = len(self.lon) + +#-- + def balance(self,N): + """ + Return indices of observations so that each species does not have more than + N observations. This is meant to be performed with a reduced dataset. + """ + I = zeros(self.lon.shape).astype(bool) + random.seed(32768) # so that we get the same permutation + + for f in (self.fdu,self.fss,self.fcc,self.fsu): + + J = f>0.5 # all obs for which species dominate + n = len(self.lon[J]) # no. obs for this species + P = random.permutation(n) # randomize obs for this species + m = min(n,N) # keep this many + + K = I[J] + K[P[0:m]] = True + I[J] = K + + return I + +#--- + def reduce(self,I): + """ + Reduce observations according to index I. + """ + for name in self.giantList: + q = self.__dict__[name] + print "{} Reducing "+name,q.shape + self.__dict__[name] = q[I] + + self.nobs = len(self.lon) + + + def getCoxMunk(self,filename='/nobackup/NNR/Misc/coxmunk_lut.npz',channel=550.): + """ + Returns ocean albedo as a function of wind speed from look up table. + """ + + # Get precomputed albedo LUT + # -------------------------- + lut = NPZ(filename) + + # Trimmed wind speed + # ------------------ + w10m = self.wind.copy() + w10m[w10m<0] = 0 + w10m[w10m>50.] = 50. + + j = list(lut.channels).index(channel) + + # Interpolate albedo + # ------------------ + albedo = zeros(len(w10m)) + albedo[:] = interp(w10m,lut.speed,lut.albedo[:,j]) + + self.CoxMunkLUT = albedo + + + def calcCoxMunk(self,windFile,channels=[470. ,550. ,660. ,870. ,1200.,1600.,2100.],npzFile=None): + """ + Calls VLIDORT wrapper to calculate CoxMunk Bidirectional Surface Reflectance. + """ + import VLIDORT_BRDF_ABC_ + + channeli = [200,250,300,337,400,488,515,550,633,694,860,1060,1300,1536,1800,2000,2250,2500] + # refractive index + mr = [1.396,1.362,1.349,1.345,1.339,1.335,1.334,1.333,1.332,1.331,1.329,1.326, + 1.323,1.318,1.312,1.306,1.292,1.261] + + wind = np.load(windFile) + u10m = wind['u10m'].astype('float64') + v10m = wind['v10m'].astype('float64') + sza = self.SolarZenith.astype('float64') + vza = self.SensorZenith.astype('float64') + + # needs to be photon travel direction. raa = 0 is forward scattering + saa = self.SolarAzimuth.astype('float64') + saa = saa + 180.0 + saa[saa>=360.0] = saa[saa>=360.0] - 360.0 + raa = np.abs(self.SensorAzimuth - saa).astype('float64') + + if type(channels) is list: + channels = np.array(channels) + + try: + some_object_iterator = iter(channels) + except: + channels = np.array([channels]) + + strch = [str(int(ch)) for ch in channels] + m = interp(channels,channeli,mr) + + albedos, rc = VLIDORT_BRDF_ABC_.coxmunk(1,channels,u10m,v10m,m,sza,raa,vza,-999,1) + + self.sample = CX_ALBEDO(strch) + for i,ch in enumerate(strch): + self.sample.__dict__['CxAlbedo' + ch] = albedos[:,i] + + + if npzFile is not None: + savez(npzFile,**self.sample.__dict__) + +#--- + def speciate(self,aer_x,FineMode=False,Verbose=False): + """ + Use GAAS to derive fractional composition. + """ + + self.sampleFile(aer_x,onlyVars=('TOTEXTTAU', + 'DUEXTTAU', + 'SSEXTTAU', + 'BCEXTTAU', + 'OCEXTTAU', + 'SUEXTTAU', + ),Verbose=Verbose) + s = self.sample + I = (s.TOTEXTTAU<=0) + s.TOTEXTTAU[I] = 1.E30 + self.fdu = s.DUEXTTAU / s.TOTEXTTAU + self.fss = s.SSEXTTAU / s.TOTEXTTAU + self.fbc = s.BCEXTTAU / s.TOTEXTTAU + self.foc = s.OCEXTTAU / s.TOTEXTTAU + self.fcc = self.fbc + self.foc + self.fsu = s.SUEXTTAU / s.TOTEXTTAU + + if FineMode: + TOTEXTTAU = s.TOTEXTTAU[:] + self.sampleFile(aer_x,onlyVars=('DUEXTTFM','SSEXTTFM'),Verbose=Verbose) + self.fduf = s.DUEXTTFM / TOTEXTTAU + self.fssf = s.SSEXTTFM / TOTEXTTAU + + del self.sample + +#--- + def sampleG5(self,gas_x=None,avk_x=None,int_x=None,slv_x=None,ext_Nc=None): + """ + Sample key parameters from GAAS files. + """ + from gfio import GFIOHandle + + if gas_x is not None: + self.sampleFile(gas_x,onlyVars=('AODANA',)) + self.tau_550 = self.sample.AODANA[:] + + if avk_x is not None: + tyme = self.tyme[:] + self.tyme = getSyn(tyme) + self.sampleFile(avk_x,onlyVars=('AOD',)) + self.avk = self.sample.AOD[:] + self.tyme[:] = tyme[:] + + if int_x is not None: + try: + self.sampleFile(int_x,onlyVars=('TQV',)) # As in file spec + self.tpw = self.sample.TQV[:] + except: + self.sampleFile(int_x,onlyVars=('TPW',)) # Larry's name + self.tpw = self.sample.TPW[:] + + if slv_x is not None: + self.sampleFile(slv_x,onlyVars=('U10M','V10M')) + self.wind = sqrt(self.sample.U10M[:]**2 + self.sample.V10M[:]**2) + + if ext_Nc is not None: + self.sampleFile(ext_Nc,onlyVars=('taod',)) + self.tau_660 = self.sample.taod[:,5] # 660 + + + del self.sample + + + def sampleMERRA(self,slv_x='tavg1_2d_slv_Nx',aer_x='tavg1_2d_aer_Nx', + FineMode=False,npzFile=None,Verbose=False): + self.sampleFile(slv_x,onlyVars=('U10M','V10M'), Verbose=Verbose) + self.u10m = self.sample.U10M + self.v10m = self.sample.V10M + self.wind = sqrt(self.sample.U10M[:]**2 + self.sample.V10M[:]**2) + + del self.sample + + self.speciate(aer_x,FineMode=FineMode,Verbose=Verbose) + + if npzFile is not None: + if FineMode: + savez(npzFile,wind=self.wind,u10m=self.u10m,v10m=self.v10m, + fdu=self.fdu,fss=self.fss,fcc=self.fcc,fsu=self.fsu, + fduf=self.fduf,fssf=self.fssf) + else: + savez(npzFile,wind=self.wind,u10m=self.u10m,v10m=self.v10m, + fdu=self.fdu,fss=self.fss,fcc=self.fcc,fsu=self.fsu) + + def sampleMCD43C(self,npzFile=None,Verbose=False): + from mcd43c import MCD43C + + brdf = MCD43C() + brdf.sample(self,Verbose=Verbose) + + if npzFile is not None: + savez(npzFile,**self.brdf.__dict__) + +#--- + def sampleFile(self, inFile, npzFile=None, onlyVars=None, Verbose=False, clmYear=None): + """ + Interpolates all variables of inFile and optionally + save them to file *npzFile* + """ + from gfio import GFIO, GFIOctl, GFIOHandle + + # Instantiate grads and open file + # ------------------------------- + name, ext = os.path.splitext(inFile) + if ext in ( '.nc4', '.nc', '.hdf'): + fh = GFIO(inFile) # open single file + if fh.lm == 1: + timeInterp = False # no time interpolation in this case + else: + raise ValueError, "cannot handle files with more tha 1 time, use ctl instead" + else: + fh = GFIOctl(inFile) # open timeseries + timeInterp = True # perform time interpolation + + self.sample = GFIOHandle(inFile) + if onlyVars is None: + onlyVars = fh.vname + + nt = self.lon.shape + lons = self.lon + lats = self.lat + + if clmYear is None: + tymes = self.tyme + else: + tymes = array([t + timedelta(days=365*(clmYear-t.year)) for t in self.tyme]) + + print 'trange',tymes.min(),tymes.max() + # Loop over variables on file + # --------------------------- + for v in onlyVars: + if Verbose: + print "<> Sampling ", v + if timeInterp: + var = fh.sample(v,lons,lats,tymes,Verbose=Verbose) + else: + var = fh.interp(v,lons,lats) + if len(var.shape) == 1: + self.sample.__dict__[v] = var + elif len(var.shape) == 2: + var = var.T # shape should be (nobs,nz) + self.sample.__dict__[v] = var + else: + raise IndexError, 'variable <%s> has rank = %d'%(v,len(var.shape)) + + if npzFile is not None: + savez(npzFile,**self.sample.__dict__) + + def sampleLoadz(self,npzFile): + """ + Loads sample from npz file. + """ + from grads.gahandle import GaHandle + self.sample = GaHandle(npzFile) + npz = load(npzFile) + for v in npz.keys(): + self.sample.__dict__[v] = npz[v] + +#--- + +class DT_LAND(GIANT): + def __init__(self,filename,tymemax=None): + GIANT.__init__(self,filename,xVars=xDT_LAND,tymemax=tymemax) + if self.sat == 'SNPP': + self.ident = 'vsdtl' + self.ident = self.ident + '_'+ filename.split('/')[-1].split('.')[0] + self.surface = 'land' + +class DT_OCEAN(GIANT): + def __init__(self,filename,tymemax=None): + GIANT.__init__(self,filename,xVars=xDT_OCEAN,tymemax=tymemax) + if self.sat == 'SNPP': + self.ident = 'vsdto' + self.ident = self.ident + '_' + filename.split('/')[-1].split('.')[0] + self.surface = 'ocean' + +class DB_OCEAN(GIANT): + def __init__(self,filename,tymemax=None): + GIANT.__init__(self,filename,xVars=xDB_OCEAN,tymemax=tymemax) + if self.sat == 'SNPP': + self.ident = 'vsdto' + self.ident = self.ident + '_' + filename.split('/')[-1].split('.')[0] + self.surface = 'ocean' + +class DB_LAND(GIANT): + def __init__(self,filename,tymemax=None): + GIANT.__init__(self,filename,xVars=xDB_LAND,tymemax=tymemax) + if self.sat == 'SNPP': + self.ident = 'vsdto' + self.ident = self.ident + '_' + filename.split('/')[-1].split('.')[0] + self.surface = 'land' + + +# ....................................................................................................... + +if __name__ == "__main__": + + # lnd = LAND('/nobackup/6/NNR/Training/giant_C6_10km_Terra_05Oct2015.nc') + ocn = DT_OCEAN('/nobackup/6/NNR/Training/giant_C6_10km_Terra_20150921.nc') + #xblu = DEEP('/Users/adasilva/workspace.local/Data_Analysis/C6/giant_C6_10km_9April2015.nc') + From c55d9a1137da4427c2571fac57e98b711bdab9d3 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 13 Oct 2022 10:33:52 -0400 Subject: [PATCH 09/67] PC - fix dbo variables names in giant_viirs --- src/Components/misc/obs_aod/ABC/giant_viirs.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/giant_viirs.py b/src/Components/misc/obs_aod/ABC/giant_viirs.py index 37ab9dec..91c23b74 100644 --- a/src/Components/misc/obs_aod/ABC/giant_viirs.py +++ b/src/Components/misc/obs_aod/ABC/giant_viirs.py @@ -111,9 +111,9 @@ "mean_AOD0550dpbl-o", "mean_AOD0670dpbl-o", "mean_AOD0865dpbl-o", - "mean_AOD01240dpbl-o", - "mean_AOD01640dpbl-o", - "mean_AOD02250dpbl-o", + "mean_AOD1240dpbl-o", + "mean_AOD1640dpbl-o", + "mean_AOD2250dpbl-o", "mean_cfracdpbl-o", "mean_mref0412dpbl-o", "mean_mref0488dpbl-o", From a090ccb916ec0556a9814cc9e188eace14d3df00 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 28 Oct 2022 12:12:32 -0400 Subject: [PATCH 10/67] PC - fix aliasing in giant_viirs, also added angstrom interpolation of retrieved AOD to MODIS wavelengths to make comparisons easier --- .../misc/obs_aod/ABC/giant_viirs.py | 153 +++++++++++++----- 1 file changed, 109 insertions(+), 44 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/giant_viirs.py b/src/Components/misc/obs_aod/ABC/giant_viirs.py index 91c23b74..13439ee1 100644 --- a/src/Components/misc/obs_aod/ABC/giant_viirs.py +++ b/src/Components/misc/obs_aod/ABC/giant_viirs.py @@ -129,75 +129,94 @@ ALIAS = { "Latitude" : 'lat', "Longitude" : 'lon', + "mean_AOD0340" : 'aTau340', + "mean_AOD0380" : 'aTau380', + "mean_AOD0410" : 'aTau412', "mean_AOD0440" : 'aTau440', "mean_AOD0470intrp" : 'aTau470', "mean_AOD0500" : 'aTau500', "mean_AOD0550intrp" : 'aTau550', "mean_AOD0660intrp" : 'aTau660', + "mean_AOD0670" : 'aTau670', "mean_AOD0870" : 'aTau870', - "mean_AOD2100intrp" : 'aTau2100', + "mean_AOD1020" : 'aTau1020', + "mean_AOD1640" : 'aTau1640', +# "mean_AOD2100intrp" : 'aTau2100', "mean_WaterVapor" : 'aWaterVapor', - "nval_AOD0550intrp" : 'aNcollo', - "mean_AOD0470corr-l" : 'mTau470', + "mean_AOD0480corr-l" : 'mTau480', "mean_AOD0550corr-l" : 'mTau550', - "mean_AOD0660corr-l" : 'mTau660', - "mean_AOD2100corr-l" : 'mTau2100', - "mean_mref0412-l" : 'mRef412', - "mean_mref0443-l" : 'mRef440', - "mean_mref0470-l" : 'mRef470', - "mean_mref0550-l" : 'mRef550', - "mean_mref0660-l" : 'mRef660', - "mean_mref0745-l" : 'mRef745', - "mean_mref0870-l" : 'mRef870', - "mean_mref1200-l" : 'mRef1200', - "mean_mref1600-l" : 'mRef1600', - "mean_mref2100-l" : 'mRef2100', - "mean_surfre0470-l" : 'mSre470', - "mean_surfre0660-l" : 'mSre660', - "mean_surfre2100-l" : 'mSre2100', + "mean_AOD0670corr-l" : 'mTau670', + "mean_AOD2250corr-l" : 'mTau2250', +# "mean_mref0412-l" : 'mRef412', +# "mean_mref0443-l" : 'mRef440', + "mean_mref0480-l" : 'mRef480', +# "mean_mref0550-l" : 'mRef550', + "mean_mref0670-l" : 'mRef670', +# "mean_mref0745-l" : 'mRef745', +# "mean_mref0870-l" : 'mRef870', +# "mean_mref1200-l" : 'mRef1200', +# "mean_mref1600-l" : 'mRef1600', + "mean_mref2250-l" : 'mRef2250', + "mean_surfre0480-l" : 'mSre480', + "mean_surfre0670-l" : 'mSre670', + "mean_surfre2250-l" : 'mSre2250', "mean_acfrac-l" : 'cloud', "mode_QA-l" : 'qa', - "nval_AOD0550corr-l" : 'mNcollo', - "mean_AOD0470ea-o" : 'mTau470', + "mean_AOD0480ea-o" : 'mTau480', "mean_AOD0550ea-o" : 'mTau550', - "mean_AOD0660ea-o" : 'mTau660', - "mean_AOD0870ea-o" : 'mTau870', - "mean_AOD1200ea-o" : 'mTau1200', + "mean_AOD0670ea-o" : 'mTau670', + "mean_AOD0860ea-o" : 'mTau860', + "mean_AOD1240ea-o" : 'mTau1240', "mean_AOD1600ea-o" : 'mTau1600', - "mean_AOD2100ea-o" : 'mTau2100', - "mean_mref0470-o" : 'mRef470', + "mean_AOD2250ea-o" : 'mTau2250', + "mean_mref0480-o" : 'mRef470', "mean_mref0550-o" : 'mRef550', - "mean_mref0660-o" : 'mRef660', - "mean_mref0870-o" : 'mRef870', - "mean_mref1200-o" : 'mRef1200', + "mean_mref0670-o" : 'mRef670', + "mean_mref0860-o" : 'mRef860', + "mean_mref1240-o" : 'mRef1240', "mean_mref1600-o" : 'mRef1600', - "mean_mref2100-o" : 'mRef2100', - "mean_wspeed-o" : 'speed', + "mean_mref2250-o" : 'mRef2250', "mean_acfrac-o" : 'cloud', "mode_QAavg-o" : 'qa', - "nval_AOD0550ea-o" : 'mNcollo', "mean_AOD0412dpbl-l" : 'mTau412', - "mean_AOD0470dpbl-l" : 'mTau470', + "mean_AOD0488dpbl-l" : 'mTau488', "mean_AOD0550dpbl-l" : 'mTau550', - "mean_AOD0660dpbl-l" : 'mTau660', + "mean_AOD0670dpbl-l" : 'mTau670', + "mean_AOD2250dpbl-l" : 'mTau2250', "mean_cfracdpbl-l" : 'cloud', "mean_mref0412dpbl-l" : 'mRef412', - "mean_mref0470dpbl-l" : 'mRef470', - "mean_mref0660dpbl-l" : 'mRef660', + "mean_mref0488dpbl-l" : 'mRef488', + "mean_mref0550dpbl-l" : 'mRef550', + "mean_mref0670dpbl-l" : 'mRef670', + "mean_mref0865dpbl-l" : 'mRef865', + "mean_mref1240dpbl-l" : 'mRef1240', + "mean_mref1640dpbl-l" : 'mRef1640', + "mean_mref2250dpbl-l" : 'mRef2250', "mean_surfre0412dpbl-l" : 'mSre412', - "mean_surfre0470dpbl-l" : 'mSre470', - "mean_surfre0660dpbl-l" : 'mSre660', + "mean_surfre0488dpbl-l" : 'mSre488', + "mean_surfre0670dpbl-l" : 'mSre670', "mode_QAdpbl-l" : 'qa', - "nval_AOD0550dpbl-l" : 'mNcollo', - - "mean_dtdbAOD0550" : 'mTau550comb', - "cval_cldpixdistavg" : 'clDist' - - + "mean_AOD0488dpbl-o" : 'mTau488', + "mean_AOD0550dpbl-o" : 'mTau550', + "mean_AOD0670dpbl-o" : 'mTau670', + "mean_AOD0865dpbl-o" : 'mTau865', + "mean_AOD1240dpbl-o" : 'mTau1240', + "mean_AOD1640dpbl-o" : 'mTau1640', + "mean_AOD2250dpbl-o" : 'mTau2250', + "mean_cfracdpbl-o" : 'cloud', + "mean_mref0412dpbl-o" : 'mRef412', + "mean_mref0488dpbl-o" : 'mRef488', + "mean_mref0550dpbl-o" : 'mRef550', + "mean_mref0670dpbl-o" : 'mRef670', + "mean_mref865dpbl-o" : 'mRef865', + "mean_mref1240dpbl-o" : 'mRef1240', + "mean_mref1640dpbl-o" : 'mRef1640', + "mean_mref2250dpbl-o" : 'mRef2250', + "mode_QAdpbl-o" : 'qa', } MISSING = 1.E20 @@ -294,7 +313,7 @@ def __init__ (self,filename,xVars=(),only_good=True,tymemax=None): # Record number of observations # ----------------------------- - self.nobs = len(self.lon) + self.nobs = len(self.lon) #-- def balance(self,N): @@ -573,6 +592,14 @@ def __init__(self,filename,tymemax=None): self.ident = self.ident + '_'+ filename.split('/')[-1].split('.')[0] self.surface = 'land' + # Angstrom interpolate retrieved AOD to + # nominal MODIS wavelengths (470, 550, 660) + # for ease of comparison + # VIIRS channels = 480., 550., 670., 2250. + # ------------------------------------------------ + self.mTau470 = aodInterpAngs(470.,self.mTau480,self.mTau550,480.,550.) + self.mTau660 = aodInterpAngs(660.,self.mTau550,self.mTau670,550.,670.) + class DT_OCEAN(GIANT): def __init__(self,filename,tymemax=None): GIANT.__init__(self,filename,xVars=xDT_OCEAN,tymemax=tymemax) @@ -581,6 +608,16 @@ def __init__(self,filename,tymemax=None): self.ident = self.ident + '_' + filename.split('/')[-1].split('.')[0] self.surface = 'ocean' + # Angstrom interpolate retrieved AOD to + # nominal MODIS wavelengths (470, 550, 660, 870, 1200, 1600, 2100) + # for ease of comparison + # VIIRS channels = 480., 550., 670., 860., 1240., 1600., 2250. + # ------------------------------------------------ + self.mTau470 = aodInterpAngs(470.,self.mTau480,self.mTau550,480.,550.) + self.mTau660 = aodInterpAngs(660.,self.mTau550,self.mTau670,550.,670.) + self.mTau870 = aodInterpAngs(870.,self.mTau670,self.mTau860,670.,860.) + self.mTau2100 = aodInterpAngs(2100.,self.mTau1600,self.mTau2250,1600.,2250.) + class DB_OCEAN(GIANT): def __init__(self,filename,tymemax=None): GIANT.__init__(self,filename,xVars=xDB_OCEAN,tymemax=tymemax) @@ -589,6 +626,16 @@ def __init__(self,filename,tymemax=None): self.ident = self.ident + '_' + filename.split('/')[-1].split('.')[0] self.surface = 'ocean' + # Angstrom interpolate retrieved AOD to + # nominal MODIS wavelengths (470, 550, 660, 870, 1200, 1600, 2100) + # for ease of comparison + # VIIRS channels = 488., 550., 670., 865., 1240., 1640., 2250. + # ------------------------------------------------ + self.mTau470 = aodInterpAngs(470.,self.mTau488,self.mTau550,488.,550.) + self.mTau660 = aodInterpAngs(660.,self.mTau550,self.mTau670,550.,670.) + self.mTau870 = aodInterpAngs(870.,self.mTau670,self.mTau865,670.,865.) + self.mTau2100 = aodInterpAngs(2100.,self.mTau1640,self.mTau2250,1640.,2250.) + class DB_LAND(GIANT): def __init__(self,filename,tymemax=None): GIANT.__init__(self,filename,xVars=xDB_LAND,tymemax=tymemax) @@ -597,7 +644,25 @@ def __init__(self,filename,tymemax=None): self.ident = self.ident + '_' + filename.split('/')[-1].split('.')[0] self.surface = 'land' + # Angstrom interpolate retrieved AOD to + # nominal MODIS wavelengths (412, 470, 550, 660) + # for ease of comparison + # VIIRS channels = 412, 488, 550, 670 + # ------------------------------------------------ + self.mTau470 = aodInterpAngs(470.,self.mTau412,self.mTau488,412.,488.) + self.mTau660 = aodInterpAngs(660.,self.mTau550,self.mTau670,550.,670.) + +#--- +def aodInterpAngs(lambda_,tau1,tau2,lambda1,lambda2): + """ + Angstrom-interpolated AOD. + """ + I = (tau1 > 0.0) & (tau2 > 0.0) + angstrom = -np.log(tau1[I]/tau2[I])/np.log(lambda1/lambda2) + tau = -9999. * np.ones(len(tau1)) + tau[I] = tau2[I] * (lambda2/lambda_)**angstrom + return tau # ....................................................................................................... if __name__ == "__main__": From 90db6c6cf7b9b4166e1ebb74de14f9682ae4a88e Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 28 Oct 2022 15:12:50 -0400 Subject: [PATCH 11/67] PC - fix variable naming in abc_viirs --- src/Components/misc/obs_aod/ABC/abc_viirs.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/abc_viirs.py b/src/Components/misc/obs_aod/ABC/abc_viirs.py index a0aae9bf..ff613540 100644 --- a/src/Components/misc/obs_aod/ABC/abc_viirs.py +++ b/src/Components/misc/obs_aod/ABC/abc_viirs.py @@ -340,7 +340,7 @@ def __init__ (self,fname, (self.mTau480 > -0.01) &\ (self.mTau550 > -0.01) &\ (self.mTau670 > -0.01) &\ - (self.mTau870 > -0.01) &\ + (self.mTau860 > -0.01) &\ (self.mRef480 > 0.0) &\ (self.mRef550 > 0.0) &\ (self.mRef670 > 0.0) &\ @@ -428,7 +428,7 @@ def __init__ (self,fname, (self.aTau550 > -0.01) &\ (self.aTau660 > -0.01) &\ (self.aTau870 > -0.01) &\ - (self.mTau480 > -0.01) &\ + (self.mTau488 > -0.01) &\ (self.mTau550 > -0.01) &\ (self.mTau670 > -0.01) &\ (self.mTau865 > -0.01) &\ @@ -598,15 +598,15 @@ def __init__ (self, fname, (self.aTau550 > -0.01) & \ (self.aTau660 > -0.01) & \ (self.mTau412 > -0.01) & \ - (self.mTau480 > -0.01) & \ + (self.mTau488 > -0.01) & \ (self.mTau550 > -0.01) & \ (self.mTau670 > -0.01) & \ - (self.mTau2250 > -0.01) & \ (self.cloud<cloud_thresh) & \ (self.cloud >= 0) & \ (self.ScatteringAngle<170.) & \ (self.mRef412 > 0) & \ (self.mRef488 > 0) & \ + (self.mRef550 > 0) & \ (self.mRef670 > 0) & \ (self.mRef865 > 0) & \ (self.mRef1240 > 0) & \ From 27ccace85885ac1d33b622808297ee6bca949209 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 28 Oct 2022 15:14:05 -0400 Subject: [PATCH 12/67] PC - fixed missing variables in DB products and some variable naming in giant_viirs --- src/Components/misc/obs_aod/ABC/giant_viirs.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/giant_viirs.py b/src/Components/misc/obs_aod/ABC/giant_viirs.py index 13439ee1..ba43a044 100644 --- a/src/Components/misc/obs_aod/ABC/giant_viirs.py +++ b/src/Components/misc/obs_aod/ABC/giant_viirs.py @@ -93,6 +93,7 @@ "mean_cfracdpbl-l", "mean_mref0412dpbl-l", "mean_mref0488dpbl-l", + "mean_mref0550dpbl-l", "mean_mref0670dpbl-l", "mean_mref0865dpbl-l", "mean_mref1240dpbl-l", @@ -117,6 +118,7 @@ "mean_cfracdpbl-o", "mean_mref0412dpbl-o", "mean_mref0488dpbl-o", + "mean_mref0550dpbl-o", "mean_mref0670dpbl-o", "mean_mref0865dpbl-o", "mean_mref1240dpbl-o", @@ -171,7 +173,7 @@ "mean_AOD1240ea-o" : 'mTau1240', "mean_AOD1600ea-o" : 'mTau1600', "mean_AOD2250ea-o" : 'mTau2250', - "mean_mref0480-o" : 'mRef470', + "mean_mref0480-o" : 'mRef480', "mean_mref0550-o" : 'mRef550', "mean_mref0670-o" : 'mRef670', "mean_mref0860-o" : 'mRef860', @@ -212,7 +214,7 @@ "mean_mref0488dpbl-o" : 'mRef488', "mean_mref0550dpbl-o" : 'mRef550', "mean_mref0670dpbl-o" : 'mRef670', - "mean_mref865dpbl-o" : 'mRef865', + "mean_mref0865dpbl-o" : 'mRef865', "mean_mref1240dpbl-o" : 'mRef1240', "mean_mref1640dpbl-o" : 'mRef1640', "mean_mref2250dpbl-o" : 'mRef2250', From 260c7aa3fb8cec2e3c802adffaf5ff5438def68d Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 31 Oct 2022 14:52:03 -0400 Subject: [PATCH 13/67] PC - fixed typo in writeg function of vx04.py --- src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index 71649253..85bf2273 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -779,7 +779,7 @@ def writeg(self,filename=None,dir='.',expid=None,refine=8,res=None, f.write('count_tau', nymd, nhms, binobscnt3d(self.lon,self.lat,aod,im,jm,MISSING) ) f.write('count_tau_', nymd, nhms, - binsobscnt3d(self.lon,self.lat,aod_,im,jm,MISSING) ) + binobscnt3d(self.lon,self.lat,aod_,im,jm,MISSING) ) f.write('cloud', nymd, nhms, binobs2d(self.lon,self.lat,self.cloud,im,jm,MISSING) ) From 8f91a1228c7bae90879291c93d1c99e3e0b97896 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Tue, 1 Nov 2022 18:11:55 -0400 Subject: [PATCH 14/67] PC - DT bug fix. need iGood to filter masked values as well --- src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index 85bf2273..0fc8656e 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -256,9 +256,9 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, # Determine index of "good" observations # -------------------------------------- if self.algo == 'DT_LAND': - self.iGood = self.Land_Ocean_Quality_Flag == BEST + self.iGood = (self.Land_Ocean_Quality_Flag == BEST) & (~self.Corrected_Optical_Depth_Land.mask[:,1]) elif self.algo == 'DT_OCEAN': - self.iGood = self.Land_Ocean_Quality_Flag > BAD + self.iGood = (self.Land_Ocean_Quality_Flag > BAD) & (~self.Effective_Optical_Depth_Average_Ocean.mask[:,1]) elif self.algo == 'DB_LAND': self.iGood = self.Aerosol_Optical_Thickness_QA_Flag_Land > BAD # for now elif self.algo == 'DB_OCEAN': From b55547d09939e1baf8827da3b1c2c1e5039bc58f Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Wed, 2 Nov 2022 10:23:49 -0400 Subject: [PATCH 15/67] PC - for viirs sometimes the aod is angstrom interpolated and not availalbe (e.g. at 470). make sure during validation plotting those points are excluded. --- src/Components/misc/obs_aod/ABC/abc_c6_aux.py | 20 +++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/abc_c6_aux.py b/src/Components/misc/obs_aod/ABC/abc_c6_aux.py index 1cbfe173..906924ed 100644 --- a/src/Components/misc/obs_aod/ABC/abc_c6_aux.py +++ b/src/Components/misc/obs_aod/ABC/abc_c6_aux.py @@ -487,7 +487,11 @@ def make_plots_angstrom(mxd,expid,ident,I=None): name = 'm'+mxd.Target[t][1:] if name in mxd.__dict__: original = log(mxd.__dict__[name][I]+0.01) - fig = _plotKDE(targets[:,t],original,y_label='Original MODIS') + + # protect against some of the othe wavelengths having negative values + ii = mxd.__dict__[name][I] > -0.01 + + fig = _plotKDE(targets[:,t][ii],original[ii],y_label='Original MODIS') title("Log("+mxd.Target[t][1:]+"+0.01)- "+ident) savefig(outdir+"/"+expid+"."+ident+"_kde-"+name[1:]+'.png') plt.close(fig) @@ -556,9 +560,12 @@ def make_plots_angstrom(mxd,expid,ident,I=None): print 'orig t,wav',t,name[4:] wav = float(name[4:]) oo = mxd.__dict__[name][I] + 0.01 # add 0.01 to handle negatives - tt = np.exp(targets[:,t]) # keep + 0.01 to handle negatives - AEo = -1.*np.log(refo/oo)/np.log(refwav/wav) - AEt = -1.*np.log(reft/tt)/np.log(refwav/wav) + # protect against interpolated wavelengths that might have -9999 + ii = oo > 0 + oo = oo[ii] + tt = np.exp(targets[:,t][ii]) # keep + 0.01 to handle negatives + AEo = -1.*np.log(refo[ii]/oo)/np.log(refwav/wav) + AEt = -1.*np.log(reft[ii]/tt)/np.log(refwav/wav) fig = _plotKDE(AEt,AEo,y_label='Original MODIS',x_bins=bins,y_bins=bins) title("AE 550/"+name[3:]) savefig(outdir+"/"+expid+"."+ident+"_kde-AE"+name[4:]+'.png') @@ -604,6 +611,11 @@ def make_error_pdfs(mxd,Input,expid,ident,K=None,I=None,Title=None,netfileRoot=N results = np.log(tau + 0.01) original = mxd.__dict__[name][I[0]] + # Protect against -999 in interpolated values + ii = original > -0.01 + original = original[ii] + targets = targets[ii] + results = results[ii] if mxd.laod: original = log(original + 0.01) From ebebe9b8e8622095bde3568e42ecb47ceb5d880d Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Wed, 2 Nov 2022 16:30:33 -0400 Subject: [PATCH 16/67] PC - add interpolated aods to gianlist so they are reduced with the other variables. these aren't used except for testing/validation --- src/Components/misc/obs_aod/ABC/giant_viirs.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/Components/misc/obs_aod/ABC/giant_viirs.py b/src/Components/misc/obs_aod/ABC/giant_viirs.py index ba43a044..2b733e69 100644 --- a/src/Components/misc/obs_aod/ABC/giant_viirs.py +++ b/src/Components/misc/obs_aod/ABC/giant_viirs.py @@ -602,6 +602,8 @@ def __init__(self,filename,tymemax=None): self.mTau470 = aodInterpAngs(470.,self.mTau480,self.mTau550,480.,550.) self.mTau660 = aodInterpAngs(660.,self.mTau550,self.mTau670,550.,670.) + self.giantList += ['mTau470','mTau660'] + class DT_OCEAN(GIANT): def __init__(self,filename,tymemax=None): GIANT.__init__(self,filename,xVars=xDT_OCEAN,tymemax=tymemax) @@ -618,7 +620,9 @@ def __init__(self,filename,tymemax=None): self.mTau470 = aodInterpAngs(470.,self.mTau480,self.mTau550,480.,550.) self.mTau660 = aodInterpAngs(660.,self.mTau550,self.mTau670,550.,670.) self.mTau870 = aodInterpAngs(870.,self.mTau670,self.mTau860,670.,860.) - self.mTau2100 = aodInterpAngs(2100.,self.mTau1600,self.mTau2250,1600.,2250.) + self.mTau2100 = aodInterpAngs(2100.,self.mTau1600,self.mTau2250,1600.,2250.) + + self.giantList +=['mTau470','mTau660','mTau870','mTau2100'] class DB_OCEAN(GIANT): def __init__(self,filename,tymemax=None): @@ -638,6 +642,8 @@ def __init__(self,filename,tymemax=None): self.mTau870 = aodInterpAngs(870.,self.mTau670,self.mTau865,670.,865.) self.mTau2100 = aodInterpAngs(2100.,self.mTau1640,self.mTau2250,1640.,2250.) + self.giantList += ['mTau470','mTau660','mTau870','mTau2100'] + class DB_LAND(GIANT): def __init__(self,filename,tymemax=None): GIANT.__init__(self,filename,xVars=xDB_LAND,tymemax=tymemax) @@ -654,6 +660,8 @@ def __init__(self,filename,tymemax=None): self.mTau470 = aodInterpAngs(470.,self.mTau412,self.mTau488,412.,488.) self.mTau660 = aodInterpAngs(660.,self.mTau550,self.mTau670,550.,670.) + self.giantList += ['mTau470','mTau660'] + #--- def aodInterpAngs(lambda_,tau1,tau2,lambda1,lambda2): """ From 8926d0b15617d258108d43633e470f979db605fa Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Wed, 2 Nov 2022 16:31:22 -0400 Subject: [PATCH 17/67] PC - fix db ocean qa filter to only be BEST. anything below that includes low confidence inland/turbid waters --- src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index 0fc8656e..acc13f4f 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -262,7 +262,7 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, elif self.algo == 'DB_LAND': self.iGood = self.Aerosol_Optical_Thickness_QA_Flag_Land > BAD # for now elif self.algo == 'DB_OCEAN': - self.iGood = self.Aerosol_Optical_Thickness_QA_Flag_Ocean > BAD + self.iGood = self.Aerosol_Optical_Thickness_QA_Flag_Ocean == BEST else: raise ValueError, 'invalid algorithm (very strange)' From e240e95cd41c008e9201cd6514526089ad29e539 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 17 Nov 2022 13:36:01 -0500 Subject: [PATCH 18/67] PC - when doing angstrom interpolation, keep obs where AOD = 0 --- src/Components/misc/obs_aod/ABC/giant_viirs.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/Components/misc/obs_aod/ABC/giant_viirs.py b/src/Components/misc/obs_aod/ABC/giant_viirs.py index 2b733e69..b612cdec 100644 --- a/src/Components/misc/obs_aod/ABC/giant_viirs.py +++ b/src/Components/misc/obs_aod/ABC/giant_viirs.py @@ -672,6 +672,9 @@ def aodInterpAngs(lambda_,tau1,tau2,lambda1,lambda2): angstrom = -np.log(tau1[I]/tau2[I])/np.log(lambda1/lambda2) tau = -9999. * np.ones(len(tau1)) tau[I] = tau2[I] * (lambda2/lambda_)**angstrom + + I = (tau1 == 0.0) & (tau2 == 0.0) + tau[I] = 0.0 return tau # ....................................................................................................... From 5e020acc6f7dc739c783dce9dc6ae9805a229fb2 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 5 Dec 2022 12:01:43 -0500 Subject: [PATCH 19/67] PC - vx04 fix bug in DT calculation of RAA --- src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index acc13f4f..0b42884b 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -302,11 +302,16 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, # this kind of seems to match the DB RAA # the DT sensor and solar azimuth angles seem off # I can't find a DB definition for RAA so this is close enough - raa = self.SolarZenith - self.SensorZenith - ii = raa <0 - raa[ii] = raa[ii] + 180. - ii = raa < 0 - raa[ii] = raa[ii]*-1. + saa = self.SolarAzimuth + vaa = self.SensorAzimuth + i = saa < 0 + saa[i] = 360. + saa[i] + i = vaa < 0 + vaa[i] = 360. + vaa[i] + raa = saa - vaa + 180. + i = raa > 180. + raa[i] = 360. - raa[i] + self.RelativeAzimuth = raa self.SDS += ('RelativeAzimuth',) From b815eb1a90700300a39956529e3619defcc38942 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 5 Dec 2022 14:54:26 -0500 Subject: [PATCH 20/67] PC - vx04 add some more ancillary variables to DB varlist --- src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index 0b42884b..55cf4895 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -66,7 +66,10 @@ 'Algorithm_Flag_Land', 'Aerosol_Type_Land', 'Number_Of_Pixels_Used_Land', - 'Number_Valid_Pixels'), + 'Number_Valid_Pixels', + 'TOA_NDVI', + 'Total_Column_Ozone', + 'Cell_Average_Elevation_Land'), DB_OCEAN = ('Aerosol_Optical_Thickness_550_Ocean', 'Spectral_Aerosol_Optical_Thickness_Ocean', 'Spectral_TOA_Reflectance_Ocean', @@ -78,7 +81,9 @@ 'Number_Of_Pixels_Used_Ocean', 'Number_Valid_Pixels', 'Wind_Speed', - 'Wind_Direction') + 'Wind_Direction', + 'Cell_Average_Elevation_Ocean', + 'Total_Column_Ozone') ) # NOTE: DEEP BLUE does not have cloud information in their files. @@ -121,13 +126,16 @@ Aerosol_Optical_Thickness_550_Land = 'aod550', Spectral_Surface_Reflectance = 'sfc_reflectance', Spectral_TOA_Reflectance_Land = 'reflectance', + TOA_NDVI = 'NDVI', Spectral_Single_Scattering_Albedo_Land = 'ssa', - Algorithm_Flag_Land = 'atype', + Algorithm_Flag_Land = 'aflag', + Aerosol_Type_Land = 'atype', Angstrom_Exponent_Land = 'angstrom', Spectral_Aerosol_Optical_Thickness_Ocean = 'aod', Aerosol_Optical_Thickness_550_Ocean = 'aod550', Spectral_TOA_Reflectance_Ocean = 'reflectance', - Algorithm_Flag_Ocean = 'atype', + Algorithm_Flag_Ocean = 'aflag', + Aerosol_Type_Ocean = 'atype', Angstrom_Exponent_Ocean = 'angstrom', Number_Of_Pixels_Used_Ocean = 'Number_Of_Pixels_Used', Number_Of_Pixels_Used_Land = 'Number_Of_Pixels_Used', From f50cfbfee3c60bd9a0ef1a39640a36acfbf2f985 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 5 Dec 2022 14:58:25 -0500 Subject: [PATCH 21/67] PC- brdf.py add option to give raa directly. suppots DB viirs product that only report relative azimuth angles --- src/Components/misc/obs_aod/ABC/brdf.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/brdf.py b/src/Components/misc/obs_aod/ABC/brdf.py index bb06f7e8..e2e9ca74 100644 --- a/src/Components/misc/obs_aod/ABC/brdf.py +++ b/src/Components/misc/obs_aod/ABC/brdf.py @@ -7,7 +7,7 @@ from numpy import arctan, sqrt, arccos import numpy as np -def rtlsReflectance(Kiso,Kgeo,Kvol,sza,vza,saa,vaa): +def rtlsReflectance(Kiso,Kgeo,Kvol,sza,vza,saa,vaa,raa=None): """ Lucht et al. (2000) IEEE Transactions on Geoscience and Remote Sensing """ @@ -16,7 +16,8 @@ def rtlsReflectance(Kiso,Kgeo,Kvol,sza,vza,saa,vaa): h_b = 2.0 b_r = 1.0 - raa = vaa - saa # raa = 0 is back scattering direction + if raa is None: + raa = vaa - saa # raa = 0 is back scattering direction # degrees to radians raa = raa*pi/180.0 From 2aed852a6bf387c3a9603979e52b8af2ee0cb613 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 5 Dec 2022 17:27:30 -0500 Subject: [PATCH 22/67] PC - vx04 modify for new path for VIIRS files --- src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index 55cf4895..94afccb3 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -894,7 +894,7 @@ def granules ( path, algo, sat, syn_time, coll='011', nsyn=8, verbose=False ): # Get product name # ----------------- Algo = algo.split('_')[0] - prod = 'AER{}_{}'.format(Algo,sat) + prod = 'AER{}'.format(Algo) # Determine synoptic time range # ----------------------------- @@ -909,8 +909,8 @@ def granules ( path, algo, sat, syn_time, coll='011', nsyn=8, verbose=False ): while t < t2: if t >= t1: doy = t.timetuple()[7] - basen = "%s/%s/Level2/%s/%04d/%03d/AER%s_L2_VIIRS_%s.A%04d%03d.%02d%02d.%s.*.nc"\ - %(path,coll,prod,t.year,doy,Algo,sat,t.year,doy,t.hour,t.minute,coll) + basen = "%s/%s/%s/%s/Level2/%04d/%03d/AER%s_L2_VIIRS_%s.A%04d%03d.%02d%02d.%s.*.nc"\ + %(path,prod,sat,coll,t.year,doy,Algo,sat,t.year,doy,t.hour,t.minute,coll) try: filen = glob(basen)[0] Granules += [filen,] From 89d49c30f6c1f6cf2b4e9b2a4d73d4f2d01d7145 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Wed, 7 Dec 2022 12:51:26 -0500 Subject: [PATCH 23/67] PC - add some additional ancillary variables to DB read --- src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index 94afccb3..3d23d68c 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -69,7 +69,8 @@ 'Number_Valid_Pixels', 'TOA_NDVI', 'Total_Column_Ozone', - 'Cell_Average_Elevation_Land'), + 'Cell_Average_Elevation_Land', + 'Precipitable_Water'), DB_OCEAN = ('Aerosol_Optical_Thickness_550_Ocean', 'Spectral_Aerosol_Optical_Thickness_Ocean', 'Spectral_TOA_Reflectance_Ocean', @@ -83,7 +84,8 @@ 'Wind_Speed', 'Wind_Direction', 'Cell_Average_Elevation_Ocean', - 'Total_Column_Ozone') + 'Total_Column_Ozone', + 'Precipitable_Water') ) # NOTE: DEEP BLUE does not have cloud information in their files. @@ -143,6 +145,7 @@ Aerosol_Optical_Thickness_QA_Flag_Ocean = 'qa_flag', Land_Ocean_Quality_Flag = 'qa_flag', Scan_Start_Time = 'Time', + Cell_Average_Elevation_Land = 'pixel_elevation', ) BAD, MARGINAL, GOOD, BEST = ( 0, 1, 2, 3 ) # DT QA marks From f0623ba04d425b93003b35154d3709ab0e7308fe Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 20 Jan 2023 12:48:40 -0500 Subject: [PATCH 24/67] PC - changes to giant viirs reader to include additional variables: ozone, water vapor, algorithm type, assumed aerosol type, and NDVI --- .../misc/obs_aod/ABC/giant_viirs.py | 21 +++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/src/Components/misc/obs_aod/ABC/giant_viirs.py b/src/Components/misc/obs_aod/ABC/giant_viirs.py index b612cdec..1f94ffc8 100644 --- a/src/Components/misc/obs_aod/ABC/giant_viirs.py +++ b/src/Components/misc/obs_aod/ABC/giant_viirs.py @@ -19,6 +19,7 @@ "Longitude", "SolarZenith", "SensorZenith", + "RelativeAzimuth", "ScatteringAngle", "GlintAngle") @@ -103,8 +104,14 @@ "mean_surfre0488dpbl-l", "mean_surfre0670dpbl-l", "mode_QAdpbl-l", + "mode_algflgdpbl-l", # "cval_cldpixdistavg", "nval_AOD0550dpbl-l", + "mean_NDVIdpbl-l", + "mean_Total_Column_Ozonedpbl-l", + "mode_atypedpbl-l", + "mean_PrecipitableWaterdpbl-l", + "mean_Elevationdpbl-l", ) xDB_OCEAN = ("mean_AOD0488dpbl-o", @@ -125,7 +132,11 @@ "mean_mref1640dpbl-o", "mean_mref2250dpbl-o", "mode_QAdpbl-o", + "mode_algflgdpbl-o", # "cval_cldpixdistavg", + "mean_Total_Column_Ozonedpbl-o", + "mode_atypedpbl-o", + "mean_PrecipitableWaterdpbl-o", ) ALIAS = { @@ -201,6 +212,12 @@ "mean_surfre0488dpbl-l" : 'mSre488', "mean_surfre0670dpbl-l" : 'mSre670', "mode_QAdpbl-l" : 'qa', + "mode_algflgdpbl-l" : 'algflag', + "mode_atypedpbl-l" : 'atype', + "mean_Total_Column_Ozonedpbl-l" : 'colO3', + "mean_NDVIdpbl-l" : 'ndvi', + "mean_PrecipitableWaterdpbl-l" : "water", + "mean_Elevationdpbl-l" : "pixel_elevation", "mean_AOD0488dpbl-o" : 'mTau488', "mean_AOD0550dpbl-o" : 'mTau550', @@ -219,6 +236,10 @@ "mean_mref1640dpbl-o" : 'mRef1640', "mean_mref2250dpbl-o" : 'mRef2250', "mode_QAdpbl-o" : 'qa', + "mode_algflgdpbl-o" : 'algflag', + "mode_atypedpbl-o" : 'atype', + "mean_Total_Column_Ozonedpbl-o" : 'colO3', + "mean_PrecipitableWaterdpbl-o" : 'water', } MISSING = 1.E20 From 9ab11e24a7a508b131af928356367de0cd84a1e9 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 20 Jan 2023 12:56:06 -0500 Subject: [PATCH 25/67] PC - fix in viirs reader for how cloud fraction is calcualted for VIIRS DB. add number of suitable and unsuitable pixels to variables read. --- src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index 3d23d68c..764fdac3 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -67,6 +67,7 @@ 'Aerosol_Type_Land', 'Number_Of_Pixels_Used_Land', 'Number_Valid_Pixels', + 'Unsuitable_Pixel_Fraction_Land_Ocean', 'TOA_NDVI', 'Total_Column_Ozone', 'Cell_Average_Elevation_Land', @@ -81,6 +82,7 @@ 'Aerosol_Type_Ocean', 'Number_Of_Pixels_Used_Ocean', 'Number_Valid_Pixels', + 'Unsuitable_Pixel_Fraction_Land_Ocean', 'Wind_Speed', 'Wind_Direction', 'Cell_Average_Elevation_Ocean', @@ -139,8 +141,9 @@ Algorithm_Flag_Ocean = 'aflag', Aerosol_Type_Ocean = 'atype', Angstrom_Exponent_Ocean = 'angstrom', - Number_Of_Pixels_Used_Ocean = 'Number_Of_Pixels_Used', - Number_Of_Pixels_Used_Land = 'Number_Of_Pixels_Used', + Number_Of_Pixels_Used_Ocean = 'npixels_used', + Number_Of_Pixels_Used_Land = 'npixels_used', + Number_Valid_Pixels = 'npixels_valid', Aerosol_Optical_Thickness_QA_Flag_Land = 'qa_flag', Aerosol_Optical_Thickness_QA_Flag_Ocean = 'qa_flag', Land_Ocean_Quality_Flag = 'qa_flag', @@ -311,8 +314,9 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, self.SDS += ('GlintAngle',) elif 'DT' in self.algo: # this kind of seems to match the DB RAA - # the DT sensor and solar azimuth angles seem off - # I can't find a DB definition for RAA so this is close enough + # the DT sensor and solar azimuth angles have a bug + # this needs to be rechecked with updated version of DT dataset + # DB definition for RAA is 'Gordon Convention', RAA=0 is back scattering saa = self.SolarAzimuth vaa = self.SensorAzimuth i = saa < 0 @@ -378,7 +382,7 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, # Create a pseudo cloud fraction for Deep Blue if Algo == 'DB': - self.cloud = self.Number_Of_Pixels_Used.astype(float)/self.Number_Valid_Pixels.astype(float) + self.cloud = 1. - self.npixels_used.astype(float)/self.npixels_valid.astype(float) #--- def _readList(self,List): From 6f41084f8e82d1700af298a81dc6487f7a2a46b9 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 23 Jan 2023 10:22:40 -0500 Subject: [PATCH 26/67] PC - changes to abc_aux_viirs and nn.py to fix issues that arise from new version sklearn when you have 1 target --- src/Components/misc/obs_aod/ABC/abc_c6_aux.py | 16 +++++++++++----- src/Components/misc/obs_aod/ABC/nn.py | 4 +++- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/abc_c6_aux.py b/src/Components/misc/obs_aod/ABC/abc_c6_aux.py index 906924ed..bdb1356d 100644 --- a/src/Components/misc/obs_aod/ABC/abc_c6_aux.py +++ b/src/Components/misc/obs_aod/ABC/abc_c6_aux.py @@ -346,6 +346,8 @@ def make_plots(mxd,expid,ident,I=None): # ------------------------- # mxd.plotKDE(I=I,figfile=expid+"."+ident+"_kde-"+mxd.Target[0][1:]+"-corrected.png") targets = mxd.getTargets(I) + if mxd.nTarget == 1: + targets.shape = targets.shape + (1,) results = mxd.eval(I) # loop through targets for t in range(mxd.nTarget): @@ -595,7 +597,9 @@ def make_error_pdfs(mxd,Input,expid,ident,K=None,I=None,Title=None,netfileRoot=N name = tname if name in mxd.__dict__: if K is None: - targets = mxd.getTargets(I[0])[:,t] + targets = mxd.getTargets(I[0]) + if mxd.nTarget > 1: + targets = targets[:,t] inputs = mxd.getInputs(I[0],Input=Input) knet = mxd.loadnet(netfileRoot+'_Tau.net') @@ -1108,12 +1112,14 @@ def TestStats(mxd,K,C): # regression[*][0:2] = slope, intercept, r-value # out.shape = [ntestobs,nTarget] out, reg = mxd.test(iprint=False) - - mxd.nnr.slope[k,c,:] = reg[:][0] - mxd.nnr.intercept[k,c,:] = reg[:][1] - mxd.nnr.R[k,c,:] = reg[:][2] + reg = np.array(reg) + mxd.nnr.slope[k,c,:] = reg[:,0] + mxd.nnr.intercept[k,c,:] = reg[:,1] + mxd.nnr.R[k,c,:] = reg[:,2] targets = mxd.getTargets(mxd.iTest) + if mxd.nTarget == 1: + targets.shape = targets.shape + (1,) # get other NNR STATS mxd.nnr.rmse[k,c,:] = rmse(out,targets) diff --git a/src/Components/misc/obs_aod/ABC/nn.py b/src/Components/misc/obs_aod/ABC/nn.py index 1961fe38..2fb69d57 100644 --- a/src/Components/misc/obs_aod/ABC/nn.py +++ b/src/Components/misc/obs_aod/ABC/nn.py @@ -238,7 +238,9 @@ def plotScat(self,iTarget=0,bins=None,I=None,figfile=None): """ if I is None: I = self.iTest # Testing data by default results = self.eval(I)[:,iTarget] - targets = self.getTargets(I)[:,iTarget] + targets = self.getTargets(I) + if self.nTarget > 1: + targets = targets[:,iTarget] original = log(self.__dict__['m'+self.Target[iTarget][1:]][I] + 0.01) if bins == None: bins = arange(-5., 1., 0.1 ) From 2ec86d7ee56093064ffad98a884fb318c876b2cc Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 23 Jan 2023 12:43:56 -0500 Subject: [PATCH 27/67] PC - changes to nn, abc_aux to deal with new version of sklearn, which handles the kfolds differently --- src/Components/misc/obs_aod/ABC/abc_c6_aux.py | 9 ++++++--- src/Components/misc/obs_aod/ABC/abc_viirs.py | 6 +++--- src/Components/misc/obs_aod/ABC/nn.py | 2 +- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/abc_c6_aux.py b/src/Components/misc/obs_aod/ABC/abc_c6_aux.py index bdb1356d..0b9f773a 100644 --- a/src/Components/misc/obs_aod/ABC/abc_c6_aux.py +++ b/src/Components/misc/obs_aod/ABC/abc_c6_aux.py @@ -639,7 +639,10 @@ def make_error_pdfs(mxd,Input,expid,ident,K=None,I=None,Title=None,netfileRoot=N for k,iTest in enumerate(I): mxd.iTest = iTest - targets.append(mxd.getTargets(mxd.iTest)[:,t]) + targets_ = mxd.getTargets(mxd.iTest) + if mxd.nTarget > 1: + targets_ = targets_[:,t] + targets.append(targets_) inputs = mxd.getInputs(mxd.iTest,Input=Input) @@ -652,7 +655,7 @@ def make_error_pdfs(mxd,Input,expid,ident,K=None,I=None,Title=None,netfileRoot=N tau = base_tau_t*np.exp(-1.*np.log(wav/mxd.AE_base_wav)*targets[k]) targets[k] = np.log(tau + 0.01) - base_tau_r = np.exp(knet(inputs)[:,mxd.base_wav_i]) - 0.01 + base_tau_r = np.exp(knet(inputs)[:,mxd.AE_base_wav_i]) - 0.01 tau = base_tau_r*np.exp(-1.*np.log(wav/mxd.AE_base_wav)*results[k]) results[k] = np.log(tau + 0.01) @@ -1186,7 +1189,7 @@ def SummaryPDFs(mxdx,mxdx2=None,varnames=['mRef870','mSre470'],doInt=False): I = [Irange] else: I = [] - for iTrain, iTest in mxdx.kf: + for iTrain, iTest in mxdx.kf.split(arange(np.sum(mxdx.iValid))): I.append(iValid[iTest]) I1, I2, I3, I4 = get_Iquartiles(mxdx,I=I) diff --git a/src/Components/misc/obs_aod/ABC/abc_viirs.py b/src/Components/misc/obs_aod/ABC/abc_viirs.py index ff613540..5fd1843f 100644 --- a/src/Components/misc/obs_aod/ABC/abc_viirs.py +++ b/src/Components/misc/obs_aod/ABC/abc_viirs.py @@ -614,7 +614,7 @@ def __init__ (self, fname, (self.mRef2250 > 0) & \ (self.mSre412 > 0.0) & \ (self.mSre488 > 0.0) & \ - (self.mSre670 > 0.0) + (self.mSre670 > 0.0) # Filter by additional variables @@ -675,7 +675,7 @@ def _train(mxd,expid,c): mxd.savenet(outdir+"/"+expid+'_Tau.net') else: k = 1 - for iTrain, iTest in mxd.kf: + for iTrain, iTest in mxd.kf.split(arange(np.sum(mxd.iValid))): I = arange(mxd.nobs) iValid = I[mxd.iValid] mxd.iTrain = iValid[iTrain] @@ -720,7 +720,7 @@ def _test(mxd,expid,c,plotting=True): make_plots(mxd,expid,ident,I=mxd.iTest) else: k = 1 - for iTrain, iTest in mxd.kf: + for iTrain, iTest in mxd.kf.split(arange(np.sum(mxd.iValid))): I = arange(mxd.nobs) iValid = I[mxd.iValid] mxd.iTrain = iValid[iTrain] diff --git a/src/Components/misc/obs_aod/ABC/nn.py b/src/Components/misc/obs_aod/ABC/nn.py index 2fb69d57..070a41be 100644 --- a/src/Components/misc/obs_aod/ABC/nn.py +++ b/src/Components/misc/obs_aod/ABC/nn.py @@ -150,7 +150,7 @@ def kfold(self,K=3): Only data with an iValid Q/C flag is considered. """ n = self.lon.size - self.kf = KFold(np.sum(self.iValid), n_folds=K, shuffle=True, random_state=n) + self.kf = KFold(n_splits=K, shuffle=True, random_state=n) def getInputs(self,I,Input=None): From 26322fd080c07c50c970b7734706b5290bb5903a Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 23 Jan 2023 13:04:17 -0500 Subject: [PATCH 28/67] PC - changes to abc_viirs to calculate brdf from given raa (VIIRS DB doesn't report sensor and solar azimuths, only raa) --- src/Components/misc/obs_aod/ABC/abc_viirs.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Components/misc/obs_aod/ABC/abc_viirs.py b/src/Components/misc/obs_aod/ABC/abc_viirs.py index 5fd1843f..71cd5a5b 100644 --- a/src/Components/misc/obs_aod/ABC/abc_viirs.py +++ b/src/Components/misc/obs_aod/ABC/abc_viirs.py @@ -235,7 +235,8 @@ def setBRDF(self): Kgeo = brdf.__dict__['Rgeo'+ch] self.__dict__[name] = rtlsReflectance(Kiso,Kgeo,Kvol, self.SolarZenith,self.SensorZenith, - self.SolarAzimuth,self.SensorAzimuth) + None,None, + raa=self.RelativeAzimuth) self.giantList.append(name) From 9c4c1d564dde3215bead225084f3bb620ca50f04 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 23 Jan 2023 15:02:14 -0500 Subject: [PATCH 29/67] PC - changes for abc_viirs and giant_viirs to make the definition of species dominated an optional input. default if fraction species > 0.5 --- src/Components/misc/obs_aod/ABC/abc_viirs.py | 7 +++++-- src/Components/misc/obs_aod/ABC/giant_viirs.py | 5 +++-- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/abc_viirs.py b/src/Components/misc/obs_aod/ABC/abc_viirs.py index 71cd5a5b..bf0829fd 100644 --- a/src/Components/misc/obs_aod/ABC/abc_viirs.py +++ b/src/Components/misc/obs_aod/ABC/abc_viirs.py @@ -53,7 +53,8 @@ def setupNN(self,retrieval,expid, 'cloud', 'albedo','fdu','fcc','fsu' ], Input_const = None, Target = ['aTau550',], - K=None): + K=None, + f_balance=0.50): self.retrieval = retrieval @@ -111,8 +112,10 @@ def setupNN(self,retrieval,expid, # Balance the dataset before splitting # No aerosol type should make up more that 35% # of the total number of obs + # f_balance is the fraction that defines whether it's 'dominated' # -------------------------------------- - self.iValid = self.balance(int(self.nobs*0.35)) + self.f_balance = f_balance + self.iValid = self.balance(int(self.nobs*0.35),f_balance=f_balance) # Flatten Input_nnr into one list # ------------------------------- diff --git a/src/Components/misc/obs_aod/ABC/giant_viirs.py b/src/Components/misc/obs_aod/ABC/giant_viirs.py index 1f94ffc8..c443ce90 100644 --- a/src/Components/misc/obs_aod/ABC/giant_viirs.py +++ b/src/Components/misc/obs_aod/ABC/giant_viirs.py @@ -339,7 +339,7 @@ def __init__ (self,filename,xVars=(),only_good=True,tymemax=None): self.nobs = len(self.lon) #-- - def balance(self,N): + def balance(self,N,frac=0.50): """ Return indices of observations so that each species does not have more than N observations. This is meant to be performed with a reduced dataset. @@ -349,7 +349,8 @@ def balance(self,N): for f in (self.fdu,self.fss,self.fcc,self.fsu): - J = f>0.5 # all obs for which species dominate + J = f>frac # all obs for which species dominate + J = J & ~I # only obs that haven't already been selected n = len(self.lon[J]) # no. obs for this species P = random.permutation(n) # randomize obs for this species m = min(n,N) # keep this many From dffef545075dbd5f1cd6ba30aa3d0c675564524e Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 23 Jan 2023 15:47:55 -0500 Subject: [PATCH 30/67] PC - add option to filter by DB land algorithm flag --- src/Components/misc/obs_aod/ABC/abc_viirs.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/Components/misc/obs_aod/ABC/abc_viirs.py b/src/Components/misc/obs_aod/ABC/abc_viirs.py index bf0829fd..936f3bb7 100644 --- a/src/Components/misc/obs_aod/ABC/abc_viirs.py +++ b/src/Components/misc/obs_aod/ABC/abc_viirs.py @@ -570,7 +570,8 @@ def __init__ (self, fname, cloud_thresh=0.70, aFilter=None, NDVI=False, - tymemax=None): + tymemax=None, + algflag=None): """ Initializes the AOD Bias Correction (ABC) for the MODIS Land algorithm. @@ -585,10 +586,17 @@ def __init__ (self, fname, laod --- if True, targets are log-transformed AOD, log(Tau+0.01) tymemax --- truncate the data record in the giant file at tymemax. set to None to read entire data record. + algflag --- DB Land algorithm flag number - this should be a list. Allows for selecting multiple algorithms. + None - don't filter, use all pixels + 0 - hybrid (heterogenous surface) + 1 - vegetated surface + 2 - bright surface + 3 - mixed """ self.verbose = verbose self.laod = laod + self.algflag = algflag DB_LAND.__init__(self,fname,tymemax=tymemax) # initialize superclass @@ -619,6 +627,11 @@ def __init__ (self, fname, (self.mSre412 > 0.0) & \ (self.mSre488 > 0.0) & \ (self.mSre670 > 0.0) + if algflag is not None: + fValid = np.zeros(self.iValid.shape).astype(bool) + for aflag in algflag: + fValid = fValid | (self.algflag == aflag) + self.iValid = self.iValid & fValid # Filter by additional variables From 65e074923aa27d148e3a06bb30a9c1bbb9cc5eb2 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 23 Jan 2023 15:48:51 -0500 Subject: [PATCH 31/67] PC - bug fix. wrong keyword name for balance subroutine --- src/Components/misc/obs_aod/ABC/abc_viirs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Components/misc/obs_aod/ABC/abc_viirs.py b/src/Components/misc/obs_aod/ABC/abc_viirs.py index 936f3bb7..cbccd1f1 100644 --- a/src/Components/misc/obs_aod/ABC/abc_viirs.py +++ b/src/Components/misc/obs_aod/ABC/abc_viirs.py @@ -115,7 +115,7 @@ def setupNN(self,retrieval,expid, # f_balance is the fraction that defines whether it's 'dominated' # -------------------------------------- self.f_balance = f_balance - self.iValid = self.balance(int(self.nobs*0.35),f_balance=f_balance) + self.iValid = self.balance(int(self.nobs*0.35),frac=f_balance) # Flatten Input_nnr into one list # ------------------------------- From f89afc911391dd15a28fcdde83d98008c9b40ff1 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 23 Jan 2023 16:28:13 -0500 Subject: [PATCH 32/67] PC - fix but in giant_virrs that gave the DB Land the wrong ident code --- src/Components/misc/obs_aod/ABC/giant_viirs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Components/misc/obs_aod/ABC/giant_viirs.py b/src/Components/misc/obs_aod/ABC/giant_viirs.py index c443ce90..dc2b935e 100644 --- a/src/Components/misc/obs_aod/ABC/giant_viirs.py +++ b/src/Components/misc/obs_aod/ABC/giant_viirs.py @@ -670,7 +670,7 @@ class DB_LAND(GIANT): def __init__(self,filename,tymemax=None): GIANT.__init__(self,filename,xVars=xDB_LAND,tymemax=tymemax) if self.sat == 'SNPP': - self.ident = 'vsdto' + self.ident = 'vsdbl' self.ident = self.ident + '_' + filename.split('/')[-1].split('.')[0] self.surface = 'land' From 543a82bb01d393f16218b354f74b8ec5b5df392e Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 23 Jan 2023 17:11:35 -0500 Subject: [PATCH 33/67] PC - fix naming of output KDE plots --- src/Components/misc/obs_aod/ABC/abc_viirs.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/abc_viirs.py b/src/Components/misc/obs_aod/ABC/abc_viirs.py index cbccd1f1..f13f2ea2 100644 --- a/src/Components/misc/obs_aod/ABC/abc_viirs.py +++ b/src/Components/misc/obs_aod/ABC/abc_viirs.py @@ -769,9 +769,9 @@ def _test(mxd,expid,c,plotting=True): TestStats(mxd,k-1,c) if plotting: if mxd.angstrom: - make_plots_angstrom(mxd,expid,'.k={}'.format(str(k)),I=mxd.iTest) + make_plots_angstrom(mxd,expid,ident+'.k={}'.format(str(k)),I=mxd.iTest) else: - make_plots(mxd,expid,'.k={}'.format(str(k)),I=mxd.iTest) + make_plots(mxd,expid,ident+'.k={}'.format(str(k)),I=mxd.iTest) k = k + 1 #--------------------------------------------------------------------- From c29a7f13fd2e582d50a2ba625862f3bd92c7f193 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 6 Feb 2023 17:13:08 -0500 Subject: [PATCH 34/67] PC - add option to angstrom interpolate retrieved AOD to AERONET channels in vx04.py --- .../GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 111 ++++++++++++++++-- 1 file changed, 104 insertions(+), 7 deletions(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index 764fdac3..f42e8dd3 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -93,7 +93,7 @@ -# AOD Shannels +# AOD Channels CHANNELS = dict ( DT_LAND = ( 480., 550., 670., 2250.), DT_OCEAN = ( 480., 550., 670., 860., 1240., 1600., 2250. ), @@ -103,8 +103,6 @@ DB_SREF = (412., 488., 670. ), ) - - ALIAS = dict ( Longitude = 'lon', Latitude = 'lat', longitude = 'lon', @@ -186,7 +184,7 @@ class Vx04_L2(object): """ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, - only_good=True,SDS=SDS,alias=None): + only_good=True,SDS=SDS,alias=None,anet_wav=False): """ Reads individual granules or a full day of Level 2 Vx04 files present on a given *Path* and returns a single object with @@ -210,6 +208,9 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, SDS --- Variables to be read from L2 Aerosol files. Must be a dictionary with keys '{Algo}_META' and '{Algo}_{Surface}' ALIAS --- dictionary of alises for SDSs + anet_wav --- angstrom interpolate retrieved AOD to common AERONET wavelengths + e.g. 480 --> 490, 1600 --> 1640, 860 --> 870 + """ @@ -249,7 +250,7 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, Path = [Path, ] self._readList(Path) - #Protect against empty MXD04 files + #Protect against empty VX04 files # -------------------------------- if len(self.Scattering_Angle) == 0: self.nobs = 0 @@ -342,7 +343,7 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, self.nobs = self.Scattering_Angle.shape[0] self.kx = KX[self.sat+'_'+self.algo] self.ident = IDENT[self.sat+'_'+self.algo] - self.channels = CHANNELS[self.algo] + self.channels = np.array(CHANNELS[self.algo]) if Surface == 'LAND': self.sChannels = CHANNELS["{}_SREF".format(Algo)] # LAND surface reflectivity (not the same as algo) @@ -353,7 +354,90 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, elif self.algo == 'DT_OCEAN': self.rChannels = np.array([480.,550.,670.,860.,1240.,1600.,2250.]) - + # Angstrom interpolate to AERONET wavelengths for ease of comparison + # and to faciliate NNR ODS files + # DT: 480 --> 490, 860 --> 870, 1600 --> 1640 + # DB: 488 --> 490, 865 --> 870 + # ------------------------------------------------------------------ + self.anet_wav = anet_wav + if anet_wav: + if 'DT' in self.algo: + tch = 490 + + ch = 550 + ich = np.argmin(np.abs(self.channels - ch)) + aod550 = self.aod[:,ich] + ch = 480 + ich = np.argmin(np.abs(self.channels - ch)) + aod480 = self.aod[:,ich] + + aod_ = self.aodInterpAngs(tch,aod480,aod550,480,550) + self.channels[ich] = tch + self.aod[:.ich] = aod_ + + elif self.algo == 'DT_OCEAN': + tch = 870 + + ch = 670 + ich = np.argmin(np.abs(self.channels - ch)) + aod670 = self.aod[:,ich] + ch = 860 + ich = np.argmin(np.abs(self.channels - ch)) + aod860 = self.aod[:,ich] + + aod_ = self.aodInterpAngs(tch,aod670,aod860,670,860) + self.channels[ich] = tch + self.aod[:.ich] = aod_ + + tch = 1640 + + ch = 2250 + ich = np.argmin(np.abs(self.channels - ch)) + aod2250 = self.aod[:,ich] + ch = 1600 + ich = np.argmin(np.abs(self.channels - ch)) + aod1600 = self.aod[:,ich] + + aod_ = self.aodInterpAngs(tch,aod1600,aod2250,1600,2250) + self.channels[ich] = tch + self.aod[:.ich] = aod_ + + +DB_LAND = ( 412, 488, 550, 670 ), # file has 550 separate +DB_OCEAN = (488., 550., 670., 865., 1240., 1640., 2250.), + + elif 'DB' in self.algo: + tch = 490 + + ch = 550 + ich = np.argmin(np.abs(self.channels - ch)) + aod550 = self.aod[:,ich] + ch = 488 + ich = np.argmin(np.abs(self.channels - ch)) + aod488 = self.aod[:,ich] + + aod_ = self.aodInterpAngs(tch,aod488,aod550,488,550) + self.channels[ich] = tch + self.aod[:.ich] = aod_ + + elif self.algo == 'DB_OCEAN': + tch = 870 + + ch = 670 + ich = np.argmin(np.abs(self.channels - ch)) + aod670 = self.aod[:,ich] + ch = 865 + ich = np.argmin(np.abs(self.channels - ch)) + aod865 = self.aod[:,ich] + + aod_ = self.aodInterpAngs(tch,aod670,aod865,670,865) + self.channels[ich] = tch + self.aod[:.ich] = aod_ + + + + # ODS time variables + # ------------------- if syn_time == None: self.syn_time = None self.time = None @@ -384,6 +468,19 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, if Algo == 'DB': self.cloud = 1. - self.npixels_used.astype(float)/self.npixels_valid.astype(float) + +#--- + def aodInterpAngs(lambda_,tau1,tau2,lambda1,lambda2): + """ + Angstrom-interpolated AOD. + lambda_ = target wavelength for AOD output + """ + I = (tau1>0) & (tau2>0) + angstrom = -np.log(tau1[I]/tau2[I])/np.log(lambda1/lambda2) + tau = MISSING * np.ones(len(tau1)) + tau[I] = tau2[I] * (lambda2/lambda_)**angstrom + return tau + #--- def _readList(self,List): """ From 84b498311fca6eef2e6fc5c9dfe37bb5c976568c Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Wed, 8 Feb 2023 12:12:17 -0500 Subject: [PATCH 35/67] PC - took out old comments --- src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index f42e8dd3..f1d0c097 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -402,10 +402,6 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, self.channels[ich] = tch self.aod[:.ich] = aod_ - -DB_LAND = ( 412, 488, 550, 670 ), # file has 550 separate -DB_OCEAN = (488., 550., 670., 865., 1240., 1640., 2250.), - elif 'DB' in self.algo: tch = 490 From c207a3ad1308be22f3094eb803d5ea102375f633 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Wed, 8 Feb 2023 21:14:34 -0500 Subject: [PATCH 36/67] PC - initial commit of viirs_l2a.py file --- src/Applications/GAAS_App/viirs_l2a.py | 67 ++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) create mode 100755 src/Applications/GAAS_App/viirs_l2a.py diff --git a/src/Applications/GAAS_App/viirs_l2a.py b/src/Applications/GAAS_App/viirs_l2a.py new file mode 100755 index 00000000..a86974c5 --- /dev/null +++ b/src/Applications/GAAS_App/viirs_l2a.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python +# -W ignore::DeprecationWarning + +""" + + Simple wrapper script to parse Prep config file and create ODS with VIIRS NNR aerosol retrievals. + + February 2011. + arlindo.dasilva@nasa.gov + + Based on modis_l2a.py + 2023 + patricia.castellanos@nasa.gov +""" + +from os import system +from optparse import OptionParser +from MAPL import Config + +if __name__ == "__main__": + +# Parse command line options +# -------------------------- + parser = OptionParser(usage="Usage: %prog prep_config_file isotime", + version='viirs_l2a-1.0.0' ) + parser.add_option("-n", "--dryrun", + action="store_true", dest="dryrun", + help="Dry run.") + + (options, args) = parser.parse_args() + + if len(args) == 2: + prep_config, isotime = args + else: + parser.error("must have 2 arguments: prep_config_filename isotime") + + # Parse prep config + # ----------------- + cf = Config(prep_config,delim=' = ') + + Options = " --expid=" + cf('VIIRS_L2A_EXPID') + \ + " --l2_dir=" + cf('VIIRS_L2A_L2_DIR') + \ + " --res=" + cf('VIIRS_L2A_RESOLUTION') + \ + " --dir=" + cf('VIIRS_L2A_OUT_DIR') + \ + " --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') + + if cf('VIIRS_L2A_OVERWRITE').upper() == 'YES': Options += " --force" + if cf('VIIRS_L2A_VERBOSE').upper() == 'YES': Options += " -v" + + # Generate products + # ----------------- + i = 0 + Coll = cf('VIIRS_L2A_COLLECTION').split(',') + for ident in cf('VIIRS_L2A_IDENTS').split(','): + coll = Coll[i] + cmd = "vx04_l2a.py %s --collection=%s %s %s "%(Options,Coll[i],ident,isotime) + print cmd + if not options.dryrun: + if system(cmd): + raise ValueError, "vx04_l2a.py failed for %s on %s "%(ident,isotime) + + i += 1 + + From b7dfb0616ea64dca817615ac68f8153b15334a93 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Wed, 8 Feb 2023 21:14:54 -0500 Subject: [PATCH 37/67] PC - initial commit of vx04l2a.py file --- src/Applications/GAAS_App/vx04_l2a.py | 241 ++++++++++++++++++++++++++ 1 file changed, 241 insertions(+) create mode 100755 src/Applications/GAAS_App/vx04_l2a.py diff --git a/src/Applications/GAAS_App/vx04_l2a.py b/src/Applications/GAAS_App/vx04_l2a.py new file mode 100755 index 00000000..134bd737 --- /dev/null +++ b/src/Applications/GAAS_App/vx04_l2a.py @@ -0,0 +1,241 @@ +#!/usr/bin/env python +# -W ignore::DeprecationWarning + +""" + A Python script to create NNR retrievals. + It uses class VX04 to directly read VIIRS Aerosol Level 2 + Files (AERDB and AERDT). + + This utility reads VIIRS Level2 files and creates an ODS file with + NNR retrievals, as well as a *gritas* type gridded output. + + 2023 Based on mxd04_l2a.py + patricia.castellanos@nasa.gov +""" + +import warnings +warnings.simplefilter('ignore',DeprecationWarning) +warnings.simplefilter('always',UserWarning) + +import os +import sys +import subprocess + +from time import clock +from optparse import OptionParser # Command-line args +from dateutil.parser import parse as isoparse +from vx04_nnr import Vx04_NNR +from MAPL import strTemplate + +Ident = dict( vsnppdto = ('SNPP04','dt_ocean'), + vsnppdtl = ('SNPP04','dt_land'), + vsnppdbo = ('SNPP04','db_ocean'), + vsnppdbl = ('SNPP04','db_land') + ) + +#--------------------------------------------------------------------- +def makethis_dir(filename): + """Creates the relevant directory if necessary.""" + path, filen = os.path.split(filename) + if path != '': + rc = os.system('mkdir -p '+path) + if rc: + raise IOError, "could not create directory "+path + +#--------------------------------------------------------------------- + +if __name__ == "__main__": + + expid = 'nnr_001' + ident = 'vsnppdbl' + +# Defaults may be platform dependent +# ---------------------------------- + if os.path.exists('/nobackup/VIIRS/Level2/'): # New calculon + l2_path = '/nobackup/VIIRS/' + out_dir = '/nobackup/NNR/VIIRS/%coll/Level%lev/%prod/Y%y4/M%m2' + nn_file = '/nobackup/NNR/Net/VIIRS/nnr_003.%ident_Tau.net' + blank_ods = '/nobackup/NNR/Misc/blank.ods' + aer_x = '/nobackup/NNR/Misc/tavg1_2d_aer_Nx' + else: # Must be somewhere else, no good defaults + out_dir = './' + l2_path = './' + nn_file = '%ident_Tau.net' + blank_ods = 'blank.ods' + aer_x = 'tavg1_2d_aer_Nx' + + out_tmpl = '%s.%prod_L%leva.%algo.%y4%m2%d2_%h2%n2z.%ext' + coll = '002' + res = 'c' + +# Parse command line options +# -------------------------- + parser = OptionParser(usage="Usage: %prog [options] ident isotime", + version='vx04_l2a-1.0.0' ) + + + parser.add_option("-x", "--expid", dest="expid", default=expid, + help="Experiment id (default=%s)"\ + %expid ) + + parser.add_option("-d", "--dir", dest="out_dir", default=out_dir, + help="Output directory (default=%s)"\ + %out_dir ) + + parser.add_option("-A", "--aer_x", dest="aer_x", default=aer_x, + help="GrADS ctl for speciated AOD file (default=%s)"\ + %aer_x ) + + parser.add_option("-B", "--blank_ods", dest="blank_ods", default=blank_ods, + help="Blank ODS file name for fillers (default=%s)"\ + %blank_ods ) + + parser.add_option("-C", "--collection", dest="coll", default=coll, + help="VIIRS collection (default=%s)"\ + %coll ) + + parser.add_option("-o", "--fname", dest="out_tmpl", default=out_tmpl, + help="Output file name template (default=%s); ODS file name will be derived from it by changing extension to '.ods' and replacing 'Level3' with 'Level2'."\ + %out_tmpl ) + + parser.add_option("-L", "--l2_dir", dest="l2_path", default=l2_path, + help="Top directory for VIIRS Level 2 files (default=%s)"\ + %l2_path ) + + parser.add_option("-N", "--net", dest="nn_file", default=nn_file, + help="Neural net file template (default=%s)"\ + %nn_file ) + + parser.add_option("-r", "--res", dest="res", default=res, + help="Resolution for gridded output (default=%s)"\ + %out_tmpl ) + + 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)") + + parser.add_option("-F", "--force", + action="store_true", dest="force",default=False, + help="Overwrites output file") + + parser.add_option("-v", "--verbose", + action="store_true", dest="verbose",default=False, + help="Turn on verbosity.") + + parser.add_option("--writenpz", dest="writenpz", default=False, + help="Write an ungridded npz file in addition to ODS and gridded files (default=False)") + + (options, args) = parser.parse_args() + + if len(args) == 2: + ident, isotime = args + prod, algo = Ident[ident] + else: + parser.error("must have 3 arguments: ident, date and time") + +# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + if options.verbose: + print "" + print " VIIRS Level 2A Processing" + print " -------------------------" + print "" + t0 = clock() + +# Time variables +# -------------- + syn_time = isoparse(isotime) + nymd = str(syn_time.date()).replace('-','') + nhms = str(syn_time.time()).replace(':','') + +# Form output gridded file name +# ----------------------------- + out_tmpl = options.out_dir+'/'+options.out_tmpl + out_tmpl = out_tmpl.replace('%coll',options.coll).replace('%prod',prod).replace('%algo',algo).replace('%lev','3').replace('%ext','nc4') + out_file = strTemplate(out_tmpl,expid=options.expid,nymd=nymd,nhms=nhms) + name, ext = os.path.splitext(out_file) + if os.path.exists(out_file) and (options.force is not True): + print "vx04_l2a: Output Gridded file <%s> exists --- cannot proceed."%out_file + raise IOError, "Specify --force to overwrite existing output file." + if os.path.exists(out_file) and options.force: + os.remove(out_file) + +# Form ODS file name +# ------------------ + ods_tmpl = options.out_dir+'/'+options.out_tmpl + ods_tmpl = ods_tmpl.replace('%coll',options.coll).replace('%prod',prod).replace('%algo',algo).replace('%lev','2').replace('%ext','ods') + ods_file = strTemplate(ods_tmpl,expid=options.expid,nymd=nymd,nhms=nhms) + if os.path.exists(ods_file) and (options.force is not True): + print "vxd04_l2a: Output ODS file <%s> exists --- cannot proceed."%ods_file + raise IOError, "Specify --force to overwrite existing output file." + if os.path.exists(ods_file) and options.force: + os.remove(ods_file) + +# Aerosol composition file name +# ----------------------------- + if options.aer_x[-3:] == 'nc4': + aer_x = strTemplate(options.aer_x,expid=options.expid,nymd=nymd,nhms=nhms) + else: + aer_x = options.aer_x + + +# VIIRS Level 2 NNR Aerosol Retrievals +# ------------------------------------ + if options.verbose: + print "NNR Retrieving %s %s on "%(prod,algo.upper()),syn_time + + viirs = Vx04_NNR(options.l2_path,prod,algo.upper(),syn_time,aer_x, + coll=options.coll, + cloud_thresh=0.7, + cloudFree = 0.0, + aodmax = 1.0, + verbose=options.verbose) + if viirs.nobs < 1: + if options.verbose: + print 'WARNING: no observation for this time in file <%s>'%ods_file + + elif any(viirs.iGood) == False: + if options.verbose: + print 'WARNING: no GOOD observation for this time in file <%s>'%ods_file + viirs.nobs = 0 + + nn_file = options.nn_file.replace('%ident',ident) + viirs.apply(nn_file) + +# Write ODS +# --------- + makethis_dir(ods_file) + if viirs.nobs>0: + viirs.writeODS(ods_file,revised=True) + else: + if os.system('ods_blank.x %s %s %s %s'%(options.blank_ods,nymd,nhms,ods_file)): + warnings.warn('cannot create empty output file <%s>'%ods_file) + else: + if options.verbose: + print "[w] Wrote empty ODS file "+ods_file + +# Write gridded output file (revised channels only) +# ------------------------------------------------- + makethis_dir(out_file) + if viirs.nobs>0: + if str.isdigit(options.res): + viirs.writeg(filename=out_file,refine=int(options.res),channels=viirs.channels_) + else: + viirs.writeg(filename=out_file,res=options.res,channels=viirs.channels_) + +# Write ungridded data +# -------------------- + if options.writenpz: + name, ext = os.path.splitext(out_file) + npz_file = name.replace('Level3','Level2') + '.npz' + makethis_dir(npz_file) + viirs.write(npz_file) + +# Compress nc output unless the user disabled it +# ---------------------------------------------- + if viirs.nobs>0: + if not options.uncompressed: + if subprocess.call("n4zip " + out_file,shell=True): + warnings.warn('cannot compress output file <%s>'%out_file) + if subprocess.call("n4zip " + ods_file,shell=True): + warnings.warn('cannot compress output ods file <%s>'%ods_file) From 7f77fcaf0ee541bd29fe7d3d382f2bac1cccf896 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Wed, 8 Feb 2023 21:16:46 -0500 Subject: [PATCH 38/67] PC - add viirs drivers to GAAS Cmakelists --- src/Applications/GAAS_App/CMakeLists.txt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/Applications/GAAS_App/CMakeLists.txt b/src/Applications/GAAS_App/CMakeLists.txt index 3932fe11..57a9f914 100644 --- a/src/Applications/GAAS_App/CMakeLists.txt +++ b/src/Applications/GAAS_App/CMakeLists.txt @@ -40,6 +40,8 @@ set (PYSCRIPTS modis_l2a.py mxd04_l2a.py patmosx_l2a.py + viirs_l2a.py + vx04_l2a.py ) set (SCRIPTS @@ -76,6 +78,7 @@ set (RCFILES set (PCFFILES avhrr_l2a.pcf modis_l2a.pcf + viirs_l2a.pcf ) install ( From a35048294909f203b0bd9dce978b057e38527a88 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Wed, 8 Feb 2023 21:17:47 -0500 Subject: [PATCH 39/67] PC - initial changes to vx04_nnr to make it work for viirs --- src/Applications/GAAS_App/vx04_nnr.py | 338 ++++++++++---------------- 1 file changed, 127 insertions(+), 211 deletions(-) diff --git a/src/Applications/GAAS_App/vx04_nnr.py b/src/Applications/GAAS_App/vx04_nnr.py index e980916a..fc3c6905 100644 --- a/src/Applications/GAAS_App/vx04_nnr.py +++ b/src/Applications/GAAS_App/vx04_nnr.py @@ -1,108 +1,55 @@ """ -This module implements the VIIRS NNR AOD retrievals. +This module implements the MODIS NNR AOD retrievals. -This version works from VIIRS DT & DT Level 2 files. +This version works from MODIS MOD04/MYD04 Level 2 files. + +Modified to work for VIIRS. Feb 2023 P. Castellanos """ import os, sys -import warnings -from pyobs.vx04 import Vx04_L2, MISSING, granules, BEST +from pyobs.vx04 import Vx04_L2, MISSING, granules, SDS from ffnet import loadnet -from numpy import c_ as cat -from numpy import copy, ones, sin, cos, exp, arccos, pi, any, log import numpy as np -from pyobs.bits import BITS - -# SDS to be read in -# ------------ -SDS = dict ( - DT_META = ('longitude', 'latitude', - 'sensor_zenith_angle', 'sensor_azimuth_angle', - 'solar_zenith_angle', 'solar_azimuth_angle', - 'Scattering_Angle', 'Glint_Angle'), - DT_LAND = ( 'Mean_Reflectance_Land', - 'Corrected_Optical_Depth_Land', - 'Surface_Reflectance_Land', - 'Aerosol_Cloud_Fraction_Land', - 'Land_Ocean_Quality_Flag'), - DT_OCEAN = ( 'Effective_Optical_Depth_Average_Ocean', - 'Aerosol_Cloud_Fraction_Ocean', - 'Mean_Reflectance_Ocean', - 'Land_Ocean_Quality_Flag' ), - DB_META = ('Longitude', 'Latitude', 'Scan_Start_Time', - 'Viewing_Zenith_Angle', 'Relative_Azimuth_Angle', - 'Solar_Zenith_Angle', - 'Scattering_Angle' ), - DB_LAND = ('Aerosol_Optical_Thickness_550_Land', - 'Spectral_Aerosol_Optical_Thickness_Land', - 'Spectral_Surface_Reflectance', - 'Spectral_TOA_Reflectance_Land', - 'Aerosol_Optical_Thickness_QA_Flag_Land', - 'Algorithm_Flag_Land', - 'Number_Of_Pixels_Used_Land', - 'Number_Valid_Pixels'), - DB_OCEAN = ('Aerosol_Optical_Thickness_550_Ocean', - 'Spectral_Aerosol_Optical_Thickness_Ocean', - 'Spectral_TOA_Reflectance_Ocean', - 'Aerosol_Optical_Thickness_QA_Flag_Ocean', - 'Algorithm_Flag_Ocean', - 'Number_Of_Pixels_Used_Ocean', - 'Number_Valid_Pixels') - ) - - -# Channels for TOA reflectance -# ----------------------------- -CHANNELS = dict ( - DT_LAND = ( 480, 670, 2250), - DT_OCEAN = ( 480, 550, 670, 860, 1240, 1600, 2250 ), - DB_LAND = ( 412, 488, 550, 670, 865, 1240, 1640, 2250 ), - DB_OCEAN = ( 412, 488, 550, 670, 865, 1240, 1640, 2250 ), - ) - -SCHANNELS = dict ( - DT_LAND = ( 480, 670, 2250 ), - DB_LAND = ( 412, 488, 670 ), - ) # Translate Inputs between NNR and MODIS classes # ----------------------------------------------- -TranslateInput = dict ( mRef412 = ('reflectance',412), - mRef480 = ('reflectance',480), - mRef550 = ('reflectance',550), - mRef670 = ('reflectance',670), - mRef860 = ('reflectance',860), - mRef1240 = ('reflectance',1240), - mRef1600 = ('reflectance',1600), - mRef2250 = ('reflectance',2250), - mSre412 = ('sfc_reflectance',412), - mSre480 = ('sfc_reflectance',480), - mSre670 = ('sfc_reflectance',670), - mSre2250 = ('sfc_reflectance',2250), - ) - -for var in ( 'ScatteringAngle','GlintAngle', - 'SolarAzimuth', 'SolarZenith', - 'SensorAzimuth','SensorZenith', - 'cloud','qa_flag' ): - TranslateInput[var] = (var,) +def TranslateInput(key): + if 'mRef' in key: + prefix = 'reflectance' + channel = int(key[4:]) + output = (prefix,channel) + elif 'mSre' in key: + prefix = 'sfc_reflectance' + channel = int(key[4:]) + output = (prefix,channel) + else: + output = (key,) + + + return output + # Translate Targets between ANET and MODIS classes # ------------------------------------------------ -TranslateTarget = dict ( aTau440 = ( 'aod_', 440 ), - aTau470 = ( 'aod_', 470 ), - aTau500 = ( 'aod_', 500 ), - aTau550 = ( 'aod_', 550 ), - aTau660 = ( 'aod_', 660 ), - aTau870 = ( 'aod_', 870 ), - ) +def TranslateTarget(key): + if 'aTau' in key: + prefix = 'aod_' + channel = int(key[4:]) + output = (prefix, channel) + elif 'aAE' in key: + prefix = 'aod_' + channel = int(key[3:]) + output = (prefix,channel) + + return output + class Vx04_NNR(Vx04_L2): """ This class extends VIIRS by adding methods for producing NNR AOD retrievals based on the Neural Net model developed - with class *abc_c6*. + with class *abc_viirs*. """ def __init__(self,l2_path,sat,algo,syn_time,aer_x, @@ -111,30 +58,36 @@ def __init__(self,l2_path,sat,algo,syn_time,aer_x, scat_thresh=170.0, cloudFree=None, aodmax=1.0, - coll='011',verbose=0): + coll='002', + nsyn=8, + verbose=0): """ Contructs a VX04 object from VIIRS Aerosol Level 2 granules. On input, - l2_path --- top directory for the VIIRS Level 2 files; - it must have subdirectories AERDB_inst and AERDT_inst. - sat --- either *SNPP* (Suomi-NPP) or *NOAA20* (NOAA-20) - algo --- aerosol algorithm: DT_LAND, DT_OCEAN, DB_LAND, or DB_OCEAN + l2_path --- top directory for the VIIRS Level 2 files; + it must have subdirectories AERDB and/or AERDT. + sat --- satellite, either SNPP or NOAAXX (e.g. NOAA20) + algo --- Algorithm: DT_LAND, DT_OCEAN, DB_LAND or DB_OCEAN syn_time --- synoptic time + aer_x --- GEOS control file of aerosol collection - cloud_tresh --- cloud fraction treshhold - cloudFree --- cloud fraction threshhold for assuring no cloud contaminations when aod is > aodmax + Optional parameters: + glint_thresh --- glint angle threshhold + scat_thresh --- scattering angle thresshold + 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 + coll --- VIIRS data collection + nsyn --- number of synoptic times The following attributes are also defined: fractions dust, sea salt, BC+OC, sulfate + aod_coarse + wind It also performs Q/C, setting attribute iGood. On, input, *cloud_thresh* is the cloud fraction limit. - - (For right now this is not implemented) - When DEEP BLUE algorithm is requested, filters for - retrievals where DARK TARGET obs are unavailable. """ self.verbose = verbose @@ -143,99 +96,42 @@ def __init__(self,l2_path,sat,algo,syn_time,aer_x, self.aodmax = aodmax # Initialize superclass - # --------------------- - Files = granules(l2_path,algo,sat,syn_time,coll=coll) - Vx04_L2.__init__(self,Files,algo,syn_time, + # set anet_wav to True so MODIS wavelengths align with AERONET + # Needed for ODS files + # ------------------------------------------------------------- + Files = granules(l2_path,algo,sat,syn_time,coll=coll,nsyn=nsyn) + Vx04_L2.__init__(self,Files,algo,syn_time=syn_time,nsyn=nsyn, only_good=True, - SDS=SDS, - Verb=verbose) + SDS=SDS, + alias=None, + Verb=verbose, + anet_wav=True) if self.nobs < 1: return # no obs, nothing to do - # Channels - # ----------------------------- - self.rChannels = CHANNELS[algo] - if algo in SCHANNELS: - self.sChannels = SCHANNELS[algo] - - - # DB Algorithm only used when Dark Target data is unavailable - # (Not currently implemented) - # -------------------------------------------------------------- - -# if algo == "DEEP": -# # Get DARK TARGET qa_flag -# self.qa_flag_lnd = BITS(self.Quality_Assurance_Land[:,0])[1:4] -# lndGood = self.qa_flag_lnd == BEST -# lndGood = lndGood & (self.cloud_lnd < cloud_thresh) -# rChannels = CHANNELS["LAND"] -# sChannels = SCHANNELS["LAND"] -# for i,c in enumerate(rChannels): -# lndGood = lndGood & (self.reflectance_lnd[:,i]>0) - -# for i,c in enumerate(sChannels): -# lndGood = lndGood & (self.sfc_reflectance_lnd[:,i]>0) - -# self.iGood = (self.qa_flag == BEST) & ~lndGood - -# # Keep only "good" observations -# # ----------------------------- -# m = self.iGood -# for sds in self.SDS: -# rank = len(self.__dict__[sds].shape) -# if rank == 1: -# self.__dict__[sds] = self.__dict__[sds][m] -# elif rank == 2: -# self.__dict__[sds] = self.__dict__[sds][m,:] -# else: -# raise IndexError, 'invalid rank=%d'%rank - -# # Reset aliases -# for sds in self.SDS: -# if sds in self.ALIAS: -# self.__dict__[self.ALIAS[sds]] = self.__dict__[sds] - - -# self.qa_flag = self.qa_flag[m] -# self.aod = self.aod[m,:] -# self.time = self.time[m] -# self.Time = self.Time[m] -# self.iGood = self.iGood[m] -# self.nobs = self.Longitude.shape[0] - -# if self.nobs < 1: -# return # no obs, nothing to do - - # Q/C # --- self.iGood = self.cloud<cloud_thresh -# if algo == "LAND": -# self.iGood = self.iGood & (self.cloud_deep<cloud_thresh) -# elif algo == "DEEP": -# self.iGood = self.iGood & (self.cloud_lnd<cloud_thresh) for i,c in enumerate(self.rChannels): self.iGood = self.iGood & (self.reflectance[:,i]>0) - if algo in SCHANNELS: + if "LAND" in algo: + self.iGood = self.iGood & (self.ScatteringAngle < scat_thresh) for i,c in enumerate(self.sChannels): self.iGood = self.iGood & (self.sfc_reflectance[:,i]>0) if "OCEAN" in algo: self.iGood = self.iGood & (self.GlintAngle > glint_thresh) - if "LAND" in algo: - self.iGood = self.iGood & (self.ScatteringAngle < scat_thresh) - - if any(self.iGood) == False: + if np.any(self.iGood) == False: print "WARNING: Strange, no good obs left to work with" return # Create attribute for holding NNR predicted AOD # ---------------------------------------------- - self.aod_ = MISSING * ones((self.nobs,len(self.channels))) + self.aod_ = MISSING * ones((self.nobs,len(self.channels))) # Make sure same good AOD is kept for gridding # -------------------------------------------- @@ -246,12 +142,12 @@ def __init__(self,l2_path,sat,algo,syn_time,aer_x, # Angle transforms: for NN calculations we work with cosine of angles # ------------------------------------------------------------------- - self.ScatteringAngle = cos(self.ScatteringAngle*pi/180.0) - self.SensorAzimuth = cos(self.SensorAzimuth*pi/180.0) - self.SensorZenith = cos(self.SensorZenith*pi/180.0) - self.SolarAzimuth = cos(self.SolarAzimuth*pi/180.0) - self.SolarZenith = cos(self.SolarZenith*pi/180.0) - self.GlintAngle = cos(self.GlintAngle*pi/180.0) + self.ScatteringAngle = np.cos(self.ScatteringAngle*np.pi/180.0) + self.SensorAzimuth = np.cos(self.SensorAzimuth*np.pi/180.0) + self.SensorZenith = np.cos(self.SensorZenith*np.pi/180.0) + self.SolarAzimuth = np.cos(self.SolarAzimuth*np.pi/180.0) + self.SolarZenith = np.cos(self.SolarZenith*np.pi/180.0) + self.GlintAngle = np.cos(self.GlintAngle*np.pi/180.0) # Get fractional composition # ------------------------------ @@ -289,6 +185,14 @@ def speciate(self,aer_x,Verbose=False): except: pass # ignore it for systems without nitrates + # Handle brown carbon + # -------------------- + try: + self.sampleFile(aer_x,onlyVars=('BRCEXTTAU',),Verbose=Verbose) + self.fcc += self.sample.BRCEXTTAU / s.TOTEXTTAU + except: + pass # ignore it for systems without brown carbon + del self.sample #--- @@ -360,22 +264,16 @@ def _getInputs(self): # ---------------- first = True for inputName in self.net.InputNames: - try: - iName = TranslateInput[inputName] - except: - iName = inputName + iName = TranslateInput(inputName) if self.verbose>0: print 'Getting NN input ',iName # Retrieve input # -------------- - if type(iName) is str: - input = self.__dict__[iName][:] - - elif len(iName) == 2: + if len(iName) == 2: name, ch = iName - if 'mSre' in inputName: # LAND, surface reflectivity + if 'mSre' in inputName: # LAND surface reflectivity k = list(self.sChannels).index(ch) # index of channel elif 'mRef' in inputName: # TOA reflectances k = list(self.rChannels).index(ch) # index of channel @@ -395,7 +293,7 @@ def _getInputs(self): inputs = input first = False else: - inputs = cat[inputs,input] + inputs = np.c_[inputs,input] # Keep only good observations # --------------------------- @@ -407,11 +305,11 @@ def apply(self,nnFile): Apply bias correction to AOD. """ - # Stop here is no good obs available + # Stop here if no good obs available # ---------------------------------- if self.nobs == 0: return # no data to work with - if any(self.iGood) == False: + if np.any(self.iGood) == False: return # no good data to work with # Load the Neural Net @@ -422,24 +320,50 @@ def apply(self,nnFile): # --------------------- targets = self.net(self._getInputs()) - - # Targets do not have to be in standard retrieval + # If target is angstrom exponent + # calculate AOD + # ------------------------------ + doAE = False + for targetName in self.net.TargetNames: + if 'AE' in targetName: + doAE = True + + if doAE: + for i,targetName in enumerate(self.net.TargetNames): + if 'Tau' in targetName: + name, base_wav = TranslateTarget[targetName] + base_wav = np.float(base_wav) + base_tau = targets[:,i] + if self.net.laod: + base_tau = exp(base_tau) - 0.01 # inverse + for i,targetName in enumerate(self.net.TargetNames): + if 'AE' in targetName: + AE = targets[:,i] + name, wav = TranslateTarget[targetName] + wav = np.float(wav) + data = base_tau*np.exp(-1.*AE*np.log(wav/base_wav)) + if self.net.laod: + targets[:,i] = np.log(data + 0.01) + else: + targets[:,i] = data + + # Targets do not have to be in VIIRS retrieval # ---------------------------------------------- for i,targetName in enumerate(self.net.TargetNames): - name, ch = TranslateTarget[targetName] + name, ch = TranslateTarget(targetName) try: k = list(self.channels).index(ch) # index of channel except: # add new target channel to end self.channels += (ch,) - self.aod = np.append(self.aod,MISSING*ones((self.nobs,1)),axis=1) - self.aod_ = np.append(self.aod_,MISSING*ones((self.nobs,1)),axis=1) + self.aod = np.append(self.aod,MISSING*np.ones((self.nobs,1)),axis=1) + self.aod_ = np.append(self.aod_,MISSING*np.ones((self.nobs,1)),axis=1) # Replace targets with unbiased values # ------------------------------------ self.channels_ = [] # channels being revised for i,targetName in enumerate(self.net.TargetNames): - name, ch = TranslateTarget[targetName] + name, ch = TranslateTarget(targetName) if self.verbose>0: if self.net.laod: print "NN Retrieving log(AOD+0.01) at %dnm "%ch @@ -448,7 +372,7 @@ def apply(self,nnFile): k = list(self.channels).index(ch) # index of channel self.channels_ = self.channels_ + [ch,] if self.net.laod: - result = exp(targets[:,i]) - 0.01 # inverse + result = np.exp(targets[:,i]) - 0.01 # inverse else: result = targets[:,i] @@ -458,22 +382,16 @@ def apply(self,nnFile): # Do extra cloud filtering if required if self.cloudFree is not None: cloudy = (self.cloud>=self.cloudFree) -# 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) contaminated = np.zeros(np.sum(self.iGood)).astype(bool) for targetName in self.net.TargetNames: - name, ch = TranslateTarget[targetName] + 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] ) for targetName in self.net.TargetNames: - name, ch = TranslateTarget[targetName] + name, ch = TranslateTarget(targetName) k = list(self.channels).index(ch) # index of channel self.__dict__[name][self.iGood,k][contaminated] = MISSING @@ -494,25 +412,23 @@ def apply(self,nnFile): from datetime import datetime - l2_path = '/nobackup/VIIRS/' - algo = 'DB_LAND' + l2_path = '/nobackup/VIIRS/AERDB/SNPP' sat = 'SNPP' - coll = '011' + algo = 'DB_LAND' + coll = '002' aer_x = '/nobackup/NNR/Misc/tavg1_2d_aer_Nx' - syn_time = datetime(2016,12,19,15,0,0) + syn_time = datetime(2012,03,01,00,0,0) - if algo == 'OCEAN': - nn_file = '/nobackup/NNR/Net/nnr_003.mydo_Tau.net' - elif algo == 'LAND': - nn_file = '/nobackup/NNR/Net/nnr_003.mydl_Tau.net' - elif algo == 'DEEP': - nn_file = '/nobackup/NNR/Net/nnr_003.mydd_Tau.net' + if algo == 'DB_LAND': + nn_file = '/nobackup/NNR/Net/VIIRS/nnr_001.vsnppdbl_Tau.net' - m = Vx04_NNR(l2_path,algo.upper(),sat,syn_time,aer_x, + m = Vx04_NNR(l2_path,sat,algo,syn_time,aer_x, coll=coll, cloud_thresh=0.7, verbose=True) m.apply(nn_file) aod = m.aod_ + + From 3c445ce9e017027b5d6a8dca9245af9835216770 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 9 Feb 2023 12:55:34 -0500 Subject: [PATCH 40/67] PC - example pcf files for processing viirs ods files --- src/Applications/GAAS_App/viirs_l2a.pcf | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) create mode 100644 src/Applications/GAAS_App/viirs_l2a.pcf diff --git a/src/Applications/GAAS_App/viirs_l2a.pcf b/src/Applications/GAAS_App/viirs_l2a.pcf new file mode 100644 index 00000000..c6c889ce --- /dev/null +++ b/src/Applications/GAAS_App/viirs_l2a.pcf @@ -0,0 +1,22 @@ +#VIIRS_L2A processing + +VIIRS_L2A_VERBOSE = YES + +VIIRS_L2A_EXPID = nnr_001 +VIIRS_L2A_IDENTS = vsnppdbl +VIIRS_L2A_COLLECTION = 002 + +VIIRS_L2A_L2_DIR = /nobackup/VIIRS +VIIRS_L2A_OUT_DIR = /nobackup/NNR/VIIRS/%coll_nnr_001/Level%lev/%prod/Y%y4/M%m2 + +VIIRS_L2A_OVERWRITE = YES +VIIRS_L2A_OUT_TEMPLATE = '%s.%prod_L%leva.%algo.%y4%m2%d2_%h2%n2z.%ext' +VIIRS_L2A_RESOLUTION = 16 + +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 + +#END VIIRS_L2A processing From 4ed3cda08f40bfca829b7b2940103d3fa2714f95 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 9 Feb 2023 12:56:16 -0500 Subject: [PATCH 41/67] PC - fix call to find viirs granules. need satellite name not product --- src/Applications/GAAS_App/vx04_l2a.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/Applications/GAAS_App/vx04_l2a.py b/src/Applications/GAAS_App/vx04_l2a.py index 134bd737..7a2e3845 100755 --- a/src/Applications/GAAS_App/vx04_l2a.py +++ b/src/Applications/GAAS_App/vx04_l2a.py @@ -27,10 +27,10 @@ from vx04_nnr import Vx04_NNR from MAPL import strTemplate -Ident = dict( vsnppdto = ('SNPP04','dt_ocean'), - vsnppdtl = ('SNPP04','dt_land'), - vsnppdbo = ('SNPP04','db_ocean'), - vsnppdbl = ('SNPP04','db_land') +Ident = dict( vsnppdto = ('SNPP','dt_ocean'), + vsnppdtl = ('SNPP','dt_land'), + vsnppdbo = ('SNPP','db_ocean'), + vsnppdbl = ('SNPP','db_land') ) #--------------------------------------------------------------------- @@ -129,7 +129,8 @@ def makethis_dir(filename): if len(args) == 2: ident, isotime = args - prod, algo = Ident[ident] + sat, algo = Ident[ident] + prod = sat + '04' else: parser.error("must have 3 arguments: ident, date and time") @@ -182,9 +183,9 @@ def makethis_dir(filename): # VIIRS Level 2 NNR Aerosol Retrievals # ------------------------------------ if options.verbose: - print "NNR Retrieving %s %s on "%(prod,algo.upper()),syn_time + print "NNR Retrieving %s %s on "%(sat,algo.upper()),syn_time - viirs = Vx04_NNR(options.l2_path,prod,algo.upper(),syn_time,aer_x, + viirs = Vx04_NNR(options.l2_path,sat,algo.upper(),syn_time,aer_x, coll=options.coll, cloud_thresh=0.7, cloudFree = 0.0, From ed81530cde981826189286abc4176f2f9b64857b Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 9 Feb 2023 12:57:33 -0500 Subject: [PATCH 42/67] PC - bug fixes to use np instead of individual numpy modules. add aliases to ndvi, column ozone, and pixel elevation --- src/Applications/GAAS_App/vx04_nnr.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/Applications/GAAS_App/vx04_nnr.py b/src/Applications/GAAS_App/vx04_nnr.py index fc3c6905..64d9deee 100644 --- a/src/Applications/GAAS_App/vx04_nnr.py +++ b/src/Applications/GAAS_App/vx04_nnr.py @@ -11,6 +11,10 @@ from ffnet import loadnet import numpy as np +ALIAS = dict( TOA_NDVI = 'ndvi', + Total_Column_Ozone = 'colO3', + Precipitable_Water = 'water', + ) # Translate Inputs between NNR and MODIS classes # ----------------------------------------------- @@ -103,7 +107,7 @@ def __init__(self,l2_path,sat,algo,syn_time,aer_x, Vx04_L2.__init__(self,Files,algo,syn_time=syn_time,nsyn=nsyn, only_good=True, SDS=SDS, - alias=None, + alias=ALIAS, Verb=verbose, anet_wav=True) @@ -131,7 +135,7 @@ def __init__(self,l2_path,sat,algo,syn_time,aer_x, # Create attribute for holding NNR predicted AOD # ---------------------------------------------- - self.aod_ = MISSING * ones((self.nobs,len(self.channels))) + self.aod_ = MISSING * np.ones((self.nobs,len(self.channels))) # Make sure same good AOD is kept for gridding # -------------------------------------------- @@ -143,9 +147,8 @@ def __init__(self,l2_path,sat,algo,syn_time,aer_x, # Angle transforms: for NN calculations we work with cosine of angles # ------------------------------------------------------------------- self.ScatteringAngle = np.cos(self.ScatteringAngle*np.pi/180.0) - self.SensorAzimuth = np.cos(self.SensorAzimuth*np.pi/180.0) + self.RelativeAzimuth = np.cos(self.RelativeAzimuth*np.pi/180.0) self.SensorZenith = np.cos(self.SensorZenith*np.pi/180.0) - self.SolarAzimuth = np.cos(self.SolarAzimuth*np.pi/180.0) self.SolarZenith = np.cos(self.SolarZenith*np.pi/180.0) self.GlintAngle = np.cos(self.GlintAngle*np.pi/180.0) @@ -331,15 +334,15 @@ def apply(self,nnFile): if doAE: for i,targetName in enumerate(self.net.TargetNames): if 'Tau' in targetName: - name, base_wav = TranslateTarget[targetName] + name, base_wav = TranslateTarget(targetName) base_wav = np.float(base_wav) base_tau = targets[:,i] if self.net.laod: - base_tau = exp(base_tau) - 0.01 # inverse + base_tau = np.exp(base_tau) - 0.01 # inverse for i,targetName in enumerate(self.net.TargetNames): if 'AE' in targetName: AE = targets[:,i] - name, wav = TranslateTarget[targetName] + name, wav = TranslateTarget(targetName) wav = np.float(wav) data = base_tau*np.exp(-1.*AE*np.log(wav/base_wav)) if self.net.laod: @@ -355,7 +358,7 @@ def apply(self,nnFile): k = list(self.channels).index(ch) # index of channel except: # add new target channel to end - self.channels += (ch,) + self.channels = np.append(self.channels,ch) self.aod = np.append(self.aod,MISSING*np.ones((self.nobs,1)),axis=1) self.aod_ = np.append(self.aod_,MISSING*np.ones((self.nobs,1)),axis=1) From a7ca0115463450a91d9c84a97a069e417c98e38c Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 9 Feb 2023 12:58:30 -0500 Subject: [PATCH 43/67] PC - add 490 channel to aeronet read --- src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/aeronet.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/aeronet.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/aeronet.py index 0ef723a7..2d6b074d 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/aeronet.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/aeronet.py @@ -197,7 +197,18 @@ def __init__ (self,Path,version=2,Vars=VARS,Verbose=False): aot_470 = _updAOT(aot_470,aot_470a) self.AOT_470 = aot_470[:] # update undefs with interpolated values - + + # Interpolate AOT to 490 nm if needed + # ----------------------------------- + aot_490 = self.AOT_490.copy() + aot_490a = aodInterpAngs(490.,self.AOT_443,self.AOT_500,443.,500.) + aot_490b = aodInterpAngs(490.,self.AOT_440,self.AOT_500,440.,500.) + + aot_490 = _updAOT(aot_490,aot_490a) + aot_490 = _updAOT(aot_490,aot_490b) + + self.AOT_490 = aot_490[:] # update undefs with interpolated values + # Create timedate # --------------- self.tyme = array([ isoparse('-'.join(d.split(':')[-1::-1])+'T'+t) for d,t in zip(self.Date,self.Time) ]) From e8d9b60e0028ca8b01bdee914df294dae17791b0 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 9 Feb 2023 13:01:44 -0500 Subject: [PATCH 44/67] PC - make use from KX values for ODS files. other bug fixes --- .../GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 58 ++++++++++--------- 1 file changed, 30 insertions(+), 28 deletions(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index f1d0c097..100600ee 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -130,13 +130,13 @@ Spectral_TOA_Reflectance_Land = 'reflectance', TOA_NDVI = 'NDVI', Spectral_Single_Scattering_Albedo_Land = 'ssa', - Algorithm_Flag_Land = 'aflag', + Algorithm_Flag_Land = 'algflag', Aerosol_Type_Land = 'atype', Angstrom_Exponent_Land = 'angstrom', Spectral_Aerosol_Optical_Thickness_Ocean = 'aod', Aerosol_Optical_Thickness_550_Ocean = 'aod550', Spectral_TOA_Reflectance_Ocean = 'reflectance', - Algorithm_Flag_Ocean = 'aflag', + Algorithm_Flag_Ocean = 'algflag', Aerosol_Type_Ocean = 'atype', Angstrom_Exponent_Ocean = 'angstrom', Number_Of_Pixels_Used_Ocean = 'npixels_used', @@ -155,10 +155,10 @@ translate_sat = {'Suomi-NPP': 'SNPP'} -KX = dict ( SNPP_DT_OCEAN = 301, - SNPP_DT_LAND = 302, - SNPP_DB_OCEAN = 310, - SNPP_DB_LAND = 311, +KX = dict ( SNPP_DT_OCEAN = 401, + SNPP_DT_LAND = 402, + SNPP_DB_OCEAN = 403, + SNPP_DB_LAND = 404, ) KT = dict ( AOD = 45, ) @@ -354,15 +354,28 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, elif self.algo == 'DT_OCEAN': self.rChannels = np.array([480.,550.,670.,860.,1240.,1600.,2250.]) + # Concatenate AOD channels for Deep Blue + # -------------------------------------- + if self.algo == 'DB_LAND': + try: + self.aod = np.ones((self.nobs,4)) + self.aod[:,0] = self.aod3ch[:,0] + self.aod[:,1] = self.aod3ch[:,1] + self.aod[:,2] = self.aod550[:] + self.aod[:,3] = self.aod3ch[:,2] + except: + pass # don't fuss, aod3ch may not have been read + + # Angstrom interpolate to AERONET wavelengths for ease of comparison # and to faciliate NNR ODS files - # DT: 480 --> 490, 860 --> 870, 1600 --> 1640 - # DB: 488 --> 490, 865 --> 870 + # DT: 480 --> 470, 860 --> 870, 1600 --> 1640 + # DB: 488 --> 470, 865 --> 870 # ------------------------------------------------------------------ self.anet_wav = anet_wav if anet_wav: if 'DT' in self.algo: - tch = 490 + tch = 470 ch = 550 ich = np.argmin(np.abs(self.channels - ch)) @@ -373,7 +386,7 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, aod_ = self.aodInterpAngs(tch,aod480,aod550,480,550) self.channels[ich] = tch - self.aod[:.ich] = aod_ + self.aod[:,ich] = aod_ elif self.algo == 'DT_OCEAN': tch = 870 @@ -387,7 +400,7 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, aod_ = self.aodInterpAngs(tch,aod670,aod860,670,860) self.channels[ich] = tch - self.aod[:.ich] = aod_ + self.aod[:,ich] = aod_ tch = 1640 @@ -400,10 +413,10 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, aod_ = self.aodInterpAngs(tch,aod1600,aod2250,1600,2250) self.channels[ich] = tch - self.aod[:.ich] = aod_ + self.aod[:,ich] = aod_ elif 'DB' in self.algo: - tch = 490 + tch = 470 ch = 550 ich = np.argmin(np.abs(self.channels - ch)) @@ -414,7 +427,7 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, aod_ = self.aodInterpAngs(tch,aod488,aod550,488,550) self.channels[ich] = tch - self.aod[:.ich] = aod_ + self.aod[:,ich] = aod_ elif self.algo == 'DB_OCEAN': tch = 870 @@ -428,7 +441,7 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, aod_ = self.aodInterpAngs(tch,aod670,aod865,670,865) self.channels[ich] = tch - self.aod[:.ich] = aod_ + self.aod[:,ich] = aod_ @@ -448,17 +461,6 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, self.nhms = 10000 * syn_time.hour + 100*syn_time.minute + syn_time.second self.nsyn = nsyn # number of synoptic times per day - # Concatenate AOD channels for Deep Blue - # -------------------------------------- - if self.algo == 'DB_LAND': - try: - self.aod = np.ones((self.nobs,4)) - self.aod[:,0] = self.aod3ch[:,0] - self.aod[:,1] = self.aod3ch[:,1] - self.aod[:,2] = self.aod550[:] - self.aod[:,3] = self.aod3ch[:,2] - except: - pass # don't fuss, aod3ch may not have been read # Create a pseudo cloud fraction for Deep Blue if Algo == 'DB': @@ -466,7 +468,7 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, #--- - def aodInterpAngs(lambda_,tau1,tau2,lambda1,lambda2): + def aodInterpAngs(self,lambda_,tau1,tau2,lambda1,lambda2): """ Angstrom-interpolated AOD. lambda_ = target wavelength for AOD output @@ -726,7 +728,7 @@ def writeODS(self,filename=None,dir='.',expid=None,channels=None, ks = np.arange(ns) + 1 for ch in channels: I = range(i,i+ns) - j = channels.index(ch) + j = list(channels).index(ch) # index of channel ods.ks[I] = ks ods.lat[I] = self.lat[:] ods.lon[I] = self.lon[:] From a5c80c8c9002df0ed4c7afe5437bb4805bcb3b48 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 9 Feb 2023 19:24:24 -0500 Subject: [PATCH 45/67] PC - when angstrom interpolating, protecta against case where all obs are <0 --- src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index 100600ee..92fe5770 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -473,10 +473,11 @@ def aodInterpAngs(self,lambda_,tau1,tau2,lambda1,lambda2): Angstrom-interpolated AOD. lambda_ = target wavelength for AOD output """ - I = (tau1>0) & (tau2>0) - angstrom = -np.log(tau1[I]/tau2[I])/np.log(lambda1/lambda2) tau = MISSING * np.ones(len(tau1)) - tau[I] = tau2[I] * (lambda2/lambda_)**angstrom + I = (tau1>0) & (tau2>0) + if any(I): + angstrom = -np.log(tau1[I]/tau2[I])/np.log(float(lambda1)/float(lambda2)) + tau[I] = tau2[I] * (lambda2/lambda_)**angstrom return tau #--- From b6d533a6e6f76d9e8079319598825b284e734f6d Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 18 May 2023 15:37:38 -0400 Subject: [PATCH 46/67] PC - update changelog to develop v2.0.1 --- CHANGELOG.md | 41 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 39 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1651b68f..8bfb6042 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,8 +8,44 @@ ### Changed +## [v2.0.1] - 2023-05-17 + +### Added + + +### Fixed + +- fix in trj_sampler for ICARTT files that have missing values in the location/time +- fix in py_ods that was trying to do element wise comparisons on tuples and lists + +### Changed + +- updated GMAOpyobs to v1.0.2 +- updates to NNR training code to handle angstrom exponent targets +- some other minor changes to NNR training code to be python3 compliant + +## [v2.0.0] - 2023-05-17 + +### Added + +- Added changelog enforcer + +### Changed + +- Updated to use latest components (matching GEOSgcm v10.25.0) +- Wholesale conversion to python3. Changes to all python scripts, and CMakelists with f2py compilation +- the components.yaml now imports the GMAOpyobs repository, which replaces the legacy GMAO_pyobs directory +- The GMAO_pyobs3 directory is no longer needed +- A new directory called GMAO_aeropyobs has been added to share. This includes python utilities that require MAPL or GFIO that were previously in GMAO_pyobs. +- There may still be issues related to byte-to-string conversions that require fixing, but will need to be fixed as they are encountered. +- Added GMAOpyobs v1.0.1 to components.yaml +- Updated MAPL to v2.39.9 + ## [v1.0.0] - 2022-06-17 -Converted to mepo, cmake, and v3 giant files + +### Changed + +- Converted to mepo, cmake, and v3 giant files ### Added @@ -18,12 +54,13 @@ Converted to mepo, cmake, and v3 giant files - legacy MAPL_Base/Python/MAPL into latest MAPL on develop branch ### Fixed + - Moved src/GMAO_Shared to src/Shared/GMAO_Shared - obs_aod codes now work with v3 giant files ## [0.99.0] - 2020-03-23 -Intial GitHub release. +Initial GitHub release. ### Added - Open source license From 8c755639791c3739f7618cc9589aa296ac86bf26 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 18 May 2023 15:41:52 -0400 Subject: [PATCH 47/67] updates ABC cmakelists to f2py3 --- src/Components/misc/obs_aod/ABC/CMakeLists.txt | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/CMakeLists.txt b/src/Components/misc/obs_aod/ABC/CMakeLists.txt index 269916d7..68b498ad 100644 --- a/src/Components/misc/obs_aod/ABC/CMakeLists.txt +++ b/src/Components/misc/obs_aod/ABC/CMakeLists.txt @@ -1,7 +1,7 @@ esma_set_this () -find_package(F2PY2 REQUIRED) -esma_add_f2py2_module (VLIDORT_BRDF_ABC_ +find_package(F2PY3 REQUIRED) +esma_add_f2py3_module (VLIDORT_BRDF_ABC_ SOURCES VLIDORT_BRDF_ABC_py.F90 DESTINATION lib/Python/${this} LIBRARIES VLIDORT90 @@ -37,7 +37,6 @@ set (PYFILES evaluate_cdr.py giant.py giant_viirs.py - man.py mcd43c.py nn.py reduce_avhrr.py From e4fd98db5dcd35d458efb7415069b986d23b3ea9 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 19 May 2023 11:30:28 -0400 Subject: [PATCH 48/67] PC - ran 2to3 on viirs_l2a.py --- src/Applications/GAAS_App/viirs_l2a.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Applications/GAAS_App/viirs_l2a.py b/src/Applications/GAAS_App/viirs_l2a.py index a86974c5..c9babc66 100755 --- a/src/Applications/GAAS_App/viirs_l2a.py +++ b/src/Applications/GAAS_App/viirs_l2a.py @@ -57,10 +57,10 @@ for ident in cf('VIIRS_L2A_IDENTS').split(','): coll = Coll[i] cmd = "vx04_l2a.py %s --collection=%s %s %s "%(Options,Coll[i],ident,isotime) - print cmd + print(cmd) if not options.dryrun: if system(cmd): - raise ValueError, "vx04_l2a.py failed for %s on %s "%(ident,isotime) + raise ValueError("vx04_l2a.py failed for %s on %s "%(ident,isotime)) i += 1 From 0bf7e526524079429b92cb95cbf700c7b277586c Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 19 May 2023 11:37:52 -0400 Subject: [PATCH 49/67] PC ran 2to3 on vx04_l2a.py --- src/Applications/GAAS_App/vx04_l2a.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/Applications/GAAS_App/vx04_l2a.py b/src/Applications/GAAS_App/vx04_l2a.py index 7a2e3845..0cb0a7d5 100755 --- a/src/Applications/GAAS_App/vx04_l2a.py +++ b/src/Applications/GAAS_App/vx04_l2a.py @@ -40,7 +40,7 @@ def makethis_dir(filename): if path != '': rc = os.system('mkdir -p '+path) if rc: - raise IOError, "could not create directory "+path + raise IOError("could not create directory "+path) #--------------------------------------------------------------------- @@ -137,10 +137,10 @@ def makethis_dir(filename): # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - if options.verbose: - print "" - print " VIIRS Level 2A Processing" - print " -------------------------" - print "" + print("") + print(" VIIRS Level 2A Processing") + print(" -------------------------") + print("") t0 = clock() # Time variables @@ -156,8 +156,8 @@ def makethis_dir(filename): out_file = strTemplate(out_tmpl,expid=options.expid,nymd=nymd,nhms=nhms) name, ext = os.path.splitext(out_file) if os.path.exists(out_file) and (options.force is not True): - print "vx04_l2a: Output Gridded file <%s> exists --- cannot proceed."%out_file - raise IOError, "Specify --force to overwrite existing output file." + print("vx04_l2a: Output Gridded file <%s> exists --- cannot proceed."%out_file) + raise IOError("Specify --force to overwrite existing output file.") if os.path.exists(out_file) and options.force: os.remove(out_file) @@ -167,8 +167,8 @@ def makethis_dir(filename): ods_tmpl = ods_tmpl.replace('%coll',options.coll).replace('%prod',prod).replace('%algo',algo).replace('%lev','2').replace('%ext','ods') ods_file = strTemplate(ods_tmpl,expid=options.expid,nymd=nymd,nhms=nhms) if os.path.exists(ods_file) and (options.force is not True): - print "vxd04_l2a: Output ODS file <%s> exists --- cannot proceed."%ods_file - raise IOError, "Specify --force to overwrite existing output file." + print("vxd04_l2a: Output ODS file <%s> exists --- cannot proceed."%ods_file) + raise IOError("Specify --force to overwrite existing output file.") if os.path.exists(ods_file) and options.force: os.remove(ods_file) @@ -183,7 +183,7 @@ def makethis_dir(filename): # VIIRS Level 2 NNR Aerosol Retrievals # ------------------------------------ if options.verbose: - print "NNR Retrieving %s %s on "%(sat,algo.upper()),syn_time + print("NNR Retrieving %s %s on "%(sat,algo.upper()),syn_time) viirs = Vx04_NNR(options.l2_path,sat,algo.upper(),syn_time,aer_x, coll=options.coll, @@ -193,11 +193,11 @@ def makethis_dir(filename): verbose=options.verbose) if viirs.nobs < 1: if options.verbose: - print 'WARNING: no observation for this time in file <%s>'%ods_file + print('WARNING: no observation for this time in file <%s>'%ods_file) elif any(viirs.iGood) == False: if options.verbose: - print 'WARNING: no GOOD observation for this time in file <%s>'%ods_file + print('WARNING: no GOOD observation for this time in file <%s>'%ods_file) viirs.nobs = 0 nn_file = options.nn_file.replace('%ident',ident) @@ -213,7 +213,7 @@ def makethis_dir(filename): warnings.warn('cannot create empty output file <%s>'%ods_file) else: if options.verbose: - print "[w] Wrote empty ODS file "+ods_file + print("[w] Wrote empty ODS file "+ods_file) # Write gridded output file (revised channels only) # ------------------------------------------------- From 91da9f71c0b7f1f56de5b268d128ce1b1bbacf83 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 19 May 2023 11:38:59 -0400 Subject: [PATCH 50/67] PC - ran 2to3 on vx04_nnr.py --- src/Applications/GAAS_App/vx04_nnr.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/Applications/GAAS_App/vx04_nnr.py b/src/Applications/GAAS_App/vx04_nnr.py index 64d9deee..7291cd6d 100644 --- a/src/Applications/GAAS_App/vx04_nnr.py +++ b/src/Applications/GAAS_App/vx04_nnr.py @@ -130,7 +130,7 @@ def __init__(self,l2_path,sat,algo,syn_time,aer_x, self.iGood = self.iGood & (self.GlintAngle > glint_thresh) if np.any(self.iGood) == False: - print "WARNING: Strange, no good obs left to work with" + print("WARNING: Strange, no good obs left to work with") return # Create attribute for holding NNR predicted AOD @@ -214,7 +214,7 @@ def sampleFile(self, inFile, npzFile=None, onlyVars=None, Verbose=False): if fh.lm == 1: timeInterp = False # no time interpolation in this case else: - raise ValueError, "cannot handle files with more tha 1 time, use ctl instead" + raise ValueError("cannot handle files with more tha 1 time, use ctl instead") else: fh = GFIOctl(inFile) # open timeseries timeInterp = True # perform time interpolation @@ -233,7 +233,7 @@ def sampleFile(self, inFile, npzFile=None, onlyVars=None, Verbose=False): # --------------------------- for v in onlyVars: if Verbose: - print "<> Sampling ", v + print("<> Sampling ", v) if timeInterp: var = fh.sample(v,lons,lats,tymes,Verbose=Verbose) else: @@ -246,7 +246,7 @@ def sampleFile(self, inFile, npzFile=None, onlyVars=None, Verbose=False): var = var.T # shape should be (nobs,nz) self.sample.__dict__[v] = var else: - raise IndexError, 'variable <%s> has rank = %d'%(v,len(var.shape)) + raise IndexError('variable <%s> has rank = %d'%(v,len(var.shape))) if npzFile is not None: savez(npzFile,**self.sample.__dict__) @@ -270,7 +270,7 @@ def _getInputs(self): iName = TranslateInput(inputName) if self.verbose>0: - print 'Getting NN input ',iName + print('Getting NN input ',iName) # Retrieve input # -------------- @@ -288,7 +288,7 @@ def _getInputs(self): input = self.__dict__[name][:] else: - raise ValueError, "strange, len(iName)=%d"%len(iName) + raise ValueError("strange, len(iName)=%d"%len(iName)) # Concatenate Inputs # ------------------ @@ -369,9 +369,9 @@ def apply(self,nnFile): name, ch = TranslateTarget(targetName) if self.verbose>0: if self.net.laod: - print "NN Retrieving log(AOD+0.01) at %dnm "%ch + print("NN Retrieving log(AOD+0.01) at %dnm "%ch) else: - print "NN Retrieving AOD at %dnm "%ch + print("NN Retrieving AOD at %dnm "%ch) k = list(self.channels).index(ch) # index of channel self.channels_ = self.channels_ + [ch,] if self.net.laod: @@ -421,7 +421,7 @@ def apply(self,nnFile): coll = '002' aer_x = '/nobackup/NNR/Misc/tavg1_2d_aer_Nx' - syn_time = datetime(2012,03,01,00,0,0) + syn_time = datetime(2012,0o3,0o1,00,0,0) if algo == 'DB_LAND': nn_file = '/nobackup/NNR/Net/VIIRS/nnr_001.vsnppdbl_Tau.net' From 6c995002ab05a74a40b25f0c918b74708d320324 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 19 May 2023 11:56:22 -0400 Subject: [PATCH 51/67] PC - ran 2to3 on abc_c6_aux.py --- src/Components/misc/obs_aod/ABC/abc_c6_aux.py | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/abc_c6_aux.py b/src/Components/misc/obs_aod/ABC/abc_c6_aux.py index 0b9f773a..f7460cf6 100644 --- a/src/Components/misc/obs_aod/ABC/abc_c6_aux.py +++ b/src/Components/misc/obs_aod/ABC/abc_c6_aux.py @@ -173,7 +173,7 @@ def SummarizeCombinations(mxd,Input_nnr,yrange=None,sortname='rmse'): nblocks = [len(group) for group in mxd.combgroups] nblocks.insert(0,0) - print "MASTERLIST", masterlist + print("MASTERLIST", masterlist) #-------------- # Default sort by mean RMSE of first target @@ -393,7 +393,7 @@ def make_plots(mxd,expid,ident,I=None): for t in range(mxd.nTarget): name = mxd.Target[t][1:] if name != refname: - print 't,wav',t,name[3:] + print('t,wav',t,name[3:]) wav = float(name[3:]) tt = np.exp(targets[:,t]) # keep the + 0.01 to handle negatives in MODIS data rr = np.exp(results[:,t]) # keep the + 0.01 to handle negatives in MODIS data @@ -419,7 +419,7 @@ def make_plots(mxd,expid,ident,I=None): name = 'm'+mxd.Target[t][1:] if name in mxd.__dict__: if name != refname: - print 'orig t,wav',t,name[4:] + print('orig t,wav',t,name[4:]) wav = float(name[4:]) oo = mxd.__dict__[name][I] + 0.01 # add 0.01 to handle negatives tt = np.exp(targets[:,t]) # keep + 0.01 to handle negatives @@ -528,7 +528,7 @@ def make_plots_angstrom(mxd,expid,ident,I=None): else: name = mxd.Target[t][1:] if name != refname: - print 't,wav',t,name[3:] + print('t,wav',t,name[3:]) wav = float(name[3:]) tt = np.exp(targets[:,t]) # keep the + 0.01 to handle negatives in MODIS data rr = np.exp(results[:,t]) # keep the + 0.01 to handle negatives in MODIS data @@ -559,7 +559,7 @@ def make_plots_angstrom(mxd,expid,ident,I=None): name = 'm'+mxd.Target[t][1:] if name in mxd.__dict__: if name != refname: - print 'orig t,wav',t,name[4:] + print('orig t,wav',t,name[4:]) wav = float(name[4:]) oo = mxd.__dict__[name][I] + 0.01 # add 0.01 to handle negatives # protect against interpolated wavelengths that might have -9999 @@ -667,8 +667,8 @@ def make_error_pdfs(mxd,Input,expid,ident,K=None,I=None,Title=None,netfileRoot=N mod04RMSE.append(rmse(original[k],targets[k])) nnrRMSE.append(rmse(results[k],targets[k])) - print 'mod04RMSE',mod04RMSE - print 'nnrRMSE',nnrRMSE + print('mod04RMSE',mod04RMSE) + print('nnrRMSE',nnrRMSE) mod04RMSE = np.mean(mod04RMSE) nnrRMSE = np.mean(nnrRMSE) @@ -812,9 +812,9 @@ def make_error_pdfs_int(mxd,Input,expid,ident,K=None,I=None,Title=None,netfileRo dbmod04RMSE.append(rmse(dboriginal[k],targets[k])) nnrRMSE.append(rmse(results[k],targets[k])) - print 'mod04RMSE',mod04RMSE - print 'dbmod04RMSE',dbmod04RMSE - print 'nnrRMSE',nnrRMSE + print('mod04RMSE',mod04RMSE) + print('dbmod04RMSE',dbmod04RMSE) + print('nnrRMSE',nnrRMSE) mod04RMSE = np.mean(mod04RMSE) dbmod04RMSE = np.mean(dbmod04RMSE) nnrRMSE = np.mean(nnrRMSE) @@ -989,10 +989,10 @@ def make_error_pdfs_dbdt(mxd,mxd2,Input,expid,ident,K=None,I=None,Title=None, nnrRMSE.append(rmse(results[k],targets[k])) nnrRMSE2.append(rmse(results2[k],targets[k])) - print 'mod04RMSE',mod04RMSE - print 'dbmod04RMSE',dbmod04RMSE - print 'nnrRMSE',nnrRMSE - print 'nnrRMSE2',nnrRMSE2 + print('mod04RMSE',mod04RMSE) + print('dbmod04RMSE',dbmod04RMSE) + print('nnrRMSE',nnrRMSE) + print('nnrRMSE2',nnrRMSE2) mod04RMSE = np.mean(mod04RMSE) dbmod04RMSE = np.mean(dbmod04RMSE) nnrRMSE = np.mean(nnrRMSE) From 0337da097739a87c34449ecb27c358c80f2d0687 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 19 May 2023 11:57:16 -0400 Subject: [PATCH 52/67] PC - brought in changes from develop --- src/Components/misc/obs_aod/ABC/abc_c6_aux.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/abc_c6_aux.py b/src/Components/misc/obs_aod/ABC/abc_c6_aux.py index f7460cf6..9710cff0 100644 --- a/src/Components/misc/obs_aod/ABC/abc_c6_aux.py +++ b/src/Components/misc/obs_aod/ABC/abc_c6_aux.py @@ -34,17 +34,17 @@ def boxplot_imshow(data,plottype,blocks,masterlist,title,filename, params = {'mathtext.default': 'regular' } plt.rcParams.update(params) - if plottype is 'box': + if plottype == 'box': bp = plt.boxplot(data,showfliers=False,showbox=True,whis='range', whiskerprops={'linestyle':'-'}) plt.setp(bp['boxes'], color='black') plt.setp(bp['whiskers'], color='black') plt.plot([0,ncomb+0.5],[np.median(data[:,0]),np.median(data[:,0])],color='b',ls='--',zorder=5) - elif plottype is 'scatter': + elif plottype == 'scatter': scat = np.mean(data,axis=0) plt.plot(np.arange(ncomb)+1,scat,'rD') plt.plot([0,ncomb+0.5],[np.mean(data[:,0]),np.mean(data[:,0])],color='b',ls='--',zorder=5,markersize=5) - elif plottype is 'errorbar': + elif plottype == 'errorbar': scat = np.mean(data,axis=0) yerr_max = np.abs(np.max(data,axis=0)-scat) yerr_min = np.abs(np.min(data,axis=0)-scat) @@ -1115,6 +1115,7 @@ def TestStats(mxd,K,C): # regression[*][0:2] = slope, intercept, r-value # out.shape = [ntestobs,nTarget] out, reg = mxd.test(iprint=False) + reg = np.array(reg) mxd.nnr.slope[k,c,:] = reg[:,0] mxd.nnr.intercept[k,c,:] = reg[:,1] From b4d1ce724a822bc8a660b30ef0f0fe6b42b07d7f Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 19 May 2023 12:37:44 -0400 Subject: [PATCH 53/67] PC - ran 2to3 on abc_viirs.py --- src/Components/misc/obs_aod/ABC/abc_viirs.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/abc_viirs.py b/src/Components/misc/obs_aod/ABC/abc_viirs.py index f13f2ea2..3ae10f87 100644 --- a/src/Components/misc/obs_aod/ABC/abc_viirs.py +++ b/src/Components/misc/obs_aod/ABC/abc_viirs.py @@ -25,6 +25,7 @@ from abc_c6_aux import make_plots, make_plots_angstrom, TestStats, SummaryPDFs from brdf import rtlsReflectance from mcd43c import BRDF +from functools import reduce # ------ @@ -258,13 +259,13 @@ def outlierRemoval(self,outliers): if outliers > 0.: d = log(self.mTau550[self.iValid]+0.01) - log(self.aTau550[self.iValid]+0.01) if self.verbose>0: - print "Outlier removal: %d sig_d = %f nGood=%d "%(-1,std(d),d.size) + print("Outlier removal: %d sig_d = %f nGood=%d "%(-1,std(d),d.size)) for iter in range(3): iValid = (abs(d)<outliers*std(d)) self.iValid[self.iValid] = iValid d = log(self.mTau550[self.iValid]+0.01) - log(self.aTau550[self.iValid]+0.01) if self.verbose>0: - print "Outlier removal: %d sig_d = %f nGood=%d "%(iter,std(d),d.size) + print("Outlier removal: %d sig_d = %f nGood=%d "%(iter,std(d),d.size)) def angleTranform(self): # Angle transforms: for NN work we work with cosine of angles @@ -681,9 +682,9 @@ def _train(mxd,expid,c): Input = mxd.comblist[c] Target = mxd.Target - print "-"*80 - print "--> nHidden = ", nHidden - print "--> Inputs = ", Input + print("-"*80) + print("--> nHidden = ", nHidden) + print("--> Inputs = ", Input) n = cpu_count() kwargs = {'nproc' : n} @@ -721,7 +722,7 @@ def _test(mxd,expid,c,plotting=True): if found: break if not found: - print '{} not found. Need to train this combinatin of inputs'.format(netFile) + print('{} not found. Need to train this combinatin of inputs'.format(netFile)) raise else: invars = mxd.comblist[0] @@ -751,13 +752,13 @@ def _test(mxd,expid,c,plotting=True): netFile = outdir+"/"+".".join(invars)+'.k={}_Tau.net'.format(str(k)) mxd.loadnet(netFile) found = True - print 'found file',netFile + print('found file',netFile) break except: pass if not found: - print '{} not found. Need to train this combinatin of inputs'.format(netFile) + print('{} not found. Need to train this combinatin of inputs'.format(netFile)) raise else: invars = mxd.comblist[0] From 9740beaabfaeccff1ddc0aec7701dce0a6bc4c94 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 19 May 2023 12:40:25 -0400 Subject: [PATCH 54/67] PC - ran 2to3 on giant_viirs.py --- src/Components/misc/obs_aod/ABC/giant_viirs.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/giant_viirs.py b/src/Components/misc/obs_aod/ABC/giant_viirs.py index dc2b935e..16b837c1 100644 --- a/src/Components/misc/obs_aod/ABC/giant_viirs.py +++ b/src/Components/misc/obs_aod/ABC/giant_viirs.py @@ -280,9 +280,9 @@ def __init__ (self,filename,xVars=(),only_good=True,tymemax=None): # Read in variables # ----------------- - print 'filename ',filename + print('filename ',filename) nc = Dataset(filename) - Alias = self.ALIAS.keys() + Alias = list(self.ALIAS.keys()) self.giantList =[] for name in Names: data = nc.variables[name][:] @@ -301,7 +301,7 @@ def __init__ (self,filename,xVars=(),only_good=True,tymemax=None): # ---------------- # new files have an ISO_DateTime variable nc = Dataset(filename) - if 'ISO_DateTime' in nc.variables.keys(): + if 'ISO_DateTime' in list(nc.variables.keys()): try: iso = nc.variables['ISO_DateTime'][:] self.tyme = array([isoparse(''.join(array(t))) for t in iso]) @@ -368,7 +368,7 @@ def reduce(self,I): """ for name in self.giantList: q = self.__dict__[name] - print "{} Reducing "+name,q.shape + print("{} Reducing "+name,q.shape) self.__dict__[name] = q[I] self.nobs = len(self.lon) @@ -557,7 +557,7 @@ def sampleFile(self, inFile, npzFile=None, onlyVars=None, Verbose=False, clmYear if fh.lm == 1: timeInterp = False # no time interpolation in this case else: - raise ValueError, "cannot handle files with more tha 1 time, use ctl instead" + raise ValueError("cannot handle files with more tha 1 time, use ctl instead") else: fh = GFIOctl(inFile) # open timeseries timeInterp = True # perform time interpolation @@ -575,12 +575,12 @@ def sampleFile(self, inFile, npzFile=None, onlyVars=None, Verbose=False, clmYear else: tymes = array([t + timedelta(days=365*(clmYear-t.year)) for t in self.tyme]) - print 'trange',tymes.min(),tymes.max() + print('trange',tymes.min(),tymes.max()) # Loop over variables on file # --------------------------- for v in onlyVars: if Verbose: - print "<> Sampling ", v + print("<> Sampling ", v) if timeInterp: var = fh.sample(v,lons,lats,tymes,Verbose=Verbose) else: @@ -591,7 +591,7 @@ def sampleFile(self, inFile, npzFile=None, onlyVars=None, Verbose=False, clmYear var = var.T # shape should be (nobs,nz) self.sample.__dict__[v] = var else: - raise IndexError, 'variable <%s> has rank = %d'%(v,len(var.shape)) + raise IndexError('variable <%s> has rank = %d'%(v,len(var.shape))) if npzFile is not None: savez(npzFile,**self.sample.__dict__) @@ -603,7 +603,7 @@ def sampleLoadz(self,npzFile): from grads.gahandle import GaHandle self.sample = GaHandle(npzFile) npz = load(npzFile) - for v in npz.keys(): + for v in list(npz.keys()): self.sample.__dict__[v] = npz[v] #--- From 3006103217fa5fb83933207756abb3534c8bb4c6 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 19 May 2023 12:51:47 -0400 Subject: [PATCH 55/67] pc - ran2to3 on nn.py --- src/Components/misc/obs_aod/ABC/nn.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/nn.py b/src/Components/misc/obs_aod/ABC/nn.py index 070a41be..4b82867b 100644 --- a/src/Components/misc/obs_aod/ABC/nn.py +++ b/src/Components/misc/obs_aod/ABC/nn.py @@ -86,8 +86,8 @@ def train (self,Input=None,Target=None,nHidden=200,maxfun=2550,biases=True, # ----- bounds = [bounds]*self.net.get_params()['conec'].shape[0] if self.verbose>0: - print "Starting training with %s inputs and %s targets"\ - %(str(inputs.shape),str(targets.shape)) + print("Starting training with %s inputs and %s targets"\ + %(str(inputs.shape),str(targets.shape))) self.net.train_tnc(inputs,targets, maxfun=maxfun,bounds=bounds,**kwargs) # self.net.train_bfgs(inputs,targets, maxfun=maxfun) @@ -160,22 +160,22 @@ def getInputs(self,I,Input=None): Returns: inputs """ if self.verbose: - print " " - print " Feature Min Max" - print " ------------------ ------- -------" + print(" ") + print(" Feature Min Max") + print(" ------------------ ------- -------") if Input==None: Input = self.Input inputs = self.__dict__[Input[0]][I] if self.verbose: - print "%20s %8.4f %8.4f"%(Input[0],inputs.min(),inputs.max()) + print("%20s %8.4f %8.4f"%(Input[0],inputs.min(),inputs.max())) for var in Input[1:]: q = self.__dict__[var][I] inputs = cat[inputs,q] if self.verbose: - print "%20s %8.4f %8.4f"%(var,q.min(),q.max()) + print("%20s %8.4f %8.4f"%(var,q.min(),q.max())) if self.verbose: - print " ------------------ ------- -------" - print "" + print(" ------------------ ------- -------") + print("") if len(inputs.shape) == 1: inputs.shape = (inputs.shape[0],1) @@ -277,7 +277,7 @@ def _plotKDE(x_values,y_values,x_bins=None,y_bins=None, Nx = len(x_bins) Ny = len(y_bins) - print "Evaluating 2D kernel on grid with (Nx,Ny)=(%d,%d) ..."%(Nx,Ny) + print("Evaluating 2D kernel on grid with (Nx,Ny)=(%d,%d) ..."%(Nx,Ny)) kernel = stats.kde.gaussian_kde(_cat2(x_values,y_values)) X, Y = meshgrid(x_bins,y_bins) # each has shape (Ny,Nx) Z = kernel(_cat2(X,Y)) # shape is (Ny*Nx) From 0fd0e9a5154c973c5e38b198d77dad0cc143609d Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 19 May 2023 12:52:24 -0400 Subject: [PATCH 56/67] PC - add bounds variables fix in nn.py from develop --- src/Components/misc/obs_aod/ABC/nn.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Components/misc/obs_aod/ABC/nn.py b/src/Components/misc/obs_aod/ABC/nn.py index 4b82867b..dca3aea2 100644 --- a/src/Components/misc/obs_aod/ABC/nn.py +++ b/src/Components/misc/obs_aod/ABC/nn.py @@ -84,7 +84,7 @@ def train (self,Input=None,Target=None,nHidden=200,maxfun=2550,biases=True, # Train # ----- - bounds = [bounds]*self.net.get_params()['conec'].shape[0] + bounds = [bounds]*self.net.conec.shape[0] if self.verbose>0: print("Starting training with %s inputs and %s targets"\ %(str(inputs.shape),str(targets.shape))) From c522dbe99b4a090415b71ae31b52289414a69ebb Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 19 May 2023 13:02:33 -0400 Subject: [PATCH 57/67] PC - ran2to3 on aeronet.py --- .../GMAO_Shared/GMAO_pyobs/pyobs/aeronet.py | 38 +++++++++---------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/aeronet.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/aeronet.py index 2d6b074d..bb1f2f28 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/aeronet.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/aeronet.py @@ -11,11 +11,11 @@ from dateutil.parser import parse as isoparse from glob import glob -from npz import NPZ +from .npz import NPZ MISSING = -999. -BAD, MARGINAL, GOOD, BEST = range(4) +BAD, MARGINAL, GOOD, BEST = list(range(4)) VARS = ( 'AERONET_Site', 'Longitude', @@ -94,7 +94,7 @@ def __init__ (self,Path,version=2,Vars=VARS,Verbose=False): # ---------------------- if type(Path) is ListType: if len(Path) == 0: - print "WARNING: Empty AERONET object created" + print("WARNING: Empty AERONET object created") return else: if Path[-4:] == '.npz': # Special handling of npz files @@ -119,11 +119,11 @@ def __init__ (self,Path,version=2,Vars=VARS,Verbose=False): try: self.__dict__[var] = concatenate(self.__dict__[var]) except: - print "Failed concatenating "+var + print("Failed concatenating "+var) # Make aliases # ------------ - Alias = ALIAS.keys() + Alias = list(ALIAS.keys()) for var in self.Vars: if var in Alias: self.__dict__[ALIAS[var]] = self.__dict__[var] @@ -218,7 +218,7 @@ def __init__ (self,Path,version=2,Vars=VARS,Verbose=False): Locations = {} for loc in self.Location: Locations[loc] = 1 - self.Stations = Locations.keys() + self.Stations = list(Locations.keys()) # By default all is good if coordinates are ok # -------------------------------------------- @@ -239,7 +239,7 @@ def _readList(self,List): if os.path.isdir(item): self._readDir(item) elif os.path.isfile(item): self._readGranule(item) else: - print "%s is not a valid file or directory, ignoring it"%item + print("%s is not a valid file or directory, ignoring it"%item) #--- def _readDir(self,dir): """Recursively, look for files in directory.""" @@ -248,7 +248,7 @@ def _readDir(self,dir): if os.path.isdir(path): self._readDir(path) elif os.path.isfile(path): self._readGranule(path) else: - print "%s is not a valid file or directory, ignoring it"%item + print("%s is not a valid file or directory, ignoring it"%item) #--- def _readGranule(self,filename): @@ -272,7 +272,7 @@ def _readGranule(self,filename): i += 1 if self.columns == None: - raise ValueError, "Cannot find Column header" + raise ValueError("Cannot find Column header") # Read relevant columns from AERONET granule # ---------------------------------------- @@ -283,7 +283,7 @@ def _readGranule(self,filename): try: i = self.columns.index(name) except: - raise ValueError, "cannot find <%s> in file <%s>"%(name,filename) + raise ValueError("cannot find <%s> in file <%s>"%(name,filename)) self.iVars += (i,) if name=='Date': self.formats += ('S10',) @@ -325,7 +325,7 @@ def writeNPZ(self,npzFile,I=None): Writes out a NPZ file with the relevant variables. """ Vars = dict() - Nicknames = ALIAS.values() + Nicknames = list(ALIAS.values()) for name in self.__dict__: if name in Nicknames: continue # alias do not get reduced @@ -373,7 +373,7 @@ def writeODS(self, syn_tyme, filename=None, dir='.', expid='aeronet', nsyn=8): ods = ODS(nobs=nobs, kx=KX, kt=KT['AOD']) - ods.ks[:] = range(1,1+nobs) + ods.ks[:] = list(range(1,1+nobs)) ods.lat[:] = self.lat[I] ods.lon[:] = self.lon[I] ods.qch[:] = zeros(nobs).astype('int') @@ -386,8 +386,8 @@ def writeODS(self, syn_tyme, filename=None, dir='.', expid='aeronet', nsyn=8): ods.xm[:] = self.AOT_675[I] if self.verb: - print "[w] Writing file <"+filename+"> with %d observations at %dZ"%\ - (ods.nobs,nhms/10000) + print("[w] Writing file <"+filename+"> with %d observations at %dZ"%\ + (ods.nobs,nhms/10000)) ods.write(filename,nymd,nhms,nsyn=nsyn) @@ -493,7 +493,7 @@ def writeGridded(self, syn_tyme, binobs2d(self.lon[I],self.lat[I],self.AOT_550[I],im,jm,MISSING) ) if self.verb: - print "[w] Wrote file "+filename+" at %02dZ"%(nhms/10000) + print("[w] Wrote file "+filename+" at %02dZ"%(nhms/10000)) return filename @@ -530,7 +530,7 @@ def bin_hourly(self,x,stn): n = int(0.5+(tf-t0).total_seconds()/60.) t = array([t0 + i*dt for i in range(n)]) - print t0, tf + print(t0, tf) return t @@ -577,10 +577,10 @@ def retrieve( filename='aeronet.csv', %(filename,request) if verbose: - print cmd + print(cmd) if os.system(cmd): - raise ValueError, "Cannot retrieve request <%cmd>" + raise ValueError("Cannot retrieve request <%cmd>") #....................................................................................... @@ -637,7 +637,7 @@ def granules( tyme, def __test__(): - t = datetime(2008,07,01) + t = datetime(2008,0o7,0o1) one_hour = timedelta(seconds=60*60) # synoptic dt = 3 hours Files = granules(t,bracket='left') From 65589c260c20ab140b56cb912fe04e7bbc811f43 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 19 May 2023 13:04:03 -0400 Subject: [PATCH 58/67] PC - add chanegs from develop, fixing byte read and no longer use listType --- .../GMAO_Shared/GMAO_pyobs/pyobs/aeronet.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/aeronet.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/aeronet.py index bb1f2f28..5ffe3719 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/aeronet.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/aeronet.py @@ -4,7 +4,6 @@ import os import sys -from types import * from numpy import loadtxt, ones, zeros, savez, pi, log, concatenate, \ arange, savez, shape, array, linspace from datetime import datetime, timedelta @@ -22,7 +21,7 @@ 'Latitude', 'Elevation', 'Date', - 'Time', + 'Time', 'AOT_1640' , 'AOT_1020', 'AOT_870', @@ -86,13 +85,13 @@ def __init__ (self,Path,version=2,Vars=VARS,Verbose=False): self.Vars[i] = 'Site_Latitude' if 'Elevation' in self.Vars: i = self.Vars.index('Elevation') - self.Vars[i] = 'Site_Elevation' + self.Vars[i] = 'Site_Elevation' self.Vars = tuple(self.Vars) # Past is string or list # ---------------------- - if type(Path) is ListType: + if isinstance(Path, (list, tuple)): if len(Path) == 0: print("WARNING: Empty AERONET object created") return @@ -164,7 +163,7 @@ def __init__ (self,Path,version=2,Vars=VARS,Verbose=False): aot_670 = _updAOT(aot_670,self.AOT_667) # close enough aot_670 = _updAOT(aot_670,self.AOT_675) # close enough - self.AOT_670 = aot_670[:] # update undefs with interpolated values + self.AOT_670 = aot_670[:] # update undefs with interpolated values # Interpolate AOT to 660 nm # ----------------------------------- @@ -180,7 +179,7 @@ def __init__ (self,Path,version=2,Vars=VARS,Verbose=False): aot_660 = _updAOT(aot_660,aot_660b) aot_660 = _updAOT(aot_660,aot_660a) - self.AOT_660 = aot_660[:] # update undefs with interpolated values + self.AOT_660 = aot_660[:] # update undefs with interpolated values # Interpolate AOT to 470 nm if needed # ----------------------------------- @@ -286,11 +285,11 @@ def _readGranule(self,filename): raise ValueError("cannot find <%s> in file <%s>"%(name,filename)) self.iVars += (i,) if name=='Date': - self.formats += ('S10',) + self.formats += ('U10',) elif name=='Time': - self.formats += ('S8',) + self.formats += ('U8',) elif name=='AERONET_Site': - self.formats += ('S20',) + self.formats += ('U20',) else: self.converters[i] = _convert2Float self.formats += ('f4',) From 01492694c12cc6c84ba23ecb7261e2ffddceeaf3 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Fri, 19 May 2023 13:12:10 -0400 Subject: [PATCH 59/67] PC - ran 2to3 on vx04.py --- .../GMAO_Shared/GMAO_pyobs/pyobs/vx04.py | 72 +++++++++---------- 1 file changed, 36 insertions(+), 36 deletions(-) diff --git a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py index 92fe5770..b5573182 100644 --- a/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py +++ b/src/Shared/GMAO_Shared/GMAO_pyobs/pyobs/vx04.py @@ -24,7 +24,7 @@ pass from netCDF4 import Dataset -from bits import BITS +from .bits import BITS #--- @@ -215,7 +215,7 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, """ if algo not in ('DT_LAND', 'DT_OCEAN', 'DB_LAND', 'DB_OCEAN'): - raise ValueError, "invalid algorithm "+algo+" --- must be DT_LAND, DT_OCEAN, DB_LAND, DB_OCEAN" + raise ValueError("invalid algorithm "+algo+" --- must be DT_LAND, DT_OCEAN, DB_LAND, DB_OCEAN") # Initially are lists of numpy arrays for each granule # ------------------------------------------------ @@ -244,7 +244,7 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, if type(Path) is list: if len(Path) == 0: self.nobs = 0 - print "WARNING: Empty Vx04_L2 object created" + print("WARNING: Empty Vx04_L2 object created") return else: Path = [Path, ] @@ -254,7 +254,7 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, # -------------------------------- if len(self.Scattering_Angle) == 0: self.nobs = 0 - print "WARNING: Empty MxD04_L2 object created" + print("WARNING: Empty MxD04_L2 object created") return # Make each attribute a single numpy array @@ -265,7 +265,7 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, try: self.__dict__[sds] = np.ma.concatenate(self.__dict__[sds]) except: - print "Failed concatenating "+sds + print("Failed concatenating "+sds) # Determine index of "good" observations @@ -279,7 +279,7 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, elif self.algo == 'DB_OCEAN': self.iGood = self.Aerosol_Optical_Thickness_QA_Flag_Ocean == BEST else: - raise ValueError, 'invalid algorithm (very strange)' + raise ValueError('invalid algorithm (very strange)') # Keep only "good" observations @@ -293,13 +293,13 @@ def __init__ (self,Path,algo,syn_time=None,nsyn=8,Verb=0, elif rank == 2: self.__dict__[sds] = self.__dict__[sds][m,:] else: - raise IndexError, 'invalid rank=%d'%rank + raise IndexError('invalid rank=%d'%rank) self.iGood = self.iGood[m] # Make aliases for compatibility with older code # ---------------------------------------------- - Alias = self.ALIAS.keys() + Alias = list(self.ALIAS.keys()) for sds in self.SDS: if sds in Alias: self.__dict__[self.ALIAS[sds]] = self.__dict__[sds] @@ -494,7 +494,7 @@ def _readList(self,List): else: self._readGranuleDT(item) else: - print "%s is not a valid file or directory, ignoring it"%item + print("%s is not a valid file or directory, ignoring it"%item) #--- def _readDir(self,dir): """Recursively, look for files in directory.""" @@ -503,7 +503,7 @@ def _readDir(self,dir): if os.path.isdir(path): self._readDir(path) elif os.path.isfile(path): self._readGranule(path) else: - print "%s is not a valid file or directory, ignoring it"%item + print("%s is not a valid file or directory, ignoring it"%item) #--- def _readGranuleDB(self,filename): @@ -513,11 +513,11 @@ def _readGranuleDB(self,filename): # --------------------------------------- try: if self.verb: - print "[] Working on "+filename + print("[] Working on "+filename) nc = Dataset(filename) except: if self.verb > 2: - print "- %s: not recognized as an netCDF file"%filename + print("- %s: not recognized as an netCDF file"%filename) return # Read select variables (reshape to allow concatenation later) @@ -542,7 +542,7 @@ def _readGranuleDB(self,filename): elif len(v.shape) == 2: v = v.ravel() else: - raise IndexError, "invalid shape for SDS <%s>"%sds + raise IndexError("invalid shape for SDS <%s>"%sds) self.__dict__[sds].append(v) @@ -568,13 +568,13 @@ def _readGranuleDT(self,filename): # --------------------------------------- try: if self.verb: - print "[] Working on "+filename + print("[] Working on "+filename) nc = Dataset(filename) data = nc.groups['geophysical_data'] loc = nc.groups['geolocation_data'] except: if self.verb > 2: - print "- %s: not recognized as an netCDF file"%filename + print("- %s: not recognized as an netCDF file"%filename) return # Read select variables (reshape to allow concatenation later) @@ -605,7 +605,7 @@ def _readGranuleDT(self,filename): elif len(v.shape) == 2: v = v.ravel() else: - raise IndexError, "invalid shape for SDS <%s>"%sds + raise IndexError("invalid shape for SDS <%s>"%sds) self.__dict__[sds].append(v) @@ -631,7 +631,7 @@ def reduce(self,I): """ Reduce observations according to index I. """ - Nicknames = self.ALIAS.values() + Nicknames = list(self.ALIAS.values()) for name in self.__dict__: if name in Nicknames: continue # alias do not get reduced @@ -641,7 +641,7 @@ def reduce(self,I): # print "{} Reducing "+name self.__dict__[name] = q[I] - Alias = self.ALIAS.keys() + Alias = list(self.ALIAS.keys()) for sds in self.SDS: if sds in Alias: self.__dict__[self.ALIAS[sds]] = self.__dict__[sds] # redefine aliases @@ -691,7 +691,7 @@ def write(self,filename=None,dir='.',expid=None,Verb=1): reflectance = self.reflectance) if Verb >=1: - print "[w] Wrote file "+filename + print("[w] Wrote file "+filename) #--- def writeODS(self,filename=None,dir='.',expid=None,channels=None, @@ -702,7 +702,7 @@ def writeODS(self,filename=None,dir='.',expid=None,channels=None, """ if self.syn_time == None: - raise ValuError, "synoptic time missing, cannot write ODS" + raise ValuError("synoptic time missing, cannot write ODS") # Stop here if no good obs available # ---------------------------------- @@ -728,7 +728,7 @@ def writeODS(self,filename=None,dir='.',expid=None,channels=None, i = 0 ks = np.arange(ns) + 1 for ch in channels: - I = range(i,i+ns) + I = list(range(i,i+ns)) j = list(channels).index(ch) # index of channel ods.ks[I] = ks ods.lat[I] = self.lat[:] @@ -761,7 +761,7 @@ def writeODS(self,filename=None,dir='.',expid=None,channels=None, ods_ = ods.select(qcx=0) if Verb >=1: - print "[w] Writing file <"+filename+"> with %d observations"%ods_.nobs + print("[w] Writing file <"+filename+"> with %d observations"%ods_.nobs) ods_.write(filename,self.nymd,self.nhms,nsyn=8,ftype='pre_anal') @@ -905,7 +905,7 @@ def writeg(self,filename=None,dir='.',expid=None,refine=8,res=None, # pass if Verb >=1: - print "[w] Wrote file "+filename + print("[w] Wrote file "+filename) #--- def addVar(self,ga,expr='mag(u10m,v10m)',vname='wind',clmYear=None,tight=True): @@ -1018,13 +1018,13 @@ def granules ( path, algo, sat, syn_time, coll='011', nsyn=8, verbose=False ): filen = glob(basen)[0] Granules += [filen,] if verbose: - print " [x] Found "+filen + print(" [x] Found "+filen) except: pass t += dt if len(Granules) == 0: - print "WARNING: no %s collection %s granules found for"%(algo,coll), syn_time + print("WARNING: no %s collection %s granules found for"%(algo,coll), syn_time) return Granules @@ -1037,25 +1037,25 @@ def print_stats(name,x=None): x = name name = 'mean,stdv,rms,min,25%,median,75%,max: ' if name == '__header__': - print '' + print('') n = (80 - len(x))/2 - print n * ' ' + x - print n * ' ' + len(x) * '-' - print '' - print ' Name mean stdv rms min 25% median 75% max' - print ' --------- ------- ------- ------- ------- ------- ------- ------- -------' + print(n * ' ' + x) + print(n * ' ' + len(x) * '-') + print('') + print(' Name mean stdv rms min 25% median 75% max') + print(' --------- ------- ------- ------- ------- ------- ------- ------- -------') elif name == '__sep__': - print ' --------- ------- ------- ------- ------- ------- ------- ------- -------' + print(' --------- ------- ------- ------- ------- ------- ------- ------- -------') elif name == '__footer__': - print ' --------- ------- ------- ------- ------- ------- ------- ------- -------' - print '' + print(' --------- ------- ------- ------- ------- ------- ------- ------- -------') + print('') else: ave = x.mean() std = x.std() rms = np.sqrt(ave*ave+std*std) prc = prctile(x) - print '%10s %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f '%\ - (name,ave,std,rms,prc[0],prc[1],prc[2],prc[3],prc[4]) + print('%10s %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f '%\ + (name,ave,std,rms,prc[0],prc[1],prc[2],prc[3],prc[4])) #-- From 02358cfd3af045a9182864bf6b0271b55aa10617 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 25 May 2023 15:20:43 -0400 Subject: [PATCH 60/67] PC- fix in giant_viirs.py to deal with reading datetimes character arrays from netcdf. need to convert to unicode characters. --- src/Components/misc/obs_aod/ABC/giant_viirs.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/Components/misc/obs_aod/ABC/giant_viirs.py b/src/Components/misc/obs_aod/ABC/giant_viirs.py index 16b837c1..2c906dcc 100644 --- a/src/Components/misc/obs_aod/ABC/giant_viirs.py +++ b/src/Components/misc/obs_aod/ABC/giant_viirs.py @@ -293,7 +293,10 @@ def __init__ (self,filename,xVars=(),only_good=True,tymemax=None): # new files use masked arrays # convert everythong to regular array filling with -9999.0 # make sure _fill_value is -9999.0 - self.__dict__[name] = np.array(data) + if data.dtype == np.dtype('S1'): + self.__dict__[name] = np.array(data).astype(str) + else: + self.__dict__[name] = np.array(data) self.giantList.append(name) nc.close() @@ -303,7 +306,7 @@ def __init__ (self,filename,xVars=(),only_good=True,tymemax=None): nc = Dataset(filename) if 'ISO_DateTime' in list(nc.variables.keys()): try: - iso = nc.variables['ISO_DateTime'][:] + iso = nc.variables['ISO_DateTime'][:].astype(str) self.tyme = array([isoparse(''.join(array(t))) for t in iso]) except: # old file only have Date and Time variables From 08ee5d1fc9960ee2bf85c478d42e704a48ef0d28 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 25 May 2023 15:22:23 -0400 Subject: [PATCH 61/67] PC - update GMAOpyobs to 1.0.3 which contains the vx04.py VIIRS reader --- components.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components.yaml b/components.yaml index 5878e035..2392c54f 100644 --- a/components.yaml +++ b/components.yaml @@ -28,5 +28,5 @@ MAPL: GMAOpyobs: local: ./src/Shared/GMAO_Shared/GMAO_pyobs@ remote: ../GMAOpyobs.git - tag: v1.0.2 + tag: v1.0.3 develop: develop From 79d34b78157a487bef2952d6bd1af493611fc647 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 25 May 2023 16:23:39 -0400 Subject: [PATCH 62/67] PC - change #! to python3 in viirs_l2a and vx04_l2a in GAAS_App --- src/Applications/GAAS_App/viirs_l2a.py | 4 ++-- src/Applications/GAAS_App/vx04_l2a.py | 6 ++---- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/Applications/GAAS_App/viirs_l2a.py b/src/Applications/GAAS_App/viirs_l2a.py index c9babc66..50c63376 100755 --- a/src/Applications/GAAS_App/viirs_l2a.py +++ b/src/Applications/GAAS_App/viirs_l2a.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 # -W ignore::DeprecationWarning """ @@ -15,7 +15,7 @@ from os import system from optparse import OptionParser -from MAPL import Config +from MAPL.config import Config if __name__ == "__main__": diff --git a/src/Applications/GAAS_App/vx04_l2a.py b/src/Applications/GAAS_App/vx04_l2a.py index 0cb0a7d5..8af5ce67 100755 --- a/src/Applications/GAAS_App/vx04_l2a.py +++ b/src/Applications/GAAS_App/vx04_l2a.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 # -W ignore::DeprecationWarning """ @@ -21,11 +21,10 @@ import sys import subprocess -from time import clock from optparse import OptionParser # Command-line args from dateutil.parser import parse as isoparse from vx04_nnr import Vx04_NNR -from MAPL import strTemplate +from MAPL.config import strTemplate Ident = dict( vsnppdto = ('SNPP','dt_ocean'), vsnppdtl = ('SNPP','dt_land'), @@ -141,7 +140,6 @@ def makethis_dir(filename): print(" VIIRS Level 2A Processing") print(" -------------------------") print("") - t0 = clock() # Time variables # -------------- From 1d97db4ca52f6f298b321598f4f5d19164145454 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 25 May 2023 16:25:05 -0400 Subject: [PATCH 63/67] PC - fix import in pyods __ini__.py to use relative imports. needed for python3 --- src/Shared/GMAO_Shared/GMAO_ods/pyods/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) mode change 100644 => 100755 src/Shared/GMAO_Shared/GMAO_ods/pyods/__init__.py diff --git a/src/Shared/GMAO_Shared/GMAO_ods/pyods/__init__.py b/src/Shared/GMAO_Shared/GMAO_ods/pyods/__init__.py old mode 100644 new mode 100755 index 3475606a..6f3ef655 --- a/src/Shared/GMAO_Shared/GMAO_ods/pyods/__init__.py +++ b/src/Shared/GMAO_Shared/GMAO_ods/pyods/__init__.py @@ -5,8 +5,8 @@ """ from types import * -from pyods_ import * # Fortran extension using f2py -from odsmeta import * # useful constants +from .pyods_ import * # Fortran extension using f2py +from .odsmeta import * # useful constants import numpy as np __VERSION__ = '1.0.3' From 95c8a29d17dcf6d760e6d9d1d05e986cb508200a Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Thu, 25 May 2023 16:32:26 -0400 Subject: [PATCH 64/67] PC - update changelog for VIIRS updates and for v2.0.2 release --- CHANGELOG.md | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8bfb6042..f91c3201 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,22 @@ ### Changed +## [v2.0.2] - 2023-05-25 + +### Added + +- VIIRS NNR training code - new VIIRS giant reader +- New code in GAAS_App to generate NNR ODS files from VIIRS obs + +### Fixed + +- fix in NNR code for new syntax of the Kfold generator in sklearn +- fix in NNR testing code to protect against cases where VIIRS standard product may not retrieve multiple wavelengths +- fix pyods __ini__.py to use relative imports. needed for python3 +### Changed + +- update GMAOpyobs to v1.0.3 + ## [v2.0.1] - 2023-05-17 ### Added From 578c82c17e902fa483dcdd974b7abca7c9001fe7 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 5 Jun 2023 15:41:00 -0400 Subject: [PATCH 65/67] PC - add VIIRS aerosol products to the GMAO_ods kx list --- src/Shared/GMAO_Shared/GMAO_ods/kx_list.rc | 8 ++++++++ src/Shared/GMAO_Shared/GMAO_ods/pyods_.F90 | 8 ++++++++ 2 files changed, 16 insertions(+) diff --git a/src/Shared/GMAO_Shared/GMAO_ods/kx_list.rc b/src/Shared/GMAO_Shared/GMAO_ods/kx_list.rc index 19a657ff..d20e518d 100644 --- a/src/Shared/GMAO_Shared/GMAO_ods/kx_list.rc +++ b/src/Shared/GMAO_Shared/GMAO_ods/kx_list.rc @@ -341,6 +341,14 @@ kx_list:: 330 GOES-17 Aerosol (Dark Target Land Algorithm) 331 Himawari Aerosol (Dark Target Ocean Algorithm) 332 Himawari Aerosol (Dark Target Land Algorithm) + 333 VIIRS SNPP Aerosol (Deep Blue Land Algorithm) + 334 VIIRS SNPP Aerosol (Deep Blue Ocean Algorithm) + 335 VIIRS SNPP Aerosol (Dark Target Land Algorithm) + 336 VIIRS SNPP Aerosol (Dark Target Ocean Algorithm) + 337 VIIRS NOAA-20 Aerosol (Deep Blue Land Algorithm) + 338 VIIRS NOAA-20 Aerosol (Deep Blue Ocean Algorithm) + 339 VIIRS NOAA-20 Aerosol (Dark Target Land Algorithm) + 340 VIIRS NOAA-20 Aerosol (Dark Target Ocean Algorithm) :: # Notes: diff --git a/src/Shared/GMAO_Shared/GMAO_ods/pyods_.F90 b/src/Shared/GMAO_Shared/GMAO_ods/pyods_.F90 index f4051132..e24d66b8 100644 --- a/src/Shared/GMAO_Shared/GMAO_ods/pyods_.F90 +++ b/src/Shared/GMAO_Shared/GMAO_ods/pyods_.F90 @@ -630,6 +630,14 @@ subroutine pyods_putAll(filename, ftype, nymd,nhms, nsyn, nobs, & ods%meta%kx_names(330) = 'GOES-17 Aerosol (Dark Target Land Algorithm)' ods%meta%kx_names(331) = 'Himawari Aerosol (Dark Target Ocean Algorithm)' ods%meta%kx_names(332) = 'Himawari Aerosol (Dark Target Land Algorithm)' + ods%meta%kx_names(333) = 'VIIRS SNPP Aerosol (Deep Blue Land Algorithm)' + ods%meta%kx_names(334) = 'VIIRS SNPP Aerosol (Deep Blue Ocean Algorithm)' + ods%meta%kx_names(335) = 'VIIRS SNPP Aerosol (Dark Target Land Algorithm)' + ods%meta%kx_names(336) = 'VIIRS SNPP Aerosol (Dark Target Ocean Algorithm)' + ods%meta%kx_names(337) = 'VIIRS NOAA-20 Aerosol (Deep Blue Land Algorithm)' + ods%meta%kx_names(338) = 'VIIRS NOAA-20 Aerosol (Deep Blue Ocean Algorithm)' + ods%meta%kx_names(339) = 'VIIRS NOAA-20 Aerosol (Dark Target Land Algorithm)' + ods%meta%kx_names(340) = 'VIIRS NOAA-20 Aerosol (Dark Target Ocean Algorithm)' call ods_put (filename, ftype, nymd, nhms, ods, rc) if ( rc .ne. 0 ) return From db8d19f87adc63e85df4d4b4a0f00aaa1322894c Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Mon, 5 Jun 2023 15:42:19 -0400 Subject: [PATCH 66/67] update changelog with VIIRS kx addition to GMAO_ods --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index f91c3201..ef06e717 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,7 @@ - fix pyods __ini__.py to use relative imports. needed for python3 ### Changed +- add VIIRS aerosol products to GMAO_ods kx list - update GMAOpyobs to v1.0.3 ## [v2.0.1] - 2023-05-17 From 4fb8453b88f7a069c3c85779fd0bb67129e00fda Mon Sep 17 00:00:00 2001 From: Patricia Castellanos <patricia.castellanos@nasa.gov> Date: Wed, 7 Jun 2023 15:50:08 -0400 Subject: [PATCH 67/67] update changelog for v2.0.2 release. includes all the updates to do viirs nnr training. this will migrate over to AeroML. --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ef06e717..327bfa41 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -23,7 +23,7 @@ ### Changed - add VIIRS aerosol products to GMAO_ods kx list -- update GMAOpyobs to v1.0.3 +- update GMAOpyobs to v1.0.4 ## [v2.0.1] - 2023-05-17