diff --git a/components.yaml b/components.yaml index 8b8b0461..550245b4 100644 --- a/components.yaml +++ b/components.yaml @@ -48,7 +48,7 @@ GMAOpyobs: AeroML: local: ./src/Shared/@AeroML remote: ../AeroML.git - tag: rt1.0.0-cak_tag0 + tag: v1.1.0 develop: develop sparse: ./config/AeroML.sparse diff --git a/src/Applications/GAAS_App/CMakeLists.txt b/src/Applications/GAAS_App/CMakeLists.txt index 79ef8578..bd956b1b 100644 --- a/src/Applications/GAAS_App/CMakeLists.txt +++ b/src/Applications/GAAS_App/CMakeLists.txt @@ -17,24 +17,18 @@ esma_add_library(${this} ecbuild_add_executable(TARGET ana_aod.x SOURCES ana_aod.F LIBS ${this} GMAO_gfio_r4 Chem_Base Chem_Shared) ecbuild_add_executable(TARGET mpana_aod.x SOURCES mpana_aod.F90 LIBS ${this} GMAO_gfio_r4 Chem_Base Chem_Shared) -set(pythonscripts modis_l2a.py mxd04_l2a.py avhrr_l2a.py patmosx_l2a.py aod_data.py) +set(pythonscripts aod_data.py) install(PROGRAMS ${pythonscripts} get_aero_obs.csh run_gaas_ana.csh DESTINATION bin) -install( - FILES mxd04_nnr.py avhrr_nnr.py - DESTINATION lib/Python - ) - if (EXTENDED_SOURCE) set_target_properties (${this} PROPERTIES COMPILE_FLAGS ${EXTENDED_SOURCE}) set_target_properties (ana_aod.x PROPERTIES COMPILE_FLAGS ${EXTENDED_SOURCE}) endif() file(GLOB rc_files *.rc) -file(GLOB pcf_files *.pcf) file(GLOB tmpl_files *.tmpl) install ( - FILES ${rc_files} ${pcf_files} ${tmpl_files} + FILES ${rc_files} ${tmpl_files} DESTINATION etc ) diff --git a/src/Applications/GAAS_App/aeronet_all.py b/src/Applications/GAAS_App/aeronet_all.py deleted file mode 100755 index 23a39ea6..00000000 --- a/src/Applications/GAAS_App/aeronet_all.py +++ /dev/null @@ -1,86 +0,0 @@ -#!/usr/bin/env python3 -""" -Splits AERONET into synoptic chunks. -""" - -import os -import sys - -from datetime import date, datetime, timedelta - -from pyobs.aeronet import AERONET_L2, granules - - -#--------------------------------------------------------------------- -def makethis_dir(path): - """Creates the relevant directory if necessary.""" - if path != '': - rc = os.system('mkdir -p '+path) - if rc: - raise IOError("could not create directory "+path) - -if __name__ == "__main__": - - RootDir = '/nobackup/AERONET/MERRA-2' - ods_blank = '/nobackup/NNR/Misc/blank.ods' - - oneday = timedelta(seconds=24*60*60) - onehour = timedelta(seconds=60*60) - - if len(sys.argv)<2: - print("Usage:") - print(" aeronet4_all.py year1 [year2]") - sys.exit(1) - else: - y1 = sys.argv[1] - if len(sys.argv)>2: - y2 = sys.argv[2] - else: - y2 = y1 - Years = list(range(int(y1),int(y2)+1)) - - # Loop over years - # --------------- - for year in Years: - tyme0 = datetime(year,1,1) - toy = datetime(year,12,31) - tyme0 - ndays = 1 + int(toy.total_seconds()/(24*60*60)) - for doy in range(1,ndays+1): - - today = tyme0 + (doy-1) * oneday - - print('Day: ', today) - - # Read AERONET for this day - # ------------------------- - Files = granules(today,bracket='left') - a = AERONET_L2(Files,Verbose=True) - - # Output directories - # ------------------ - dirL2 = RootDir + '/Level2/ODS/Y%d/M%02d'%(today.year,today.month) - dirL3 = RootDir + '/Level3/Y%d/M%02d'%(today.year,today.month) - - for h in range(0,24,3): - - tyme = tyme0 + (doy-1) * oneday + h * onehour - - # Make sure directories exist - # --------------------------- - makethis_dir(dirL2) - makethis_dir(dirL3) - - # Write & compress the files - # -------------------------- - fnL2, nobs = a.writeODS(tyme,dir=dirL2) - fnL3 = a.writeGridded(tyme,dir=dirL3) - - # Compress daily file - # ------------------- - if a.nobs>0: - if os.system("n4zip "+fnL2+" > /dev/null"): - warnings.warn('cannot compress output ODS file <%s>'%fnL2) - if os.system("n4zip "+fnL3+" > /dev/null"): - warnings.warn('cannot compress output NC4 file <%s>'%fnL3) - else: - os.system("touch %s.empty"%fnL2) diff --git a/src/Applications/GAAS_App/avhrr_all.py b/src/Applications/GAAS_App/avhrr_all.py deleted file mode 100755 index 4085ba40..00000000 --- a/src/Applications/GAAS_App/avhrr_all.py +++ /dev/null @@ -1,46 +0,0 @@ -#!/usr/bin/env python3 -""" -Splits AVHRR into synoptic chunks. -""" - -import os -import sys - -from datetime import date, datetime, timedelta - -if __name__ == "__main__": - - oneday = timedelta(seconds=24*60*60) - onehour = timedelta(seconds=60*60) - - if len(sys.argv)<2: - print("Usage:") - print(" avhrr_all.py year1 [year2]") - sys.exit(1) - else: - y1 = sys.argv[1] - if len(sys.argv)>2: - y2 = sys.argv[2] - else: - y2 = y1 - Years = list(range(int(y1),int(y2)+1)) - - # Loop over years - # --------------- - for year in Years: - tyme0 = datetime(year,1,1) - toy = datetime(year,12,31) - tyme0 - ndays = 1 + int(toy.total_seconds()/(24*60*60)) - for doy in range(1,ndays+1): - - for h in range(0,24,3): - - tyme = tyme0 + (doy-1) * oneday + h * onehour - - nymd = tyme.year*10000 + tyme.month*100 + tyme.day - nhms = tyme.hour*10000 + tyme.minute*100 + tyme.second - - cmd = 'python avhrr_l2a.py -v asc %d %d'%(nymd,nhms) - - print(cmd) - os.system(cmd) diff --git a/src/Applications/GAAS_App/avhrr_l2a.pcf b/src/Applications/GAAS_App/avhrr_l2a.pcf deleted file mode 100644 index 867064e8..00000000 --- a/src/Applications/GAAS_App/avhrr_l2a.pcf +++ /dev/null @@ -1,25 +0,0 @@ -#AVHRR_L2A processing - -AVHRR_L2A_VERBOSE = YES - -AVHRR_L2A_EXPID = ${EXPID} -AVHRR_L2A_ORBITS = asc,des - -#--AVHRR_L2A_L2_DIR = /archive/input/dao_ops/obs/reanalysis/patmosx/Synoptic -#--AVHRR_L2A_OUT_DIR = ./Results/obs/Level%lev/%prod/Y%y4/M%m2 -AVHRR_L2A_L2_DIR = ${FVWORK} -AVHRR_L2A_OUT_DIR = ${FVWORK} - -AVHRR_L2A_OVERWRITE = YES -AVHRR_L2A_OUT_TEMPLATE = '%s.%prod_L%leva.%orb.%y4%m2%d2_%h2%n2z.%ext' -AVHRR_L2A_RESOLUTION = e - -AVHRR_L2A_WIND_FILE = ${EXPID}.gaas_bkg.sfc.%y4%m2%d2_%h2z.nc4 -AVHRR_L2A_TPW_FILE = ${EXPID}.gaas_bkg.sfc.%y4%m2%d2_%h2z.nc4 -AVHRR_L2A_AOD_FILE = ${EXPID}.gaas_bkg.sfc.%y4%m2%d2_%h2z.nc4 - -AVHRR_L2A_NN_FILE = ExtData/g5chem/x/NN/nnr_001.avhrr_Tau.net -AVHRR_L2A_BLANK_ODS = ExtData/g5chem/x/blank_syn8.ods -AVHRR_L2A_COXMUNK_LUT = ExtData/g5chem/x/avhrr.cox-munk_lut.npz - -#END AVHRR_L2A processing diff --git a/src/Applications/GAAS_App/avhrr_l2a.py b/src/Applications/GAAS_App/avhrr_l2a.py deleted file mode 100755 index 9ffbf0c3..00000000 --- a/src/Applications/GAAS_App/avhrr_l2a.py +++ /dev/null @@ -1,64 +0,0 @@ -#!/usr/bin/env python3 -# -W ignore::DeprecationWarning - -""" - - Simple wrapper script to parse Prep config file and create ODS with NNR AVHRR retrievals. - - October 2013 - arlindo.dasilva@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 nymd nhms", - version='avhrr_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) == 3: - prep_config, nymd, nhms = args - else: - parser.error("must have 3 arguments: prep_config_filename nymd nhms") - - # Parse prep config - # ----------------- - cf = Config(prep_config,delim=' = ') - - Options = " --expid=" + cf('AVHRR_L2A_EXPID') + \ - " --l2_dir=" + cf('AVHRR_L2A_L2_DIR') + \ - " --res=" + cf('AVHRR_L2A_RESOLUTION') + \ - " --dir=" + cf('AVHRR_L2A_OUT_DIR') + \ - " --fname=" + cf('AVHRR_L2A_OUT_TEMPLATE') + \ - " --net=" + cf('AVHRR_L2A_NN_FILE') + \ - " --wind=" + cf('AVHRR_L2A_WIND_FILE') + \ - " --tpw=" + cf('AVHRR_L2A_TPW_FILE') + \ - " --aod=" + cf('AVHRR_L2A_AOD_FILE') + \ - " --blank_ods=" + cf('AVHRR_L2A_BLANK_ODS') + \ - " --coxmunk=" + cf('AVHRR_L2A_COXMUNK_LUT') - - if cf('AVHRR_L2A_OVERWRITE').upper() == 'YES': Options += " --force" - if cf('AVHRR_L2A_VERBOSE').upper() == 'YES': Options += " -v" - - # Generate products - # ----------------- - i = 0 - for orbit in cf('AVHRR_L2A_ORBITS').split(','): - cmd = "patmosx_l2a.py %s %s %s %s"%(Options,orbit,nymd,nhms) - print(cmd) - if not options.dryrun: - if system(cmd): - raise ValueError("patmosx_l2a.py failed for %s on %s %s"%(orbit,nymd,nhms)) - - i += 1 - - diff --git a/src/Applications/GAAS_App/avhrr_nnr.py b/src/Applications/GAAS_App/avhrr_nnr.py deleted file mode 100644 index 0266cadc..00000000 --- a/src/Applications/GAAS_App/avhrr_nnr.py +++ /dev/null @@ -1,153 +0,0 @@ -""" -This module implements the evaluation of AVHRR Neural Net Retrievals (NNR) -based on PATMOS-X data. - -""" - -import warnings -from pyobs import avhrr, sknet -from numpy import c_ as cat -from numpy import copy, ones, sin, cos, exp, arccos, pi, any, log - -MISSING = -1.0e20 -d2r = pi / 180. - -class AVHRR_NNR(avhrr.AVHRR_L2B): - """ - This class extends the AVHRR Level 2B class by adding NNR evaluation - methods. - """ - - def __init__(self, Path='/nobackup/AVHRR/Level2/NPZ/2008/*_00?.npz', - N_bal=None,Verb=False): - - self.verbose = Verb - self.ident = 'avhrr' - self.surface = 'ocean' - - avhrr.AVHRR_L2B.__init__(self,Path,Verb=Verb) - - # Balance Observing system - # ------------------------ - if N_bal is not None: - I = self.balance(N_bal) - self.reduce(I) - - # Air mass factor - # --------------- - self.amf = (1./cos(d2r*self.SolarZenith))+(1./cos(d2r*self.SensorZenith)) - - # Glint Angle - # ----------- - RelativeAzimuth = self.SensorAzimuth # = anchor_relative_azimuth - cosGlintAngle = cos(self.SolarZenith*d2r) * cos(self.SensorZenith*d2r) + \ - sin(self.SolarZenith*d2r) * sin(self.SensorZenith*d2r) * \ - cos(RelativeAzimuth*d2r) - - # Angle transforms: for NN calculations we work with cosine of angles - # ------------------------------------------------------------------- - self.SensorAzimuth = cos(self.SensorAzimuth*d2r) - self.SensorZenith = cos(self.SensorZenith*d2r) - self.SolarAzimuth = cos(self.SolarAzimuth*d2r) - self.SolarZenith = cos(self.SolarZenith*d2r) - self.GlintAngle = cosGlintAngle - - # Sanity check - # ------------ - self.iValid = (self.tau_630 > -0.01) &\ - (self.ref_630 > 0) &\ - (self.ref_860 > 0) -# (self.tau_860 > -0.01) &\ - - # Log transforms - # -------------- - self.ltau_630 = log(self.tau_630+0.01) -# self.ltau_860 = log(self.tau_860+0.01) - self.lref_630 = log(self.ref_630) - self.lref_860 = log(self.ref_860) - - def _getInputs(self): - """ - Get Inputs for Neural Net. - """ - - # Loop over inputs - # ---------------- - first = True - for inputName in self.net.InputNames: - - if self.verbose>0: - print('Getting NN input ',inputName) - - # Retrieve input - # -------------- - input = self.__dict__[inputName][:] - self.iValid = self.iValid & (input!=MISSING) # Q/C - - # Concatenate Inputs - # ------------------ - if first: - inputs = input - first = False - else: - inputs = cat[inputs,input] - - # Keep only good observations - # --------------------------- - return inputs[self.iValid,:] - -#-- - def apply(self,nnFile='/nobackup/NNR/Net/nnr_001b.avhrr_Tau.net'): - """ - Evaluates NN retrieval. - """ - - # Load Network - # ------------ - self.net = sknet.loadnet(nnFile) - - # Stop here is no good obs available - # ---------------------------------- - if self.nobs == 0: - return # no data to work with - if any(self.iValid) == False: - return # no good data to work with - - if len(self.net.TargetNames)>1: - raise ValueError('Strange, more than one predictor') - - # Evaluate NN on inputs - # --------------------- - targets = self.net(self._getInputs()) - - name = self.net.TargetNames[0] - if self.verbose>0: - print("Evaluating NNR for <%s> with Log-AOD = "%name, self.net.laod) - - # Output is always AOD - # -------------------- - if self.net.laod: - result = exp(targets) - 0.01 # inverse - - else: - result = targets - - # Set retrieved values, possibly with UNDEFS - # ------------------------------------------ - self.__dict__[name] = MISSING * ones(self.lon.shape) - self.channels_ = [550.,] # channels being retrieved - self.__dict__[name][self.iValid] = result.ravel() - - return result - -#--- - - __call__= apply - -#--- - -if __name__ == "__main__": - - a = AVHRR_NNR(Verb=True) - - aod = a.apply() diff --git a/src/Applications/GAAS_App/misr_land.py b/src/Applications/GAAS_App/misr_land.py deleted file mode 100755 index f41311de..00000000 --- a/src/Applications/GAAS_App/misr_land.py +++ /dev/null @@ -1,139 +0,0 @@ -#!/usr/bin/env python3 -""" - -Utility to read MISR ODS and write a smaller file with data only over land. - -""" - -import sys -import os - -from numpy import sort, array - -from pyods import ODS -from pyobs import igbp - -Bright = True # if True, keep only pixels brighter than Alb_min -Alb_min = 0.15 - -vpath = '/nobackup/Emissions/Vegetation/GL_IGBP_INPE/' -apath = '/nobackup/MODIS/Level3/ALBEDO/albedo_clim.ctl' - -if __name__ == "__main__": - - def blank(nymd,nhms): - print(" - NO DATA for %8d %6d "%(nymd, nhms)) - ods = ODS(nobs=1,kt=45,kx=313) - ods.qcx[0] = 1 - ods.ks[0] = 1 - ods.lev[0] = 550 - return ods - - def doList(List): - """ - Recursively, look for files in list; list items can - be files or directories. - """ - for item in List: - if os.path.isdir(item): doDir(item) - elif os.path.isfile(item): doFile(item) - else: - print("%s is not a valid file or directory, ignoring it"%item) - - def doDir(dir): - """Recursively, look for files in directory.""" - for item in sort(os.listdir(dir)): - path = dir + os.sep + item - if os.path.isdir(path): doDir(path) - elif os.path.isfile(path): doFile(path) - else: - print("%s is not a valid file or directory, ignoring it"%item) - - - def doFile(fname): - - if fname.split('.')[-1] != 'ods': - print("[*] NOT ODS "+fname) - return - - dirn = os.path.abspath(os.path.dirname(fname)) - basen = os.path.basename(fname) - - nymd = int(fname.split('.')[-2]) - - if Bright: - landdn = dirn.replace('/ODS/','/ODS_Bright/') - landfn = landdn+'/'+basen.replace('aero_','bright_') - else: - landdn = dirn.replace('/ODS/','/ODS_Land/') - landfn = landdn+'/'+basen.replace('aero_','aland_') - - if os.path.exists(landfn): - print("[ ] Skipping "+landfn) - return - - os.system('/bin/mkdir -p '+landdn) - - print("[x] Working on "+landfn) - - for nhms in range(0,240000,30000): - - # AOD and 550 nm - # -------------- - ods = ODS(fname,nymd,nhms,only_good=True).select(kt=45,lev=558.) - - # Get vegetation type as means of filtering out water - # --------------------------------------------------- - try: - v = igbp.getDetailedVeg(ods.lon,ods.lat,vpath) - - # I = (v==15)|(v>=17)|(ods.lat<-60.) # water and ice - I = (v==15)|(v>=17) # water and ice - I = array((I==False)) # land - - if any(I): - ods = ods.__copy__(Indices=I) - else: - ods = blank(nymd,nhms) - - if Bright and any(I): - ods.addVar(ga,expr='albedo',clmYear=2000) - I = (ods.albedo>Alb_min) - if any(I): - ods = ods.__copy__(Indices=I) - else: - ods = blank(nymd,nhms) - - if any(I): - print(" + Processing %8d %6d with %5d obs left"%(nymd, nhms, len(v[I]))) - - except: - ods = blank(nymd,nhms) - - if ods.nobs == 0: - ods = blank(nymd,nhms) - - ods.write(landfn,nymd,nhms,nsyn=8) - - - os.system("n4zip "+landfn) - - #--- - - if len(sys.argv) < 2 : - print('Usage: ') - print(' %s misr_ods_files/dir_names'%os.path.basename(sys.argv[0])) - sys.exit(1) - - # Albedo dataset - # -------------- - if Bright: - from grads import GrADS - ga = GrADS(Echo=False,Window=False) - fh = ga.open(apath) - print(fh) - - # Process list of files or directories - # ------------------------------------ - doList(sys.argv[1:]) - diff --git a/src/Applications/GAAS_App/modis_l2a.pcf b/src/Applications/GAAS_App/modis_l2a.pcf deleted file mode 100644 index 8c11dee5..00000000 --- a/src/Applications/GAAS_App/modis_l2a.pcf +++ /dev/null @@ -1,25 +0,0 @@ -# MODIS_L2A processing: NNR v3 for MODIS C6 -# -# GEOS-5 FP -# -#.................................................................................. - -MODIS_L2A_VERBOSE = YES - -MODIS_L2A_EXPID = ${EXPID} -MODIS_L2A_IDENTS = mydl,mydo,mydd,modl,modo,modd -MODIS_L2A_COLLECTION = 061,061,061,061,061,061 - -MODIS_L2A_L2_DIR = ${FVWORK} -MODIS_L2A_OUT_DIR = ${FVWORK} - -MODIS_L2A_OVERWRITE = YES -MODIS_L2A_OUT_TEMPLATE = '%s.%prod_L%leva.%algo.%y4%m2%d2_%h2%n2z.%ext' -MODIS_L2A_RESOLUTION = e - -MODIS_L2A_AER_X = %s.gaas_bkg.sfc.%y4%m2%d2_%h2z.nc4 - -MODIS_L2A_NN_FILE = ExtData/g5chem/x/NN/nnr_003.%ident_Tau.net -MODIS_L2A_BLANK_ODS = ExtData/g5chem/x/blank_syn8.ods - -#END MODIS_L2A processing diff --git a/src/Applications/GAAS_App/modis_l2a.py b/src/Applications/GAAS_App/modis_l2a.py deleted file mode 100755 index 440bd65b..00000000 --- a/src/Applications/GAAS_App/modis_l2a.py +++ /dev/null @@ -1,63 +0,0 @@ -#!/usr/bin/env python3 -# -W ignore::DeprecationWarning - -""" - - Simple wrapper script to parse Prep config file and create ODS with MODIS NNR aerosol retrievals. - - February 2011. - arlindo.dasilva@nasa.gov -""" - -from os import system -from optparse import OptionParser -from MAPL.config import Config - -if __name__ == "__main__": - -# Parse command line options -# -------------------------- - parser = OptionParser(usage="Usage: %prog prep_config_file isotime", - version='modis_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('MODIS_L2A_EXPID') + \ - " --l2_dir=" + cf('MODIS_L2A_L2_DIR') + \ - " --res=" + cf('MODIS_L2A_RESOLUTION') + \ - " --dir=" + cf('MODIS_L2A_OUT_DIR') + \ - " --fname=" + cf('MODIS_L2A_OUT_TEMPLATE') + \ - " --net=" + cf('MODIS_L2A_NN_FILE') + \ - " --aer_x=" + cf('MODIS_L2A_AER_X') + \ - " --blank_ods=" + cf('MODIS_L2A_BLANK_ODS') - - if cf('MODIS_L2A_OVERWRITE').upper() == 'YES': Options += " --force" - if cf('MODIS_L2A_VERBOSE').upper() == 'YES': Options += " -v" - - # Generate products - # ----------------- - i = 0 - Coll = cf('MODIS_L2A_COLLECTION').split(',') - for ident in cf('MODIS_L2A_IDENTS').split(','): - coll = Coll[i] - cmd = "mxd04_l2a.py %s --collection=%s %s %s "%(Options,Coll[i],ident,isotime) - print(cmd) - if not options.dryrun: - if system(cmd): - raise ValueError("mxd04_l2a.py failed for %s on %s "%(ident,isotime)) - - i += 1 - - diff --git a/src/Applications/GAAS_App/mxd04_l2a.py b/src/Applications/GAAS_App/mxd04_l2a.py deleted file mode 100755 index 00e5f600..00000000 --- a/src/Applications/GAAS_App/mxd04_l2a.py +++ /dev/null @@ -1,237 +0,0 @@ -#!/usr/bin/env python3 -# -W ignore::DeprecationWarning - -""" - A Python script to create NNR retrievals. - It now uses class MXD04 to directly read MODIS Aerosol Level 2 - Files (MOD04/MYD04). - - This utility reads MODIS Level2 files and creates an ODS file with - NNR retrievals, as well as a *gritas* type gridded output. - - February 2011, revised Novembre 2016 for MODIS Collection 6. - arlindo.dasilva@nasa.gov -""" - -import warnings -warnings.simplefilter('ignore',DeprecationWarning) -warnings.simplefilter('always',UserWarning) - -import os -import sys -import subprocess - -from optparse import OptionParser # Command-line args -from dateutil.parser import parse as isoparse -from mxd04_nnr import MxD04_NNR -from MAPL.config import strTemplate - -Ident = dict( modo = ('MOD04','ocean'), - modl = ('MOD04','land'), - modd = ('MOD04','deep'), - mydo = ('MYD04','ocean'), - mydl = ('MYD04','land'), - mydd = ('MYD04','deep') - ) - -#--------------------------------------------------------------------- -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 = 'nnr3' - ident = 'modo' - -# Defaults may be platform dependent -# ---------------------------------- - if os.path.exists('/nobackup/MODIS/Level2/'): # New calculon - l2_path = '/nobackup/MODIS/Level2/' - out_dir = '/nobackup/NNR/%coll/Level%lev/%prod/Y%y4/M%m2' - nn_file = '/nobackup/NNR/Net/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 = '006' - res = 'c' - -# Parse command line options -# -------------------------- - parser = OptionParser(usage="Usage: %prog [options] ident isotime", - version='mxd04_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="MODIS 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 MODIS 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.") - - (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(" MODIS Level 2A Processing") - print(" -------------------------") - print("") - -# 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("mxd04_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("mxd04_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 - - -# MODIS Level 2 NNR Aerosol Retrievals -# ------------------------------------ - if options.verbose: - print("NNR Retrieving %s %s on "%(prod,algo.upper()),syn_time) - - modis = MxD04_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 modis.nobs < 1: - if options.verbose: - print('WARNING: no observation for this time in file <%s>'%ods_file) - - elif any(modis.iGood) == False: - if options.verbose: - print('WARNING: no GOOD observation for this time in file <%s>'%ods_file) - modis.nobs = 0 - - nn_file = options.nn_file.replace('%ident',ident) - modis.apply(nn_file) - -# Write ODS -# --------- - makethis_dir(ods_file) - if modis.nobs>0: - modis.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 modis.nobs>0: - if str.isdigit(options.res): - modis.writeg(filename=out_file,refine=int(options.res),channels=modis.channels_) - else: - modis.writeg(filename=out_file,res=options.res,channels=modis.channels_) - -# Write ungridded data -# -------------------- -# name, ext = os.path.splitext(out_file) -# npz_file = name.replace('Level3','Level2') + '.npz' -# makethis_dir(npz_file) -# modis.write(npz_file) - -# Compress nc output unless the user disabled it -# ---------------------------------------------- - if modis.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) diff --git a/src/Applications/GAAS_App/mxd04_nnr.py b/src/Applications/GAAS_App/mxd04_nnr.py deleted file mode 100644 index 5595dc20..00000000 --- a/src/Applications/GAAS_App/mxd04_nnr.py +++ /dev/null @@ -1,541 +0,0 @@ -""" -This module implements the MODIS NNR AOD retrievals. - -This version works from MODIS MOD04/MYD04 Level 2 files. - -""" -import os, sys -import warnings -from pyobs.mxd04 import MxD04_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( META = ( "Scan_Start_Time", - "Latitude", - "Longitude", - "Solar_Zenith", - "Solar_Azimuth", - "Sensor_Zenith", - "Sensor_Azimuth", - "Scattering_Angle", - "Glint_Angle"), - - LAND = ( 'Corrected_Optical_Depth_Land', - 'Mean_Reflectance_Land', - 'Surface_Reflectance_Land', - 'Cloud_Fraction_Land', - 'Quality_Assurance_Land', - 'Deep_Blue_Cloud_Fraction_Land'), - - OCEAN = ( 'Effective_Optical_Depth_Best_Ocean', - 'Mean_Reflectance_Ocean', - 'Cloud_Fraction_Ocean', - 'Quality_Assurance_Ocean'), - - DEEP = ( 'Deep_Blue_Aerosol_Optical_Depth_550_Land', - 'Deep_Blue_Spectral_Aerosol_Optical_Depth_Land', - 'Deep_Blue_Spectral_TOA_Reflectance_Land', - 'Deep_Blue_Spectral_Surface_Reflectance_Land', - 'Deep_Blue_Cloud_Fraction_Land', - 'Deep_Blue_Aerosol_Optical_Depth_550_Land_QA_Flag', - 'Mean_Reflectance_Land', - 'Surface_Reflectance_Land', - 'Aerosol_Cloud_Fraction_Land', - 'Quality_Assurance_Land')) - -ALIAS = dict( Deep_Blue_Aerosol_Optical_Depth_550_Land = 'aod550', - Mean_Reflectance_Land = 'reflectance_lnd', - Surface_Reflectance_Land = 'sfc_reflectance_lnd', - Aerosol_Cloud_Fraction_Land = 'cloud_lnd', - Quality_Assurance_Land = 'qa_lnd' ) - -# Channels for TOA reflectance -# ----------------------------- -CHANNELS = dict ( - LAND = ( 470, 550, 660, 870, 1200, 1600, 2100, 412, 440), - OCEAN = ( 470, 550, 660, 870, 1200, 1600, 2100 ), - DEEP = ( 412, 470, 660 ), - ) - -SCHANNELS = dict ( - LAND = ( 470, 660, 2100 ), - DEEP = ( 412, 470, 660 ), - ) - - -# Translate Inputs between NNR and MODIS classes -# ----------------------------------------------- -TranslateInput = dict ( mRef412 = ('reflectance',412), - mRef440 = ('reflectance',440), - mRef470 = ('reflectance',470), - mRef550 = ('reflectance',550), - mRef660 = ('reflectance',660), - mRef870 = ('reflectance',870), - mRef1200 = ('reflectance',1200), - mRef1600 = ('reflectance',1600), - mRef2100 = ('reflectance',2100), - mSre412 = ('sfc_reflectance',412), - mSre470 = ('sfc_reflectance',470), - mSre660 = ('sfc_reflectance',660), - mSre2100 = ('sfc_reflectance',2100), - ) - -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 MxD04_NNR(MxD04_L2): - """ - This class extends MODIS by adding methods for producing - NNR AOD retrievals based on the Neural Net model developed - with class *abc_c6*. - """ - - def __init__(self,l2_path,prod,algo,syn_time,aer_x, - cloud_thresh=0.70, - glint_thresh=40.0, - scat_thresh=170.0, - cloudFree=None, - aodmax=1.0, - coll='006',verbose=0): - """ - Contructs a MXD04 object from MODIS Aerosol Level 2 - granules. On input, - - l2_path --- top directory for the MODIS Level 2 files; - it must have subdirectories MOD04 and MYD04. - prod --- either *MOD04* (Terra) or *MYD04* (Aqua) - algo --- aerosol algorithm: LAND, OCEAN or DEEP (for - Deep Blue) - 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 - aod_coarse - wind - - It also performs Q/C, setting attribute iGood. On, - input, *cloud_thresh* is the cloud fraction limit. - When DEEP BLUE algorithm is requested, filters for - land 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,prod,syn_time,coll=coll) - if algo != "DEEP": - MxD04_L2.__init__(self,Files,algo,syn_time, - only_good=True, - SDS=SDS, - alias={'Deep_Blue_Cloud_Fraction_Land':'cloud_deep'}, - Verb=verbose) - else: - MxD04_L2.__init__(self,Files,algo,syn_time, - only_good=False, - SDS=SDS, - alias=ALIAS, - Verb=verbose) - - if self.nobs < 1: - return # no obs, nothing to do - - # Reorganize Reflectance Arrays - # ----------------------------- - self.rChannels = CHANNELS[algo] - if algo in SCHANNELS: - self.sChannels = SCHANNELS[algo] - - if algo == "OCEAN": - self.reflectance = self.reflectance[:,0:7] #not using 412, 443, and 745 for now - if algo == "LAND": - self.reflectance = self.reflectance[:,0:-1] #not using 745 for now - - # 3-Ch Algorithm only used when Dark Target data is unavailable - # -------------------------------------------------------------- - - 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.cloud0) - - if algo in SCHANNELS: - for i,c in enumerate(self.sChannels): - self.iGood = self.iGood & (self.sfc_reflectance[:,i]>0) - - if algo == "OCEAN": - self.iGood = self.iGood & (self.GlintAngle > glint_thresh) - - if algo != "OCEAN": - 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 - - # Special handle for brown carbon (treat it as it were organic carbon) - # ---------------------------------------------------- - try: - self.sampleFile(aer_x,onlyVars=('BREXTTAU',),Verbose=Verbose) - self.foc += self.sample.BREXTTAU / s.TOTEXTTAU - self.fcc = self.fbc + self.foc - except: - pass # ignore it for systems without brown carbon as a separate tracer - 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 or DEEP, surface reflectivity - k = list(self.sChannels).index(ch) # index of channel - elif 'mRef' in inputName: # MOD04 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 MODIS 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: - 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/MODIS/Level2/' - algo = 'DEEP' - prod = 'MOD04' - coll = '006' - aer_x = '/nobackup/NNR/Misc/tavg1_2d_aer_Nx' - - syn_time = datetime(2008,6,30,12,0,0) - 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 = MxD04_NNR(l2_path,prod,algo.upper(),syn_time,aer_x, - coll=coll, - cloud_thresh=0.7, - verbose=True) - - m.apply(nn_file) - aod = m.aod_ diff --git a/src/Applications/GAAS_App/patmosx_l2a.py b/src/Applications/GAAS_App/patmosx_l2a.py deleted file mode 100755 index 49e5a027..00000000 --- a/src/Applications/GAAS_App/patmosx_l2a.py +++ /dev/null @@ -1,254 +0,0 @@ -#!/usr/bin/env python3 -# -W ignore::DeprecationWarning - -""" - A Python script to create AVHRR Neural Net Retrievals, - - This utility reads AVHRR Level2 files (or its NPZ conterparts) and - creates an ODS file with the retrieve 550 nm AOD, as well as a - *gritas* type gridded output. - - October 2013 - arlindo.dasilva@nasa.gov -""" - -import warnings -warnings.simplefilter('ignore',DeprecationWarning) -warnings.simplefilter('always',UserWarning) - -import os -import sys - -from time import clock -from optparse import OptionParser # Command-line args -from datetime import datetime -from numpy import load - -from pyobs import NPZ -from avhrr_nnr import AVHRR_NNR -from MAPL import strTemplate - -prod = 'patmosx_v05r02' # baseline product hardwired for now - -#--------------------------------------------------------------------- -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) - -def patmosx_has_obs(fname): - """ - Returns true if PATMOSX reflectance file has something in it. - """ - return 'latitude' in list(load(fname).keys()) - -#--------------------------------------------------------------------- - -if __name__ == "__main__": - - expid = 'nnr_001' - -# Defaults may be platform dependent -# ---------------------------------- - if os.path.exists('/nobackup/AVHRR/Level2/'): # New calculon - l2_path = '/nobackup/AVHRR/Level2/Synoptic' - out_dir = '/nobackup/AVHRR/Level%lev/NNR/Y%y4/M%m2' - blank_ods = '/nobackup/NNR/Misc/blank.ods' - coxmunk_lut = '/nobackup/NNR/Misc/avhrr.cox-munk_lut.npz' - - MAerDir = '/nobackup/MERRAero' - MDir = '/nobackup/MERRA' - nn_file = '/nobackup/NNR/Net/'+expid+'.avhrr_Tau.net' - - else: # Must be somewhere else, no good defaults - out_dir = './' - l2_path = './' - blank_ods = 'blank.ods' - coxmunk_lut = 'cox-munk_lut.npz' - MAerDir = './' - MDir = './' - nn_file = expid+'.avhrr_Tau.net' - - aod_file = MAerDir + '/inst2d_hwl_x.ddf' - tpw_file = MDir + '/int_Nx' - wind_file = MDir + '/slv_Nx' - - out_tmpl = '%s.%prod_l%leva.%orb.%y4%m2%d2_%h2%n2z.%ext' - res = 'e' - -# Parse command line options -# -------------------------- - parser = OptionParser(usage="Usage: %prog [options] asc|des nymd nhms", - version='patmosx_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("-B", "--blank_ods", dest="blank_ods", default=blank_ods, - help="Blank ODS file name for fillers (default=%s)"\ - %blank_ods ) - - 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 AVHRR Level 2 files (default=%s)"\ - %l2_path ) - parser.add_option("-M", "--coxmunk", dest="coxmunk_lut", default=coxmunk_lut, - help="File name for Cox-Munk LUT (default=%s)"\ - %coxmunk_lut ) - - parser.add_option("-N", "--net", dest="nn_file", default=nn_file, - help="Neural net file template (default=%s)"\ - %nn_file ) - - parser.add_option("-A", "--aod", dest="aod_file", default=aod_file, - help="Gridded speciated AOD file (default=%s)"\ - %aod_file ) - - parser.add_option("-T", "--tpw", dest="tpw_file", default=tpw_file, - help="Gridded TPW file (default=%s)"\ - %tpw_file ) - - parser.add_option("-W", "--wind", dest="wind_file", default=wind_file, - help="Gridded Wind file (default=%s)"\ - %wind_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", - help="Do not use n4zip to compress gridded/ODS output file (default=False)") - - parser.add_option("-F", "--force", - action="store_true", dest="force", - help="Overwrites output file") - - parser.add_option("-v", "--verbose", - action="store_true", dest="verbose", - help="Turn on verbosity.") - - (options, args) = parser.parse_args() - - if len(args) == 3: - orb, nymd, nhms = args - else: - parser.error("must have 3 arguments: orbit nymd nhms") - -# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - if options.verbose: - print("") - print(" AVHRR NNR Retrievals") - print(" --------------------") - print("") - t0 = clock() - -# Time variables -# -------------- - nymd_, nhms_ = ( int(nymd), int(nhms) ) - year, month, day = (nymd_/10000, (nymd_%10000)/100, nymd_%100) - hour = nhms_/10000 - syn_tyme = datetime(year,month,day,hour,0,0) # no minute, second - -# Form output gridded file name -# ----------------------------- - out_tmpl = options.out_dir+'/'+options.out_tmpl - out_tmpl = out_tmpl.replace('%prod',prod).replace('%orb',orb).\ - 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) - -# Form ODS file name -# ------------------ - ods_tmpl = options.out_dir+'/'+options.out_tmpl - ods_tmpl = ods_tmpl.replace('%prod',prod).replace('%orb',orb).\ - 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("avhrr_l3a: Output ODS file <%s> exists --- cannot proceed."%ods_file) - raise IOError("Specify --force to overwrite existing output file.") - -# Produce Level 2 AVHRR Aerosol Retrievals -# ---------------------------------------- - if options.verbose: - print("Retrieving AVHRR AOD from %s (%sing) on "%(prod,orb),syn_tyme) - - t = syn_tyme - patmosx = options.l2_path+'/Y%d/M%02d/D%02d/%s.%s.%d%02d%02d_%02dz.npz'%\ - (t.year,t.month,t.day,prod,orb,t.year,t.month,t.day,t.hour) - - options.tpw_file = strTemplate(options.tpw_file,nymd=nymd,nhms=nhms) - options.aod_file = strTemplate(options.aod_file,nymd=nymd,nhms=nhms) - options.wind_file = strTemplate(options.wind_file,nymd=nymd,nhms=nhms) - - # Special handle empty files - # -------------------------- - if patmosx_has_obs(patmosx): - a = AVHRR_NNR(patmosx, Verb=options.verbose) - else: - a = NPZ(patmosx) - - if a.nobs < 1: - if options.verbose: - print('WARNING: no observation for this time in file <%s>'%ods_file) - - # elif any(a.iValid) == False: - elif len(a.tyme[a.iValid]) < 2: - if options.verbose: - print('WARNING: not enough GOOD observation for this time in file <%s>'%ods_file) - a.nobs = 0 - - else: - - # Attach state dependent metadata - # ------------------------------- - a.speciate(options.aod_file) - a.sampleG5(int_x=options.tpw_file,slv_x=options.wind_file) - a.getCoxMunk(options.coxmunk_lut) # oceanic albedo - - # Evaluate Neural Net - # ------------------- - a.apply(options.nn_file) - -# Write ODS -# --------- - if os.path.exists(ods_file) and options.force: - os.remove(ods_file) # remove ODS file - makethis_dir(ods_file) - if a.nobs>0: - a.writeODS(syn_tyme,ods_file,doNNR=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 -# ------------------------- - if os.path.exists(out_file) and options.force: - os.remove(out_file) # remove gridded file - makethis_dir(out_file) - if a.nobs>0: - a.writeGridded(syn_tyme,filename=out_file,res=options.res,doNNR=True,doFilter=False) - -# Compress nc output unless the user disabled it -# ---------------------------------------------- - if a.nobs>0: - if not options.uncompressed: - if os.system("n4zip "+out_file): - warnings.warn('cannot compress output file <%s>'%out_file) - if os.system("n4zip "+ods_file): - warnings.warn('cannot compress output ods file <%s>'%ods_file)