From 97e95434701f2dd1758301321f7cb203f4dbb590 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Thu, 1 Jun 2023 16:55:44 -0400 Subject: [PATCH 01/21] update sample script #! to python3 --- src/NNR/modis/sample_cxalbedo.py | 2 +- src/NNR/modis/sample_mcd43c1.py | 2 +- src/NNR/modis/sample_merra.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/NNR/modis/sample_cxalbedo.py b/src/NNR/modis/sample_cxalbedo.py index e750076..1f7b814 100755 --- a/src/NNR/modis/sample_cxalbedo.py +++ b/src/NNR/modis/sample_cxalbedo.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 """ Precalculate Cox-Munk albedo from sampled wind parameters according to giant file Requires MERRA-2 npz file to be present diff --git a/src/NNR/modis/sample_mcd43c1.py b/src/NNR/modis/sample_mcd43c1.py index ec9042a..c8b1638 100755 --- a/src/NNR/modis/sample_mcd43c1.py +++ b/src/NNR/modis/sample_mcd43c1.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 """ Sample MCD43C1 BRDF parameters according to giant file """ diff --git a/src/NNR/modis/sample_merra.py b/src/NNR/modis/sample_merra.py index 58865ad..c36b575 100755 --- a/src/NNR/modis/sample_merra.py +++ b/src/NNR/modis/sample_merra.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 """ Sample MERRA-2 according to giant file """ From 65f58f2d3351bb51a38ef6e293c542ce52775759 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Thu, 1 Jun 2023 16:56:09 -0400 Subject: [PATCH 02/21] add option to read other vars that uv10m and v10m --- src/NNR/pyabc/giant.py | 37 ++++++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 13 deletions(-) diff --git a/src/NNR/pyabc/giant.py b/src/NNR/pyabc/giant.py index 910ac50..7ce575f 100644 --- a/src/NNR/pyabc/giant.py +++ b/src/NNR/pyabc/giant.py @@ -499,24 +499,35 @@ def sampleG5(self,gas_x=None,avk_x=None,int_x=None,slv_x=None,ext_Nc=None): 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) + FineMode=False,npzFile=None,Verbose=False,onlyVars=('U10M','V10M')): + + if onlyVars is not None: + self.sampleFile(slv_x,onlyVars=onlyVars, Verbose=Verbose) + + labels = () + for varname in onlyVars: + self.__dict__[varname.lower()] = self.sample.__dict__[varname] + labels = labels + (varname.lower(),) + + if 'U10M' in onlyVars: + self.wind = sqrt(self.sample.U10M[:]**2 + self.sample.V10M[:]**2) del self.sample - self.speciate(aer_x,FineMode=FineMode,Verbose=Verbose) + if aer_x is not None: + self.speciate(aer_x,FineMode=FineMode,Verbose=Verbose) + + if aer_x is not None: + labels = labels + ('fdu','fss','fcc','fsu') + if FinMode: + labels = labels + ('fduf','fssf') + + kwds = {} + for varname in labels: + kwds[varname] = self.__dict__[varname] 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) + savez(npzFile,**kwds) def sampleMCD43C(self,npzFile=None,Verbose=False): from pyabc.mcd43c import MCD43C From 3548ba4eb820eeaa7f12956c7771c9ea3b9d6ad9 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Thu, 1 Jun 2023 16:56:48 -0400 Subject: [PATCH 03/21] script to sample tqv and to3 only from merra-2 --- src/NNR/modis/sample_merra_tqv_to3.py | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100755 src/NNR/modis/sample_merra_tqv_to3.py diff --git a/src/NNR/modis/sample_merra_tqv_to3.py b/src/NNR/modis/sample_merra_tqv_to3.py new file mode 100755 index 0000000..d5272f5 --- /dev/null +++ b/src/NNR/modis/sample_merra_tqv_to3.py @@ -0,0 +1,25 @@ +#!/usr/bin/env python3 +""" + Sample MERRA-2 according to giant file +""" + +import os, sys +from pyabc.giant import GIANT +import argparse + + +if __name__ == '__main__': + # Parse command line options + # -------------------------- + parser = argparse.ArgumentParser() + parser.add_argument("giantFile",help="giant spreadhseet file to base sampling on") + + args = parser.parse_args() + outDir = os.path.dirname(args.giantFile) + fname = os.path.basename(args.giantFile) + npzFile = outDir + '/' + fname[:-3] + '_MERRA2_TQV_TO3.npz' + + g = GIANT(args.giantFile) + g.sampleMERRA(npzFile=npzFile,Verbose=False,aer_x=None,onlyVars=('TQV','TO3')) + + From bb063781e5d360cfc4eee5c33616e2b13c34264e Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Thu, 1 Jun 2023 17:08:59 -0400 Subject: [PATCH 04/21] add feature to read tqv_to3 file --- src/NNR/pyabc/abc_c6.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/NNR/pyabc/abc_c6.py b/src/NNR/pyabc/abc_c6.py index 8831d55..91afafa 100644 --- a/src/NNR/pyabc/abc_c6.py +++ b/src/NNR/pyabc/abc_c6.py @@ -188,6 +188,7 @@ def __init__(self,fname,Albedo,coxmunk_lut=None,NDVI=False): # Get Auxiliary Data self.fnameRoot = fname[:-3] self.setWind() + self.setTQV_TO3() self.setAlbedo(Albedo,coxmunk_lut=coxmunk_lut) self.setSpec() if NDVI: @@ -200,6 +201,13 @@ def setWind(self): self.giantList.append('wind') self.Wind = '' #need this for backwards compatibility + def setTQV_TO3(self): + # Read in MERRA-2 sampled TQV and TO3 + # ------------------------------------ + for name in ('tqv','to3'): + self.__dict__[name] = load(self.fnameRoot + "_MERRA2_TQV_TO3.npz")[name] + self.giantList.append(name) + def setAlbedo(self,Albedo,coxmunk_lut=None): # Define wind speed dependent ocean albedo # ---------------------------------------- From 42c3daf76e1f59d12061f44b15317757166e04e6 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Thu, 1 Jun 2023 17:09:28 -0400 Subject: [PATCH 05/21] script to run sample tqv_to3 --- src/NNR/modis/run_sample_merra_tqv_to3.sh | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) create mode 100755 src/NNR/modis/run_sample_merra_tqv_to3.sh diff --git a/src/NNR/modis/run_sample_merra_tqv_to3.sh b/src/NNR/modis/run_sample_merra_tqv_to3.sh new file mode 100755 index 0000000..c5aa222 --- /dev/null +++ b/src/NNR/modis/run_sample_merra_tqv_to3.sh @@ -0,0 +1,20 @@ +#!/usr/bin/sh +#======================================================================= +# name - run_sample_merra.sh +# purpose - +# This script run the sample_merra.py script +# this reads the giant spreadsheet files of co-located MODIS-AERONET +# and samples MERRA-2 at these locations +# + + +# Terra +#--------- +terraFile=/nobackup/NNR/Training/061_py2/giant_C061_10km_Terra_v3.0_20211231.nc +nohup python3 -u ./sample_merra_tqv_to3.py ${terraFile} >& nohup.terra.merra.log & + +# Aqua +#---------- +aquaFile=/nobackup/NNR/Training/061_py2/giant_C061_10km_Aqua_v3.0_20211231.nc +nohup python3 -u ./sample_merra_tqv_to3.py ${aquaFile} >& nohup.aqua.merra.log & + From 8c5b87786102c574d006870023ce85c73babe1e0 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Thu, 1 Jun 2023 17:09:58 -0400 Subject: [PATCH 06/21] modify modis cmakelists to install sample tqv_to3 scripts --- src/NNR/modis/CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/NNR/modis/CMakeLists.txt b/src/NNR/modis/CMakeLists.txt index 0fa2834..908639d 100644 --- a/src/NNR/modis/CMakeLists.txt +++ b/src/NNR/modis/CMakeLists.txt @@ -9,9 +9,11 @@ set (PYSCRIPTS run_sample_cxalbedo.sh run_sample_mcd43c1.sh run_sample_merra.sh + run_sample_merra_tqv_to3.sh sample_cxalbedo.py sample_mcd43c1.py sample_merra.py + sample_merra_tqv_to3.py setup_env.sh ) From 75a23b3f364b93f6fdc913b4149b79af5a93ecfb Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Wed, 7 Jun 2023 14:32:38 -0400 Subject: [PATCH 07/21] updates to train and test on the angstrom linear fit --- src/NNR/pyabc/abc_c6.py | 44 +++++++++++-- src/NNR/pyabc/abc_c6_aux.py | 126 ++++++++++++++++++++++++++++++++++++ 2 files changed, 164 insertions(+), 6 deletions(-) diff --git a/src/NNR/pyabc/abc_c6.py b/src/NNR/pyabc/abc_c6.py index 91afafa..e27ffa5 100644 --- a/src/NNR/pyabc/abc_c6.py +++ b/src/NNR/pyabc/abc_c6.py @@ -22,7 +22,7 @@ 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 .abc_c6_aux import make_plots, make_plots_angstrom, make_plots_angstrom_fit, TestStats, SummaryPDFs from .brdf import rtlsReflectance from .mcd43c import BRDF from functools import reduce @@ -91,15 +91,19 @@ def setupNN(self,retrieval,expid, # figure out if you need to calculate angstrom exponent angstrom = False + angstrom_fit = False for tname in Target: if not angstrom: - if 'AE' in tname: + if ('AE' in tname) and ('AEfit' not in tname): angstrom = True + if not angstrom_fit: + if 'AEfit' in tname: + angstrom_fit = True self.angstrom = angstrom + self.angstrom_fit = angstrom_fit # if angstrom is being trained # find the base wavelength - # calculate angstrom with respect to the base wavelength # ------------------------------------------------------- if angstrom: # find base wavelength @@ -114,6 +118,10 @@ def setupNN(self,retrieval,expid, self.AE_base_wav = base_wav self.AE_base_wav_i = base_wav_i + # if angstrom is being trained + # calculate angstrom with respect to the base wavelength + # ------------------------------------------------------- + if angstrom: # Calculate the angstrom exponent # with respect to the base wavelength for tname in Target: @@ -124,6 +132,26 @@ def setupNN(self,retrieval,expid, AE = -1.*np.log(tau/base_tau)/np.log(wav/base_wav) self.__dict__['aAE'+wavs] = AE + # if angstrom_fit is being trained + # calculate a linear fit to the log-log + # of wavelength and AOD 440-870 + # ------------------------------------------------------- + if angstrom_fit: + wavs = ['440','470','500','550','660','870'] + wav = np.array(wavs).astype(float) + + # Calculate the angstrom exponent + # with a linear fit to log AOD + tau = [] + for w in wavs: + tau.append(self.__dict__['aTau'+w]) + + tau = np.array(tau) + fit = np.polyfit(np.log(wav),-1.*np.log(tau+0.01),1) + self.__dict__['aAEfitm'] = fit[0,:] + self.__dict__['aAEfitb'] = fit[1,:] + + # Balance the dataset before splitting # No aerosol type should make up more that 35% # of the total number of obs @@ -205,7 +233,7 @@ def setTQV_TO3(self): # Read in MERRA-2 sampled TQV and TO3 # ------------------------------------ for name in ('tqv','to3'): - self.__dict__[name] = load(self.fnameRoot + "_MERRA2_TQV_TO3.npz")[name] + self.__dict__[name] = load(self.fnameRoot + "_MERRA2_TQV_TO3.npz")[name]*0.01 self.giantList.append(name) def setAlbedo(self,Albedo,coxmunk_lut=None): @@ -1180,6 +1208,8 @@ def _test(mxd,expid,c,plotting=True): if plotting: if mxd.angstrom: make_plots_angstrom(mxd,expid,ident,I=mxd.iTest) + elif mxd.angstrom_fit: + make_plots_angstrom_fit(mxd,expid,ident,I=mxd.iTest) else: make_plots(mxd,expid,ident,I=mxd.iTest) else: @@ -1217,6 +1247,8 @@ def _test(mxd,expid,c,plotting=True): if plotting: if mxd.angstrom: make_plots_angstrom(mxd,expid,'.k={}'.format(str(k)),I=mxd.iTest) + elif mxd.angstrom_ft: + make_plots_angstrom_fit(mxd,expid,'.k={}'.format(str(k)),I=mxd.iTest) else: make_plots(mxd,expid,'.k={}'.format(str(k)),I=mxd.iTest) k = k + 1 @@ -1295,10 +1327,10 @@ def flatten_list(Input_nnr): Albedo = ['CoxMunkBRF'] K = 2 - if sat is 'Terra': + if sat == 'Terra': filename = '/nobackup/6/NNR/Training/giant_C6_10km_Terra_20150921.nc' retrieval = 'MOD_OCEAN' - if sat is 'Aqua': + if sat == 'Aqua': filename = '/nobackup/6/NNR/Training/giant_C6_10km_Aqua_20151005.nc' retrieval = 'MYD_OCEAN' diff --git a/src/NNR/pyabc/abc_c6_aux.py b/src/NNR/pyabc/abc_c6_aux.py index 887e045..fcc5707 100644 --- a/src/NNR/pyabc/abc_c6_aux.py +++ b/src/NNR/pyabc/abc_c6_aux.py @@ -573,6 +573,132 @@ def make_plots_angstrom(mxd,expid,ident,I=None): savefig(outdir+"/"+expid+"."+ident+"_kde-AE"+name[4:]+'.png') plt.close(fig) + +#--------------------------------------------------------------------- + +def make_plots_angstrom_fit(mxd,expid,ident,I=None): + outdir = mxd.plotdir + if not os.path.exists(outdir): + os.makedirs(outdir) + if I is None: + I = ones(mxd.lon.shape).astype(bool) + + # get the targets and predictions + # ------------------------------------------------ + targets = mxd.getTargets(I) + results = mxd.eval(I) + + # see if there are any AOD predictions + # plot KDE of corrected AOD + # -------------------------- + for t in range(mxd.nTarget): + tname = mxd.Target[t] + if 'Tau' in tname: + fig = _plotKDE(targets[:,t],results[:,t],y_label='NNR') + title("Log(tau"+tname[4:]+"+0.01)- "+ident) + savefig(outdir+"/"+expid+"."+ident+"_kde-Tau"+tname[4:]+'-corrected.png') + plt.close(fig) + + + # Plot KDE of predicted AE fit coefficients + for t in range(mxd.nTarget): + tname = mxd.Target[t] + if 'AEfit' in tname: + if 'AEfitm' in tname: + x_bins=np.arange(-1,3,0.1) + else: + x_bins=np.arange(-20,10,0.5) + fig = _plotKDE(targets[:,t],results[:,t],y_label='NNR',x_bins=x_bins) + title("Log-Log AOD Linear Fit Coefficient "+tname[-1]+" "+ident) + savefig(outdir+"/"+expid+"."+ident+"_kde-"+tname[1:]+'-corrected.png') + plt.close(fig) + if 'AEfitm' in tname: + AEmt = targets[:,t] + AEmr = results[:,t] + elif 'AEfitb' in tname: + AEbt = targets[:,t] + AEbr = results[:,t] + + # get the AOD from the predicted angstrom + # linear fit + # ------------------------------------------------ + nwav = 6 + wavs = ['440','470','500','550','660','870'] + wav = np.array(wavs).astype(float) + + nobs, nt = targets.shape + targets_ = np.zeros([nobs,nwav]) + results_ = np.zeros([nobs,nwav]) + for i in range(nwav): + tau = mxd.__dict__['aTau'+wavs[i]] + targets_[:,i] = np.log(tau + 0.01) + + results_[:,i] = -1.*(AEmr*np.log(wav[i]) + AEbr) + + # Plot KDE of corrected AOD from AEfit + # ------------------------------------ + # loop through targets + for t in range(nwav): + fig = _plotKDE(targets_[:,t],results_[:,t],y_label='NNR') + title("Log(tau"+wavs[t]+"+0.01)- "+ident) + savefig(outdir+"/"+expid+"."+ident+"_kde-Tau"+wavs[t]+'-fitcorrected.png') + plt.close(fig) + + + # Plot KDE of uncorrected AOD + # --------------------------- + # loop through targets + # plot if there's a corresponding MODIS retrieval + wavm = [] + for t in range(nwav): + name = 'mTau'+wavs[t] + if name in mxd.__dict__: + wavm.append(t) + original = log(mxd.__dict__[name][I]+0.01) + + # 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(tau"+wavs[t]+"+0.01)- "+ident) + savefig(outdir+"/"+expid+"."+ident+"_kde-"+name[1:]+'.png') + plt.close(fig) + + if name == 'mTau550': + if 'aTau550' in mxd.Target: + for t in range(mxd.nTarget): + tname = mxd.Target[t] + if tname == 'aTau550': + # Scatter diagram for testing + # --------------------------- + mxd.plotScat(iTarget=t,I=I,figfile=outdir+"/"+expid+"."+ident+"_scat-"+name+'.png') + + # Get AE of MODIS + #----------------- + mdata = [] + fwav = [] + for t in wavm: + name = 'mTau'+wavs[t] + fwav.append(wav[t]) + mdata.append(mxd.__dict__[name][I]) + + mdata = np.array(mdata) + fwav = np.array(fwav) + fit = np.polyfit(np.log(fwav),-1.*np.log(mdata+0.01),1) + AEm = fit[0,:] + AEb = fit[1,:] + + fig = _plotKDE(AEmt,AEm,y_label='Original Modis',x_bins=np.arange(-1,3,0.1)) + title("AE 440-870 " +ident) + savefig(outdir+"/"+expid+"."+ident+"_kde-AEfitm.png") + plt.close(fig) + + fig = _plotKDE(AEbt,AEb,y_label='Original Modis',x_bins=np.arange(-20,10,0.5)) + title("AE intercept 440-870 " +ident) + savefig(outdir+"/"+expid+"."+ident+"_kde-AEfitb.png") + plt.close(fig) + + #--------------------------------------------------------------------- def make_error_pdfs(mxd,Input,expid,ident,K=None,I=None,Title=None,netfileRoot=None, emin=-1.5,emax=2.5): From 77d061de509e3c1a21dc5faf92574030f59b86c5 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Wed, 7 Jun 2023 14:34:18 -0400 Subject: [PATCH 08/21] use develop GMAOpyobs branch. this is temporary --- components.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components.yaml b/components.yaml index 9899863..3bf0ada 100644 --- a/components.yaml +++ b/components.yaml @@ -22,5 +22,5 @@ ecbuild: GMAOpyobs: local: ./src/Shared/GMAO_Shared/GMAO_pyobs@ remote: ../GMAOpyobs.git - branch: fix/pcastell/no-types-module + branch: develop develop: develop From b97a92da46f361baa9833a0db8fcfe93be7eed33 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Wed, 7 Jun 2023 15:32:43 -0400 Subject: [PATCH 09/21] update GMAOpyobs to v1.0.4 --- CHANGELOG.md | 4 ++++ components.yaml | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4506e88..c38b44c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,8 +4,12 @@ ### Added +- ability to train on the AE linear fit to spectral AOD + ### Fixed +- updated GMAOpyobs to latest v1.0.4 + ### Changed ## [v1.0.0] 2023-06-07 diff --git a/components.yaml b/components.yaml index 3bf0ada..843e45f 100644 --- a/components.yaml +++ b/components.yaml @@ -22,5 +22,5 @@ ecbuild: GMAOpyobs: local: ./src/Shared/GMAO_Shared/GMAO_pyobs@ remote: ../GMAOpyobs.git - branch: develop + branch: v1.0.4 develop: develop From c9ea144c27cb576d5567f5e0aca5689492eb387f Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Tue, 22 Aug 2023 11:38:58 -0400 Subject: [PATCH 10/21] update GMAOpyobs to v1.0.5 --- CHANGELOG.md | 2 +- components.yaml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ea25083..98f0caf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,7 +8,7 @@ ### Fixed -- updated GMAOpyobs to latest v1.0.4 +- updated GMAOpyobs to latest v1.0.5 ### Changed diff --git a/components.yaml b/components.yaml index 843e45f..5cdc5c7 100644 --- a/components.yaml +++ b/components.yaml @@ -22,5 +22,5 @@ ecbuild: GMAOpyobs: local: ./src/Shared/GMAO_Shared/GMAO_pyobs@ remote: ../GMAOpyobs.git - branch: v1.0.4 + tag: v1.0.5 develop: develop From a91bbc23436dd99c3117b686a5764af53c3d9abc Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Wed, 6 Sep 2023 10:50:40 -0400 Subject: [PATCH 11/21] update changelog for added tqv and to3 --- CHANGELOG.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 98f0caf..03ffbeb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,13 +5,14 @@ ### Added - ability to train on the AE linear fit to spectral AOD +- use TQV and TO3 as inputs ### Fixed -- updated GMAOpyobs to latest v1.0.5 - ### Changed +- updated GMAOpyobs to latest v1.0.5 + ## [v1.0.0] 2023-06-07 ### Added From 5f542256a8481b848f32664ff4d125a5dc3539d9 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Fri, 15 Dec 2023 14:53:22 -0500 Subject: [PATCH 12/21] add sknet to this repo --- src/CMakeLists.txt | 8 ++++++++ src/pyabc/__init__.py | 5 +++++ src/pyabc/sknet.py | 48 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 61 insertions(+) create mode 100644 src/pyabc/__init__.py create mode 100644 src/pyabc/sknet.py diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 96d7246..51969a5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -6,3 +6,11 @@ # ------------------------------------- esma_add_subdirectory(Shared) esma_add_subdirectory(NNR) + +# Install 'abc' package +# --------------------- +file (GLOB abc_files pyabc/*.py) +install ( + FILES ${abc_files} + DESTINATION lib/Python/pyabc + ) diff --git a/src/pyabc/__init__.py b/src/pyabc/__init__.py new file mode 100644 index 0000000..fe013a1 --- /dev/null +++ b/src/pyabc/__init__.py @@ -0,0 +1,5 @@ +""" +ABC utilities +""" + + diff --git a/src/pyabc/sknet.py b/src/pyabc/sknet.py new file mode 100644 index 0000000..8b498d0 --- /dev/null +++ b/src/pyabc/sknet.py @@ -0,0 +1,48 @@ +""" + + Simple wrapper around FFNET to make it look like a sklearn estimator. + It also provides sklearn style cross_validation + +""" + +from ffnet import * +from sklearn.base import BaseEstimator +from sklearn.model_selection import cross_val_score + +class SKNET(ffnet,BaseEstimator): + + def predict(self,X): + """ + Evaluate the Feed-Forward Neural Network. + """ + y_ = self(X) + return y_.ravel() + + def fit(self,X,y,**kwopts): + """ + Train the Feed-Forward Neural Net. + """ + self.train_tnc(X,y,**kwopts) + + def score(self,X,y): + """ + Returns the coefficient of determination R^2 of the prediction. + + The coefficient R^2 is defined as (1 - u/v), where u is the + regression sum of squares ((y - y_pred) ** 2).sum() and v is the + residual sum of squares ((y_true - y_true.mean()) ** 2).sum(). + Best possible score is 1.0, lower values are worse. + """ + y_ = self.predict(X) + u = ((y_-y)**2).sum() + v = ((y-y.mean())**2).sum() + return (1 - u/v) + + def cross_validate(self,X,y,**kwopts): + """ + Return cross-validation scores. See cross_val_score() for + optional parameters. + scores = net.cross_validate(X,y,cv=10) + """ + scores = cross_val_score(self, X, y, **kwopts) + return scores From 62b9289a995925ce3cd5e754e684497ffee60076 Mon Sep 17 00:00:00 2001 From: Ricardo Todling Date: Sat, 13 Jan 2024 09:00:15 -0500 Subject: [PATCH 13/21] per Patricia guidance --- src/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 51969a5..9b507aa 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -12,5 +12,5 @@ file (GLOB abc_files pyabc/*.py) install ( FILES ${abc_files} - DESTINATION lib/Python/pyabc - ) + DESTINATION lib/Python/pyobs + ) From 0a8be43ea0aabc3179ebb4b426e4cd4dcb4faffa Mon Sep 17 00:00:00 2001 From: Ricardo Todling Date: Sat, 13 Jan 2024 09:02:50 -0500 Subject: [PATCH 14/21] comment on change --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 03ffbeb..c1efb98 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,8 @@ ### Fixed +- Minor fix associated with installation location of software. + ### Changed - updated GMAOpyobs to latest v1.0.5 From ca4968b28365cf7e026eb0533d9693a129784a25 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Tue, 24 Sep 2024 11:42:22 -0400 Subject: [PATCH 15/21] Update to ESMA_env v4.8.2 --- CHANGELOG.md | 1 + components.yaml | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 03ffbeb..20cbc7e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ ### Changed - updated GMAOpyobs to latest v1.0.5 +- Update to ESMA_env v4.8.2 (Fixes for RHEL8 GMAO Machines) ## [v1.0.0] 2023-06-07 diff --git a/components.yaml b/components.yaml index 5cdc5c7..21aba3e 100644 --- a/components.yaml +++ b/components.yaml @@ -5,7 +5,7 @@ AeroML: env: local: ./env@ remote: ../ESMA_env.git - tag: v4.8.0 + tag: v4.8.2 develop: main cmake: From 67e7c2cb497a1e1ec44f45eca5087bf131b1a5fb Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Wed, 6 Nov 2024 10:31:53 -0500 Subject: [PATCH 16/21] Update GitHub Actions Soon, GitHub will deprecate the upload-artifacts v3 action: https://github.blog/changelog/2024-04-16-deprecation-notice-v3-of-the-artifact-actions/ This PR updates the action to the latest version, v4. It also updates the checkout action to the latest version, v4, as well. --- .github/workflows/validate_yaml_files.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/validate_yaml_files.yml b/.github/workflows/validate_yaml_files.yml index 86df2bb..449db6e 100644 --- a/.github/workflows/validate_yaml_files.yml +++ b/.github/workflows/validate_yaml_files.yml @@ -15,7 +15,7 @@ jobs: validate-YAML: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3.3.0 + - uses: actions/checkout@v4 - id: yaml-lint name: yaml-lint uses: ibiqlik/action-yamllint@v3 @@ -24,7 +24,7 @@ jobs: format: colored config_file: .yamllint.yml - - uses: actions/upload-artifact@v3 + - uses: actions/upload-artifact@v4 if: always() with: name: yamllint-logfile From c072132a75e8bba3dabc77e5db3be837d50e4f13 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Wed, 2 Jul 2025 13:53:05 -0400 Subject: [PATCH 17/21] new cmakelists for GAAS_App --- src/Applications/GAAS_App/CMakeLists.txt | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) create mode 100644 src/Applications/GAAS_App/CMakeLists.txt diff --git a/src/Applications/GAAS_App/CMakeLists.txt b/src/Applications/GAAS_App/CMakeLists.txt new file mode 100644 index 0000000..97609ae --- /dev/null +++ b/src/Applications/GAAS_App/CMakeLists.txt @@ -0,0 +1,17 @@ +esma_set_this() + +set(pythonscripts modis_l2a.py mxd04_l2a.py avhrr_l2a.py patmosx_l2a.py) +install(PROGRAMS ${pythonscripts} DESTINATION bin) + +install( + FILES mxd04_nnr.py avhrr_nnr.py + DESTINATION lib/Python + ) + + +file(GLOB pcf_files *.pcf) + +install ( + FILES ${pcf_files} + DESTINATION etc + ) From c6a96c10954748f9468b2f8e21d76215b0ef00bd Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Wed, 2 Jul 2025 13:53:42 -0400 Subject: [PATCH 18/21] add gaas_app scripts from latest version of GEOSadas --- src/Applications/GAAS_App/aeronet_all.py | 86 ++++ src/Applications/GAAS_App/avhrr_all.py | 46 ++ src/Applications/GAAS_App/avhrr_l2a.pcf | 25 ++ src/Applications/GAAS_App/avhrr_l2a.py | 64 +++ src/Applications/GAAS_App/avhrr_nnr.py | 153 +++++++ src/Applications/GAAS_App/misr_land.py | 139 ++++++ src/Applications/GAAS_App/modis_l2a.pcf | 25 ++ src/Applications/GAAS_App/modis_l2a.py | 63 +++ src/Applications/GAAS_App/mxd04_l2a.py | 237 ++++++++++ src/Applications/GAAS_App/mxd04_nnr.py | 541 +++++++++++++++++++++++ src/Applications/GAAS_App/patmosx_l2a.py | 254 +++++++++++ 11 files changed, 1633 insertions(+) create mode 100755 src/Applications/GAAS_App/aeronet_all.py create mode 100755 src/Applications/GAAS_App/avhrr_all.py create mode 100644 src/Applications/GAAS_App/avhrr_l2a.pcf create mode 100755 src/Applications/GAAS_App/avhrr_l2a.py create mode 100644 src/Applications/GAAS_App/avhrr_nnr.py create mode 100755 src/Applications/GAAS_App/misr_land.py create mode 100644 src/Applications/GAAS_App/modis_l2a.pcf create mode 100755 src/Applications/GAAS_App/modis_l2a.py create mode 100755 src/Applications/GAAS_App/mxd04_l2a.py create mode 100644 src/Applications/GAAS_App/mxd04_nnr.py create mode 100755 src/Applications/GAAS_App/patmosx_l2a.py diff --git a/src/Applications/GAAS_App/aeronet_all.py b/src/Applications/GAAS_App/aeronet_all.py new file mode 100755 index 0000000..23a39ea --- /dev/null +++ b/src/Applications/GAAS_App/aeronet_all.py @@ -0,0 +1,86 @@ +#!/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 new file mode 100755 index 0000000..4085ba4 --- /dev/null +++ b/src/Applications/GAAS_App/avhrr_all.py @@ -0,0 +1,46 @@ +#!/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 new file mode 100644 index 0000000..867064e --- /dev/null +++ b/src/Applications/GAAS_App/avhrr_l2a.pcf @@ -0,0 +1,25 @@ +#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 new file mode 100755 index 0000000..9ffbf0c --- /dev/null +++ b/src/Applications/GAAS_App/avhrr_l2a.py @@ -0,0 +1,64 @@ +#!/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 new file mode 100644 index 0000000..0266cad --- /dev/null +++ b/src/Applications/GAAS_App/avhrr_nnr.py @@ -0,0 +1,153 @@ +""" +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 new file mode 100755 index 0000000..f41311d --- /dev/null +++ b/src/Applications/GAAS_App/misr_land.py @@ -0,0 +1,139 @@ +#!/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 new file mode 100644 index 0000000..8c11dee --- /dev/null +++ b/src/Applications/GAAS_App/modis_l2a.pcf @@ -0,0 +1,25 @@ +# 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 new file mode 100755 index 0000000..440bd65 --- /dev/null +++ b/src/Applications/GAAS_App/modis_l2a.py @@ -0,0 +1,63 @@ +#!/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 new file mode 100755 index 0000000..00e5f60 --- /dev/null +++ b/src/Applications/GAAS_App/mxd04_l2a.py @@ -0,0 +1,237 @@ +#!/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 new file mode 100644 index 0000000..5595dc2 --- /dev/null +++ b/src/Applications/GAAS_App/mxd04_nnr.py @@ -0,0 +1,541 @@ +""" +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 new file mode 100755 index 0000000..49e5a02 --- /dev/null +++ b/src/Applications/GAAS_App/patmosx_l2a.py @@ -0,0 +1,254 @@ +#!/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) From c5fcc39d980b0cabfe3bb5455bc229cfe12253cb Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Wed, 2 Jul 2025 13:59:15 -0400 Subject: [PATCH 19/21] update cmake to install new applications directory --- src/Applications/CMakeLists.txt | 1 + src/CMakeLists.txt | 1 + 2 files changed, 2 insertions(+) create mode 100644 src/Applications/CMakeLists.txt diff --git a/src/Applications/CMakeLists.txt b/src/Applications/CMakeLists.txt new file mode 100644 index 0000000..dee8ecf --- /dev/null +++ b/src/Applications/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory (GAAS_App) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 9b507aa..5f467bc 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -6,6 +6,7 @@ # ------------------------------------- esma_add_subdirectory(Shared) esma_add_subdirectory(NNR) + esma_add_subdirectory(Applications) # Install 'abc' package # --------------------- From 6a62318a81ba098ca10b0d1587efb7ce489814be Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Wed, 2 Jul 2025 15:47:09 -0400 Subject: [PATCH 20/21] update Changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 843fcf9..148337e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,7 @@ ### Added +- create an Applications folder for GAAS_App scripts that do the NNR pre-processing for the DAS - ability to train on the AE linear fit to spectral AOD - use TQV and TO3 as inputs From 6b9357754480f86fd03261a4ad835b14bbdc3195 Mon Sep 17 00:00:00 2001 From: Patricia Castellanos Date: Wed, 2 Jul 2025 15:52:48 -0400 Subject: [PATCH 21/21] update changelog for new release --- CHANGELOG.md | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 148337e..8aabc57 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,14 @@ ### Added +### Fixed + +### Changed + +## [v1.1.0] 2025-07-02 + +### Added + - create an Applications folder for GAAS_App scripts that do the NNR pre-processing for the DAS - ability to train on the AE linear fit to spectral AOD - use TQV and TO3 as inputs