From a266319e4ca2fc4e2ec8a77bc5f564468f2e4fab Mon Sep 17 00:00:00 2001 From: roarkmd Date: Thu, 22 Jan 2026 14:29:53 -0500 Subject: [PATCH 01/10] Added geoschem 14.6.2 speciation mechanism. --- cmaq2hemco/mechrc/__init__.py | 4 +- cmaq2hemco/mechrc/geoschem_14_6_2.py | 82 ++++++++++++++++++++++++++++ 2 files changed, 84 insertions(+), 2 deletions(-) create mode 100644 cmaq2hemco/mechrc/geoschem_14_6_2.py diff --git a/cmaq2hemco/mechrc/__init__.py b/cmaq2hemco/mechrc/__init__.py index f786df4..0ab3024 100644 --- a/cmaq2hemco/mechrc/__init__.py +++ b/cmaq2hemco/mechrc/__init__.py @@ -1,3 +1,3 @@ -__all__ = ['cb6r5_ae7_aq'] +__all__ = ['cb6r5_ae7_aq', 'geoschem_14_6_2'] -from . import cb6r5_ae7_aq +from . import cb6r5_ae7_aq, geoschem_14_6_2 diff --git a/cmaq2hemco/mechrc/geoschem_14_6_2.py b/cmaq2hemco/mechrc/geoschem_14_6_2.py new file mode 100644 index 0000000..c4a5c93 --- /dev/null +++ b/cmaq2hemco/mechrc/geoschem_14_6_2.py @@ -0,0 +1,82 @@ +__all__ = ['writeconfig'] + + +cq2gc = { + 'ACET': [['ACET', '']], + 'ACR': [['ACR', '']], + 'ACTA': [['ACTA', '']], + 'ALD2': [['ALD2', '']], + 'ALK4': [['ALK4', '']], + 'ALK6': [['ALK6', '']], + 'BALD': [['BALD', '']], + 'BC': [['BCPI', '70'], ['BCPO', '71']], + 'BENZ': [['BENZ', '']], + 'C2H2': [['C2H2', '']], + 'C2H4': [['C2H4', '']], + 'C2H6': [['C2H6', '']], + 'C3H8': [['C3H8', '']], + 'C4H6': [['C4H6', '']], + 'CCl4': [['CCl4', '']], + 'CH2Br2': [['CH2Br2', '']], + 'CH2Cl2': [['CH2Cl2', '']], + 'CH2O': [['CH2O', '']], + 'CH3Br': [['CH3Br', '']], + 'CH3Cl': [['CH3Cl', '']], + 'CH3I': [['CH3I', '']], + 'CH4': [['CH4', '']], + 'CHCl3': [['CHCl3', '']], + 'CL2': [['CL2', '']], + 'CO': [['CO', '']], + 'CSL': [['CSL', '']], + 'EBZ': [['EBZ', '']], + 'EOH': [['EOH', '']], + 'FURA': [['FURA', '']], + 'GLYX': [['GLYX', '']], + 'HACTA': [['HACTA', '']], + 'HCL': [['HCL', '']], + 'HCOOH': [['HCOOH', '']], + 'HONO': [['HONO', '']], + 'ISOP': [['ISOP', '']], + 'IVOC': [['IVOC', '']], + 'LIMO': [['LIMO', '']], + 'MACR': [['MACR', '']], + 'MEK': [['MEK', '']], + 'MGLY': [['MGLY', '']], + 'MOH': [['MOH', '']], + 'MTPA': [['MTPA', '']], + 'MTPO': [['MTPO', '']], + 'MVK': [['MVK', '']], + 'NAP': [['NAP', '']], + 'NH3': [['NH3', '']], + 'NH4': [['NH4', '']], + 'NIT': [['NIT', '']], + 'NO': [['NO', '115']], + 'NO2': [['NO2', '']], + 'OC': [['OCPI', '72'], ['OCPO', '73']], + 'OCS': [['OCS', '']], + 'pFe': [['pFe', '']], + 'PHEN': [['PHEN', '']], + 'PNA': [['PNA', '']], + 'PRPE': [['PRPE', '']], + 'RCHO': [['RCHO', '']], + 'RCOOH': [['RCOOH', '']], + 'ROH': [['ROH', '']], + 'SO2': [['SO2', '']], + 'SO4': [['SO4', '']], + 'STYR': [['STYR', '']], + 'SULF': [['SULF', '']], + 'TMB': [['TMB', '']], + 'TOLU': [['TOLU', '']], + 'XYLE': [['XYLE', '']], +} +# ignore special species, inventory meta variables, HAP tracers, UNK/UNR +for key in ['NH3_FERT','HFLUX', 'VOC_INV','NMOG','CH4_INV','ACROLEIN','BUTADIENE13','ETHYLBENZ','UNK','UNR']: + cq2gc[key] = [] + +# Ignore unused particulate species +for key in ['PAL', 'PCA', 'PCL', 'PH2O', 'PK', 'PSI', 'PTI','PMC', 'PMG', 'PMN', 'PMOTHR', 'PNCOM']: + cq2gc[key] = [] + +def writeconfig(outpath, year, sector, filepatt): + from .core import writeconfig as genwriteconfig + return genwriteconfig(outpath, year, sector, filepatt, cq2gc) From b3228337144bf806bbe24fff71ec6063a4c11cbe Mon Sep 17 00:00:00 2001 From: roarkmd Date: Fri, 23 Jan 2026 14:17:22 -0500 Subject: [PATCH 02/10] Improved error messaging. --- cmaq2hemco/mechrc/core.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmaq2hemco/mechrc/core.py b/cmaq2hemco/mechrc/core.py index 37cd90a..ae09562 100644 --- a/cmaq2hemco/mechrc/core.py +++ b/cmaq2hemco/mechrc/core.py @@ -17,7 +17,7 @@ def writeconfig(outpath, year, sector, filepatt, cq2gc): if os.path.exists(hcpath): break else: - raise IOError('file not found') + raise IOError(f'File not found: {hcpath}') with open(outpath, 'w') as hcf: print(sector, hcpatt, end='', flush=True) From cfd46a484e5bd58882ea9c40b3e2484d3147e58c Mon Sep 17 00:00:00 2001 From: roarkmd Date: Fri, 23 Jan 2026 14:18:29 -0500 Subject: [PATCH 03/10] Added scipy as a requirement. --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 8474e31..9710c3b 100644 --- a/setup.py +++ b/setup.py @@ -29,7 +29,7 @@ ], python_requires='>=3.6', install_requires=[ - "pyproj", "pandas", "xarray", "requests", "cmaqsatproc", "netcdf4" + "pyproj", "pandas", "xarray", "requests", "cmaqsatproc", "netcdf4", "scipy", ], extras_require={ 'aws': [ From d5ea559296be9b9feda48ce7397864b29fbf8114 Mon Sep 17 00:00:00 2001 From: roarkmd Date: Fri, 23 Jan 2026 14:20:09 -0500 Subject: [PATCH 04/10] Created finction to open files without specifying AWS bucket --- cmaq2hemco/utils.py | 34 +++++++++++++++++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/cmaq2hemco/utils.py b/cmaq2hemco/utils.py index b5c95f7..1581892 100644 --- a/cmaq2hemco/utils.py +++ b/cmaq2hemco/utils.py @@ -1,7 +1,7 @@ __all__ = [ 'plumerise_briggs', 'open_date', 'gd2hemco', 'pt2hemco', 'pt2gd', 'merge', 'to_ioapi', 'getmw', 'se_file', 'gd2matrix', 'gd2hemco_fast', - 'unitconvert', 'hemco_area', 'symlinks', 'gd_file' + 'unitconvert', 'hemco_area', 'symlinks', 'gd_file', 'open_file', ] import numpy as np @@ -227,6 +227,38 @@ def plumerise_briggs( return dzbriggs +def open_file( + date, tmpl +): + """ + Open all files for specific date + + Arguments + --------- + date : str + Date parsable by pandas.to_datetime + tmpl : str + strftime template for date file + (e.g., MCIP/12US1/GRIDCRO2D.12US1.35L.%y%m%d) + + Returns + ------- + ds : xarray.Dataset + File opened (either in memory or from disk) + """ + import pandas as pd + import io + import gzip + import cmaqsatproc as csp + global res + date = pd.to_datetime(date) + path = date.strftime(tmpl) + if path.endswith('.gz'): + bdy = io.BytesIO(gzip.open(path).read()) + f = csp.open_ioapi(bdy, engine='scipy') + else: + f = csp.open_ioapi(path, engine='scipy') + return f def open_date( date, tmpl, bucket, cache=True From a5bebfaf3b313d2072b228a1d689f4e0a37f1df0 Mon Sep 17 00:00:00 2001 From: roarkmd Date: Wed, 28 Jan 2026 00:45:55 -0500 Subject: [PATCH 05/10] Updated functions so that the speciation mechanism can be changed to a custom mechanism, other than cb6. Custom mechanisms can be uses as long as there is an appropriate mechanism molecular weight file listed in the running directory. --- cmaq2hemco/utils.py | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/cmaq2hemco/utils.py b/cmaq2hemco/utils.py index 1581892..7ba66b1 100644 --- a/cmaq2hemco/utils.py +++ b/cmaq2hemco/utils.py @@ -524,7 +524,7 @@ def gd2matrix(gf, elat, elon): return ol.set_index(['ROW', 'COL', 'lati', 'loni']) -def gd2hemco_fast(path, gf, elat, elon, verbose=0): +def gd2hemco_fast(path, gf, elat, elon, verbose=0, gc='cb6r5_ae7_aq', nr='cb6r5hap_ae7_aq'): """ Bilinear interpolation of fluxes (w/ MSFX2 factor) @@ -540,6 +540,10 @@ def gd2hemco_fast(path, gf, elat, elon, verbose=0): Edge longitudes for regular grid verbose : int Level of verbosity + gc : str + Name of gas-phase chemical mechanism. Used in getmw only. + nr : str + Name of non-reactive gas-phase mechanism (typically haps). Used in getmw on Arguments --------- @@ -579,14 +583,14 @@ def gd2hemco_fast(path, gf, elat, elon, verbose=0): tmp = (gf[dk] / qarea).interp(ROW=Y, COL=X) attrs = {k: v for k, v in gf[dk].attrs.items()} unit = attrs['units'].strip() + '/m**2' - tmp, unit = unitconvert(dk, tmp, unit=unit) + tmp, unit = unitconvert(dk, tmp, unit=unit, gc=gc, nr=nr) attrs['units'] = unit outf.addvar(dk, tmp.data, **attrs) return outf -def gd2hemco(path, gf, elat, elon, matrix=None, verbose=0): +def gd2hemco(path, gf, elat, elon, matrix=None, verbose=0, gc='cb6r5_ae7_aq', nr='cb6r5hap_ae7_aq'): """ Uses a fractional aera overlap interoplation. @@ -604,6 +608,10 @@ def gd2hemco(path, gf, elat, elon, matrix=None, verbose=0): fraction from row/col centroids to lat/lon centroids verbose : int Level of verbosity + gc : str + Name of gas-phase chemical mechanism. Used in getmw only. + nr : str + Name of non-reactive gas-phase mechanism (typically haps). Used in getmw on Arguments --------- @@ -694,14 +702,14 @@ def gd2hemco(path, gf, elat, elon, matrix=None, verbose=0): loni = gval.index.get_level_values('loni') tmp[ti, ki, lati, loni] += gval attrs = {k: v for k, v in gf[dk].attrs.items()} - tmp, unit = unitconvert(dk, tmp, attrs['units'], area=area) + tmp, unit = unitconvert(dk, tmp, attrs['units'], area=area, gc=gc, nr=nr) attrs['units'] = unit outf.addvar(dk, tmp, **attrs) return outf -def unitconvert(key, val, unit, area=None, inplace=True): +def unitconvert(key, val, unit, area=None, inplace=True, gc='cb6r5_ae7_aq', nr='cb6r5hap_ae7_aq'): """ Arguments --------- @@ -720,6 +728,11 @@ def unitconvert(key, val, unit, area=None, inplace=True): inplace : bool If True, do the unit conversion within the val array without allocating additional memory + gc : str + Name of gas-phase chemical mechanism. Used in getmw only. + nr : str + Name of non-reactive gas-phase mechanism (typically haps). Used in getmw only. + Returns ------- @@ -741,7 +754,7 @@ def unitconvert(key, val, unit, area=None, inplace=True): outunit.append('kg') elif unit in nps: try: - mw = getmw(key) + mw = getmw(key, gc=gc, nr=nr) factor *= mw outunit.append('kg') except KeyError as e: From bdccdddbbceb41838965218ed1cc1fae20ecb9d5 Mon Sep 17 00:00:00 2001 From: roarkmd Date: Wed, 28 Jan 2026 14:12:15 -0500 Subject: [PATCH 06/10] Added example running cmaq2hemco tool from files on disk (SMOKE) using GEOS Chem speciation mechanism. --- examples/smoke_geos_chem.py | 163 ++++++++++++++++++++++++++++++++++++ 1 file changed, 163 insertions(+) create mode 100755 examples/smoke_geos_chem.py diff --git a/examples/smoke_geos_chem.py b/examples/smoke_geos_chem.py new file mode 100755 index 0000000..17b84e2 --- /dev/null +++ b/examples/smoke_geos_chem.py @@ -0,0 +1,163 @@ +import pandas as pd +import numpy as np +import os +from cmaq2hemco.utils import pt2hemco, gd2hemco, gd2hemco_fast, gd_file, open_file +from cmaq2hemco.mechrc.geoschem_14_6_2 import writeconfig + +def geos_chem_mw(): + molecDct = {"CO": 28, "NH3": 17, "NO": 46, "NO2": 46, "SO2": 64, \ + "SULF": 98, "BC": 1, "OC": 1, "PMFINE": 1, "NIT": 1, \ + "SO4": 1, "PMC": 1, "HGIIGAS": 200.59,\ + "HGNRVA":200.59, "PHGI": 1, "DIESEL_PMEC": 1, \ + "DIESEL_PMFINE": 1, "DIESEL_PMNO3":1, "DIESEL_PMOC": 1, \ + "DIESEL_PMSO4": 1, "DIESEL_PMC": 1, "ACROLEIN": 56.0633, \ + "BUTADIENE13":54.0904, "VOC_INV": 1.0, \ + "CL2": 70.91,"CHROMHEX_C": 1, \ + "CHROMHEX_F": 1, "CHROMTRI_C": 1, "CHROMTRI_F": 1,\ + "NICKEL_C": 1, "NICKEL_F": 1, "BERYLLIUM_C": 1, \ + "BERYLLIUM_F": 1,"CADMIUM_C": 1, "CADMIUM_F": 1, "LEAD_C":1,\ + "LEAD_F": 1, "MANGANESE_C":1, "MANGANESE_F": 1, \ + "OXYL": 106.165, "PXYL": 106.165, "MXYL": 106.165,\ + "TOLU": 92.1384, "CL4_ETHE": 165.83, "TRIETHYLAMINE": 101.19,\ + "HEXAMETHY_DIIS": 168.2, "CHCL3": 119.3776, "CL_ETHE": 62.5,\ + "CL4_ETHANE1122": 167.85, "ETOX": 44.0526, "QUINOLINE": 129.16,\ + "ACRYLONITRILE": 53.06, "CL2_C2_12": 98.9592, \ + "BR2_C2_12": 187.86,"HYDRAZINE": 32.05, "CARBONTET": 153.8227,\ + "DICHLOROPROPENE": 110.97,"PROPDICHLORIDE": 112.9, \ + "MAL_ANHYDRIDE": 98.06, "DICHLOROBENZENE":147.0002, \ + "TOL_DIIS": 174.1561, "CL2_ME": 84.93, "CL3_ETHE": 131.3883,\ + "HCL": 36.46, "HONO": 46, "NOX": 46, "PM2_5": 1, "PM10": 1, "HFLUX": 1,\ + "NH3_FERT": 17, "PAL": 1, "PCA": 1, "PCL": 1, \ + "pFe": 1, "PH2O": 1, "PK": 1, "PMG": 1, "PMN": 1, "PMOTHR": 1, \ + "PNA": 1, "PNCOM": 1, "NH4": 1, "PSI": 1, "PTI": 1, + "ARSENIC_C": 1, "ARSENIC_F": 1, \ + "PAH_000E0": 379.00, "PAH_101E2": 268.00, "PAH_114E1": 256.00, + "PAH_176E2": 302.00, "PAH_176E3": 244.00, "PAH_176E4": 248.00, + "PAH_176E5": 228.00, "PAH_192E3": 278.00, "PAH_880E5": 196.00, + "ACETONITRILE": 41.05, "ACRYLICACID": 72.06, "ACRYLONITRILE": 53.06, + "CARBSULFIDE": 60.07, "CHLOROPRENE": 88.54, "ETHYLBENZ": 106.165, + "HEXANE": 86.175, "METHCHLORIDE": 50.49, "STYRENE": 104.15, + "XYLENES": 106.165, "NOX_INV": 46, "HF": 20.01, "NMOG": 1.0, + "VOC_BEIS": 46, "APIN": 136.234, "BPIN": 120, "SESQ": 180, "NR": 24, + "BENZOAPYRNE": 252.316, "TOG_INV": 1, "CO2": 44.01, "CO2_INV": 1, + "N2O_INV": 1, "CH4_INV": 16.042, + "ACET": 58.09, "ACR": 56.06, "ACTA": 60.06, "ALD2": 44.06, "ALK4": 58.12, + "ALK6": 100.2, "BALD": 106.12, "BENZ": 78.12, "C2H2": 26.05, "C2H4": 28.05, + "C2H6": 30.08, "C3H8": 44.11, "C4H6": 54.09, "CCl4": 153.82, "CH2Br2": 173.83, + "CH2Cl2": 84.93, "CH2O": 30.03, "CH4": 16.04, "CHCl3": 119.35, "CSL": 108.14, + "EBZ": 106.167, "EOH": 46.07, "FURA": 68.07, "GLYX": 58.04, "HCOOH": 46.03, + "ISOP": 68.13, "IVOC": 240.5, "MEK": 72.11, "MOH": 32.05, "MTPA": 136.26, + "MTPO": 136.26, "MVK": 70.09, "NAP": 128.18, "OCS": 60.07, "PHEN": 94.11, + "PRPE": 42.09, "RCHO": 58.09, "RCOOH": 74.09, "ROH": 60.11, "STYR": 104.1491, + "TMB": 106.167, "TOLU": 92.15, "UNK": 142.3, "UNR": 142.3, "XYLE": 106.18, + "CH3Br": 94.94, "CH3Cl": 50.45, "CH3I": 141.94, "HACTA": 76.0514, "LIMO": 136.26, + "MACR": 70.1, "MGLY": 72.07} + df = pd.DataFrame({"SPECIES": molecDct.keys(), "MOLWT": molecDct.values()}) + df.to_csv("cmaq_geoschem14_6_3_molwt.csv", index=False) + + + +def main(): + # Debug set to True will only process a single date + # or whatever is in the debug if satement below + debug = True + + # excluding burn and biogenics, and point fire sectors + # openburn + # beis4 + gkeys = """rwc + rail + onroad_gas + onroad_diesel + onroad_ca_adj_gas + onroad_ca_adj_diesel + np_solvents + np_oilgas + nonroad_gas + nonroad_diesel + nonpt + mexico_onroad + livestock + canada_ptdust_adj + canada_onroad + canada_og2D + canada_afdust_adj + airports + afdust_adj""".split() + + pkeys = """ptegu + canmex_point + cmv_c1c2_12 + cmv_c3_12 + pt_oilgas + ptnonipm""".split() + + #Generate speciation molecular weights + geos_chem_mw() + + # Date range to process emissions + dates = pd.date_range('2022-01-01', '2022-12-31', freq='d') + + # Define grid by edges + # Grid is latitude/longitude grid tenth degree grid + elat = np.linspace(15, 65, 501) + elon = np.linspace(-135, -50, 851) + + if debug: + dates = pd.to_datetime(['2022-01-04']) + # pkeys = pkeys[:1] + # gkeys = gkeys[:1] + print('**WARNING: in debug mode; only processing') + print(dates) + print(pkeys) + print(gkeys) + + for date in dates: + for gkey in gkeys: + outpath = ( + f'epa2022v2/{gkey}/{gkey}_{date:%Y-%m-%d}_geoschem14_6_3_2022he_geoschem.ncf' + ) + if os.path.exists(outpath): + continue + print(date, gkey) + os.makedirs(os.path.dirname(outpath), exist_ok=True) + try: + # open_gdemis could fail if date is not available (eg., only weekday) + gf = gd_file(open_file(date, tmpl='/home/mroark/git_repositories/cmaq2hemco/examples/rail/emis_mole_rail_%Y%m%d_12US1_geoschem14_6_3_2022he_geoschem.ncf')) + except Exception as e: + print(f'**WARNING:: Skipping {date} {gkey}: {e}') + continue + # using bilinear interpolation of fluxes + rgf = gd2hemco_fast(outpath, gf, elat, elon, gc='geoschem14_6_3', nr='geoschem14_6_3') + # use matrix interpolation for fractional area overlap (slow) + # rgf = gd2hemco(outpath, gf, elat, elon, matrix=matrix) + del rgf, gf + + for pkey in pkeys: + outpath = ( + f'epa2022v2/{pkey}/{pkey}_{date:%Y-%m-%d}_geoschem14_6_3_2022he_geoschem.ncf' + ) + if os.path.exists(outpath): + continue + print(date, pkey) + os.makedirs(os.path.dirname(outpath), exist_ok=True) + try: + # open_ptemis could fail if date is not available (eg., only weekday) + pf = gd_file(open_file(date, tmpl='/home/mroark/git_repositories/cmaq2hemco/examples/rail/emis_mole_rail_%Y%m%d_12US1_geoschem14_6_3_2022he_geoschem.ncf')) + #pf = gd_file(open_date(date, tmpl='sector/emis_mole_rail_%Y%m%d_12US1_geoschem14_6_3_2022he_geoschem.ncf', \ + # bucket='/home/mroark/git_repositories/cmaq2hemco/examples')) + except IOError as e: + print(f'**WARNING:: Skipping {date} {pkey}: {e}') + continue + rpf = pt2hemco(outpath, pf, elat, elon, gc='geoschem14_6_3', nr='geoschem14_6_3') # apply plume rise + del rpf, pf + + for sector in gkeys + pkeys: + hcpath = f'epa2022v2/{sector}/HEMCO_{sector}.rc' + secttmpl = f'epa2022v2/{sector}/{sector}_%Y-%m-%d_geoschem14_6_3_2022he_geoschem.ncf' + writeconfig(hcpath, 2022, sector, secttmpl) + + + +if __name__ == '__main__': + main() From a3353e3737163708b22c49be5a66a8cad3e7ab15 Mon Sep 17 00:00:00 2001 From: roarkmd Date: Wed, 4 Feb 2026 14:42:56 -0500 Subject: [PATCH 07/10] Created gd2hemco_fast_3D and added 3D capabilities to gd_file. Previously, gd_file only read the first layer of the file and silently dropped other layers if present. The mode parameter can be set to "3D" to not drop any files. gd2hemco_fast_3D is a modification of gd2hemco_fast that allows layers in 3D premerged files to undergo bilinear interpolation without dropping/ignoring layers above the first layer. For SMOKE to hemco usage, levels should be specified as the midpoints of the vertical layers. (SMOKE's default is to report VGLVLs as the endpoints (inclusive) of each model layer.) --- cmaq2hemco/utils.py | 134 ++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 128 insertions(+), 6 deletions(-) diff --git a/cmaq2hemco/utils.py b/cmaq2hemco/utils.py index 7ba66b1..85790c3 100644 --- a/cmaq2hemco/utils.py +++ b/cmaq2hemco/utils.py @@ -1,6 +1,6 @@ __all__ = [ 'plumerise_briggs', 'open_date', 'gd2hemco', 'pt2hemco', 'pt2gd', 'merge', - 'to_ioapi', 'getmw', 'se_file', 'gd2matrix', 'gd2hemco_fast', + 'to_ioapi', 'getmw', 'se_file', 'gd2matrix', 'gd2hemco_fast', 'gd2hemco_fast_3D', 'unitconvert', 'hemco_area', 'symlinks', 'gd_file', 'open_file', ] @@ -558,11 +558,13 @@ def gd2hemco_fast(path, gf, elat, elon, verbose=0, gc='cb6r5_ae7_aq', nr='cb6r5h ) clat = (elat[1:] + elat[:-1]) / 2 clon = (elon[1:] + elon[:-1]) / 2 - clev = _deflevs[:1] + clev = _deflevs[:1] #should be changed for all layers for 3D files + # Datakeys are pollutant variables datakeys = [ k for k in gf if k not in ('TFLAG', 'ti', 'ki', 'ri', 'ci', 'lon', 'lat') ] + outf = hemcofile( path, gf.time, clat, clon, lev=clev, varkeys=datakeys, attrs=gf.attrs ) @@ -708,6 +710,117 @@ def gd2hemco(path, gf, elat, elon, matrix=None, verbose=0, gc='cb6r5_ae7_aq', nr return outf +def gd2hemco_fast_3D(path, gf, elat, elon, lev=None, verbose=0, gc='cb6r5_ae7_aq', nr='cb6r5hap_ae7_aq'): + """ + Bilinear interpolation of fluxes to a regular lat/lon grid, preserving + multiple vertical layers. + + Arguments + --------- + path : str + Path to save the HEMCO file. + gf : xarray.Dataset + Gridded emissions (projected), typically with dims ('time','LAY','ROW','COL') + and having 'lon'/'lat' (in degrees) for the input grid and a 'crs' attribute. + elat : array-like + Edge latitudes for the destination regular grid (degrees_north). + elon : array-like + Edge longitudes for the destination regular grid (degrees_east). + gc : str + Name of gas-phase chemical mechanism. Used in getmw only. + nr : str + Name of non-reactive gas-phase mechanism (typically haps). Used in getmw on + lev : array-like or None + Vertical level coordinates to write to HEMCO. If None, uses _deflevs + truncated to the number of layers in gf (LAY). + verbose : int + Level of verbosity. + + Returns + ------- + outf : hemcofile + HEMCO file object (with .nc property). + """ + import pyproj + import numpy as np + import xarray as xr + + # Projected-area correction to convert totals to flux per m2 + proj = pyproj.Proj(gf.crs) + # areal_scale converts projected areas to true (on-sphere) areas + scale = proj.get_factors(gf['lon'], gf['lat']).areal_scale + qarea = gf.XCELL * gf.YCELL / scale # shape (ROW, COL) + + # Destination grid centers + clat = (elat[1:] + elat[:-1]) / 2 + clon = (elon[1:] + elon[:-1]) / 2 + + # Number of vertical layers (from input), and HEMCO levs to write + nk = gf.sizes['LAY'] if 'LAY' in gf.dims else 1 + if lev is None: + clev = _deflevs[:nk] + else: + clev = np.asarray(lev) + nk = clev.size + + # Variables to write (exclude helpers and coords) + datakeys = [ + k for k in gf + if k not in ('TFLAG', 'ti', 'ki', 'ri', 'ci', 'lon', 'lat') + ] + + # Initialize output file + outf = hemcofile( + path, gf.time, clat, clon, lev=clev, varkeys=datakeys, attrs=gf.attrs + ) + + # Add destination cell areas + area = hemco_area(elat, elon) + outf.addvar('AREA', area, units='m2', dims=('lat', 'lon')) + + # Build target coordinates in projected space for bilinear interpolation + LON, LAT = np.meshgrid(clon, clat) # destination centers (deg) + X, Y = proj(LON, LAT) # destination centers in proj + X = xr.DataArray(X, dims=('ROW', 'COL')) # use ROW/COL names for interp + Y = xr.DataArray(Y, dims=('ROW', 'COL')) + + # Prepare area factor as DataArray for safe broadcasting + qarea_da = xr.DataArray(qarea, dims=('ROW', 'COL')) + + for dk in datakeys: + dtot = float(gf[dk].sum()) + if dtot == 0: + if verbose > 0: + print(f'excluded {dk} (zero total)') + continue + if verbose > 0: + print(dk) + + da = gf[dk] + # Ensure we have a LAY dimension for consistency + if 'LAY' not in da.dims: + da = da.expand_dims(LAY=[0]) + + # Convert totals to flux per m^2 in the source grid (account MSFX2/area scaling) + flux = (da / qarea_da) + + # Bilinear interpolation to destination lat/lon grid for every time, layer + interp_flux = flux.interp(ROW=Y, COL=X) + + # Order to match HEMCO expected (time, lev, lat, lon) + interp_flux = interp_flux.transpose('time', 'LAY', 'ROW', 'COL') + + # Units: append /m**2 before converting so unitconvert won't divide by area again + attrs = {k: v for k, v in gf[dk].attrs.items()} + unit = attrs.get('units', 'unknown').strip() + unit = unit + '/m**2' + interp_flux, unit = unitconvert(dk, interp_flux, unit=unit, gc=gc, nr=nr) + attrs['units'] = unit + + # Write with explicit dims to match the 4D variable shape + outf.addvar(dk, interp_flux.data, dims=('time', 'lev', 'lat', 'lon'), **attrs) + + return outf def unitconvert(key, val, unit, area=None, inplace=True, gc='cb6r5_ae7_aq', nr='cb6r5hap_ae7_aq'): """ @@ -1232,7 +1345,7 @@ def symlinks(tmpl, dates, datetype=None, verbose=0): return links -def gd_file(ef): +def gd_file(ef, mode=None): """ Add lon/lat and rename TSTEP as time. @@ -1241,6 +1354,12 @@ def gd_file(ef): ef : xarray.Dataset ef must have the crs attribute and TSTEP coordinate variable + Optional Arguments + ------- + mode : None or "3D" + Used for files that have vertical layers. "3D" keeps all layers. + By default leyers are dropped. + Returns ------- gf : xarray.Dataset @@ -1248,9 +1367,12 @@ def gd_file(ef): removed. """ import pyproj - ef = ef.isel( - LAY=0, drop=True - ).rename(TSTEP='time') + if mode == '3D': + ef = ef.rename(TSTEP='time') + else: + ef = ef.isel( + LAY=0, drop=True + ).rename(TSTEP='time') proj = pyproj.Proj(ef.crs) Y, X = xr.broadcast(ef.ROW, ef.COL) LON, LAT = proj(X, Y, inverse=True) From 5330e06f42f4ea64158e5dc07d7ad34037fb6125 Mon Sep 17 00:00:00 2001 From: roarkmd Date: Thu, 19 Feb 2026 10:44:52 -0500 Subject: [PATCH 08/10] Added ability to open inline and stack groups files on disk. Added engine option to specify the backend used to open the netcdf file --- cmaq2hemco/utils.py | 49 +++++++++++++++++++++++++++++++++++++-------- 1 file changed, 41 insertions(+), 8 deletions(-) diff --git a/cmaq2hemco/utils.py b/cmaq2hemco/utils.py index 85790c3..77adfd7 100644 --- a/cmaq2hemco/utils.py +++ b/cmaq2hemco/utils.py @@ -1,7 +1,7 @@ __all__ = [ 'plumerise_briggs', 'open_date', 'gd2hemco', 'pt2hemco', 'pt2gd', 'merge', 'to_ioapi', 'getmw', 'se_file', 'gd2matrix', 'gd2hemco_fast', 'gd2hemco_fast_3D', - 'unitconvert', 'hemco_area', 'symlinks', 'gd_file', 'open_file', + 'unitconvert', 'hemco_area', 'symlinks', 'gd_file', 'open_gdfile', 'open_sefile', ] import numpy as np @@ -227,8 +227,8 @@ def plumerise_briggs( return dzbriggs -def open_file( - date, tmpl +def open_gdfile( + date, tmpl, engine='scipy' ): """ Open all files for specific date @@ -255,13 +255,46 @@ def open_file( path = date.strftime(tmpl) if path.endswith('.gz'): bdy = io.BytesIO(gzip.open(path).read()) - f = csp.open_ioapi(bdy, engine='scipy') + f = csp.open_ioapi(bdy, engine=engine) else: - f = csp.open_ioapi(path, engine='scipy') + f = csp.open_ioapi(path, engine=engine) return f +def open_sefile( + date, stack_groups_tmpl, inln_tmpl, engine='scipy' +): + """ + Open all files for specific date + + Arguments + --------- + date : str + Date parsable by pandas.to_datetime + stack_groups_tmpl : str + strftime template for date file + (e.g., MCIP/12US1/GRIDCRO2D.12US1.35L.%y%m%d) + inln_tmpl : str + strftime template for date file + (e.g., MCIP/12US1/GRIDCRO2D.12US1.35L.%y%m%d) + + Returns + ------- + ds : xarray.Dataset + File opened (either in memory or from disk) + """ + import pandas as pd + import cmaqsatproc as csp + date = pd.to_datetime(date) + stack_groups_path = date.strftime(stack_groups_tmpl) + inln_path = date.strftime(inln_tmpl) + sf = csp.open_ioapi(stack_groups_path, engine=engine) + ef = csp.open_ioapi(inln_path, engine=engine).isel(LAY=0, COL=0, drop=True) + ef = se_file(sf, ef) + return ef + + def open_date( - date, tmpl, bucket, cache=True + date, tmpl, bucket, cache=True, engine='scipy' ): """ Open all files for specific date @@ -307,7 +340,7 @@ def open_date( if dest.endswith('.gz'): bdy = io.BytesIO(gzip.open(dest).read()) - f = csp.open_ioapi(bdy, engine='scipy') + f = csp.open_ioapi(bdy, engine=engine) else: f = csp.open_ioapi(dest) else: @@ -320,7 +353,7 @@ def open_date( if path.endswith('.gz'): bdy = gzip.decompress() bdy = io.BytesIO(bdy) - f = csp.open_ioapi(bdy, engine='scipy') + f = csp.open_ioapi(bdy, engine=engine) return f From 1b5c92313dd447ac951f4565b22182d92eb1f116 Mon Sep 17 00:00:00 2001 From: roarkmd Date: Mon, 2 Mar 2026 16:52:43 -0500 Subject: [PATCH 09/10] Added gc and nr used to specify speciation mechanism molecular weights in the p22hemco function. Added saftey mechanism by removing default speciation mechanisms from unitconvert options. --- cmaq2hemco/utils.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/cmaq2hemco/utils.py b/cmaq2hemco/utils.py index 77adfd7..cd5fcf3 100644 --- a/cmaq2hemco/utils.py +++ b/cmaq2hemco/utils.py @@ -68,7 +68,7 @@ ) -def getmw(key, gc='cb6r5_ae7_aq', nr='cb6r5hap_ae7_aq'): +def getmw(key, gc=None, nr=None): """ Get the molecular weight (kg/mol) for a chemical mechanism species. The species may be an explicit or lumped species, so weights are mechanism @@ -359,7 +359,7 @@ def open_date( def pt2hemco( path, pf, elat, elon, ez=None, nk=11, temp_a=288.15, pres_a=101325, u=2.5, - verbose=0 + verbose=0, gc=None, nr=None ): """ Convert a point source file to a hemco-ready file @@ -478,7 +478,7 @@ def pt2hemco( tmp[vals.ti, vals.ki, vals.ri, vals.ci] = vals[dk].values attrs = {k: v for k, v in pf[dk].attrs.items()} unit = attrs['units'].strip() - tmp, unit = unitconvert(dk, tmp, unit, area=area) + tmp, unit = unitconvert(dk, tmp, unit, area=area, gc=gc, nr=nr) attrs['units'] = unit outf.addvar(dk, tmp, **attrs) @@ -557,7 +557,7 @@ def gd2matrix(gf, elat, elon): return ol.set_index(['ROW', 'COL', 'lati', 'loni']) -def gd2hemco_fast(path, gf, elat, elon, verbose=0, gc='cb6r5_ae7_aq', nr='cb6r5hap_ae7_aq'): +def gd2hemco_fast(path, gf, elat, elon, verbose=0, gc=None, nr=None): """ Bilinear interpolation of fluxes (w/ MSFX2 factor) @@ -625,7 +625,7 @@ def gd2hemco_fast(path, gf, elat, elon, verbose=0, gc='cb6r5_ae7_aq', nr='cb6r5h return outf -def gd2hemco(path, gf, elat, elon, matrix=None, verbose=0, gc='cb6r5_ae7_aq', nr='cb6r5hap_ae7_aq'): +def gd2hemco(path, gf, elat, elon, matrix=None, verbose=0, gc=None, nr=None): """ Uses a fractional aera overlap interoplation. @@ -743,7 +743,7 @@ def gd2hemco(path, gf, elat, elon, matrix=None, verbose=0, gc='cb6r5_ae7_aq', nr return outf -def gd2hemco_fast_3D(path, gf, elat, elon, lev=None, verbose=0, gc='cb6r5_ae7_aq', nr='cb6r5hap_ae7_aq'): +def gd2hemco_fast_3D(path, gf, elat, elon, lev=None, verbose=0, gc=None, nr=None): """ Bilinear interpolation of fluxes to a regular lat/lon grid, preserving multiple vertical layers. @@ -855,7 +855,7 @@ def gd2hemco_fast_3D(path, gf, elat, elon, lev=None, verbose=0, gc='cb6r5_ae7_aq return outf -def unitconvert(key, val, unit, area=None, inplace=True, gc='cb6r5_ae7_aq', nr='cb6r5hap_ae7_aq'): +def unitconvert(key, val, unit, area=None, inplace=True, gc=None, nr=None): """ Arguments --------- From 91d382f4267b32bce0aac18f6fc7e10cd86a6083 Mon Sep 17 00:00:00 2001 From: roarkmd Date: Fri, 27 Mar 2026 11:26:56 -0400 Subject: [PATCH 10/10] Changed processing of point files. Added function to remove pollutants with zero emissions. Added arguments to set variables in run script using argparse. --- examples/smoke_geos_chem.py | 169 +++++++++++++++++++----------------- 1 file changed, 91 insertions(+), 78 deletions(-) diff --git a/examples/smoke_geos_chem.py b/examples/smoke_geos_chem.py index 17b84e2..d22d51b 100755 --- a/examples/smoke_geos_chem.py +++ b/examples/smoke_geos_chem.py @@ -1,9 +1,13 @@ +#!/home/mroark/miniforge3/envs/cmaq2hemco/bin/python +import sys import pandas as pd import numpy as np import os -from cmaq2hemco.utils import pt2hemco, gd2hemco, gd2hemco_fast, gd_file, open_file +from cmaq2hemco.utils import gd2hemco_fast_3D, gd2hemco_fast, gd_file, open_gdfile, open_sefile, pt2hemco from cmaq2hemco.mechrc.geoschem_14_6_2 import writeconfig +import argparse +# Molecular weights for speciation mechanism def geos_chem_mw(): molecDct = {"CO": 28, "NH3": 17, "NO": 46, "NO2": 46, "SO2": 64, \ "SULF": 98, "BC": 1, "OC": 1, "PMFINE": 1, "NIT": 1, \ @@ -31,91 +35,84 @@ def geos_chem_mw(): "pFe": 1, "PH2O": 1, "PK": 1, "PMG": 1, "PMN": 1, "PMOTHR": 1, \ "PNA": 1, "PNCOM": 1, "NH4": 1, "PSI": 1, "PTI": 1, "ARSENIC_C": 1, "ARSENIC_F": 1, \ - "PAH_000E0": 379.00, "PAH_101E2": 268.00, "PAH_114E1": 256.00, - "PAH_176E2": 302.00, "PAH_176E3": 244.00, "PAH_176E4": 248.00, - "PAH_176E5": 228.00, "PAH_192E3": 278.00, "PAH_880E5": 196.00, - "ACETONITRILE": 41.05, "ACRYLICACID": 72.06, "ACRYLONITRILE": 53.06, - "CARBSULFIDE": 60.07, "CHLOROPRENE": 88.54, "ETHYLBENZ": 106.165, - "HEXANE": 86.175, "METHCHLORIDE": 50.49, "STYRENE": 104.15, - "XYLENES": 106.165, "NOX_INV": 46, "HF": 20.01, "NMOG": 1.0, - "VOC_BEIS": 46, "APIN": 136.234, "BPIN": 120, "SESQ": 180, "NR": 24, - "BENZOAPYRNE": 252.316, "TOG_INV": 1, "CO2": 44.01, "CO2_INV": 1, - "N2O_INV": 1, "CH4_INV": 16.042, - "ACET": 58.09, "ACR": 56.06, "ACTA": 60.06, "ALD2": 44.06, "ALK4": 58.12, - "ALK6": 100.2, "BALD": 106.12, "BENZ": 78.12, "C2H2": 26.05, "C2H4": 28.05, - "C2H6": 30.08, "C3H8": 44.11, "C4H6": 54.09, "CCl4": 153.82, "CH2Br2": 173.83, - "CH2Cl2": 84.93, "CH2O": 30.03, "CH4": 16.04, "CHCl3": 119.35, "CSL": 108.14, - "EBZ": 106.167, "EOH": 46.07, "FURA": 68.07, "GLYX": 58.04, "HCOOH": 46.03, - "ISOP": 68.13, "IVOC": 240.5, "MEK": 72.11, "MOH": 32.05, "MTPA": 136.26, - "MTPO": 136.26, "MVK": 70.09, "NAP": 128.18, "OCS": 60.07, "PHEN": 94.11, - "PRPE": 42.09, "RCHO": 58.09, "RCOOH": 74.09, "ROH": 60.11, "STYR": 104.1491, - "TMB": 106.167, "TOLU": 92.15, "UNK": 142.3, "UNR": 142.3, "XYLE": 106.18, - "CH3Br": 94.94, "CH3Cl": 50.45, "CH3I": 141.94, "HACTA": 76.0514, "LIMO": 136.26, + "PAH_000E0": 379.00, "PAH_101E2": 268.00, "PAH_114E1": 256.00, \ + "PAH_176E2": 302.00, "PAH_176E3": 244.00, "PAH_176E4": 248.00, \ + "PAH_176E5": 228.00, "PAH_192E3": 278.00, "PAH_880E5": 196.00, \ + "ACETONITRILE": 41.05, "ACRYLICACID": 72.06, "ACRYLONITRILE": 53.06, \ + "CARBSULFIDE": 60.07, "CHLOROPRENE": 88.54, "ETHYLBENZ": 106.165, \ + "HEXANE": 86.175, "METHCHLORIDE": 50.49, "STYRENE": 104.15, \ + "XYLENES": 106.165, "NOX_INV": 46, "HF": 20.01, "NMOG": 1.0, \ + "VOC_BEIS": 46, "APIN": 136.234, "BPIN": 120, "SESQ": 180, "NR": 24, \ + "BENZOAPYRNE": 252.316, "TOG_INV": 1, "CO2": 44.01, "CO2_INV": 1, \ + "N2O_INV": 1, "CH4_INV": 16.042, \ + "ACET": 58.09, "ACR": 56.06, "ACTA": 60.06, "ALD2": 44.06, "ALK4": 58.12, \ + "ALK6": 100.2, "BALD": 106.12, "BENZ": 78.12, "C2H2": 26.05, "C2H4": 28.05, \ + "C2H6": 30.08, "C3H8": 44.11, "C4H6": 54.09, "CCl4": 153.82, "CH2Br2": 173.83, \ + "CH2Cl2": 84.93, "CH2O": 30.03, "CH4": 16.04, "CHCl3": 119.35, "CSL": 108.14, \ + "EBZ": 106.167, "EOH": 46.07, "FURA": 68.07, "GLYX": 58.04, "HCOOH": 46.03, \ + "ISOP": 68.13, "IVOC": 240.5, "MEK": 72.11, "MOH": 32.05, "MTPA": 136.26, \ + "MTPO": 136.26, "MVK": 70.09, "NAP": 128.18, "OCS": 60.07, "PHEN": 94.11, \ + "PRPE": 42.09, "RCHO": 58.09, "RCOOH": 74.09, "ROH": 60.11, "STYR": 104.1491, \ + "TMB": 106.167, "TOLU": 92.15, "UNK": 142.3, "UNR": 142.3, "XYLE": 106.18, \ + "CH3Br": 94.94, "CH3Cl": 50.45, "CH3I": 141.94, "HACTA": 76.0514, "LIMO": 136.26, \ "MACR": 70.1, "MGLY": 72.07} df = pd.DataFrame({"SPECIES": molecDct.keys(), "MOLWT": molecDct.values()}) df.to_csv("cmaq_geoschem14_6_3_molwt.csv", index=False) +# Function to drop pollutants that have zero emissions +# Consider turnning off if merging sectors together post-cmaq2hemco tool +def drop_zero_values(ds): + exclude_vars = ['lon', 'lat', 'ISTACK', 'STKDM', 'STKHT', 'STKTK', 'STKVE', 'STKFLW', 'STKCNT', 'XLOCA', 'YLOCA', 'IFIP', 'LMAJOR', 'LPING'] + data_variables = [ var for var in list(ds.data_vars) if var not in exclude_vars] + for var in data_variables: + if (ds[var].values.sum() == 0.0): + ds = ds.drop_vars(var) + # Get list of non-zero vars + data_variables = [ var for var in list(ds.data_vars) if var not in exclude_vars] + ds = ds.assign_attrs({'VAR-LIST': ''.join([f"{var:<16}" for var in data_variables]), 'NVARS': np.int32(len(data_variables))}) + return ds def main(): - # Debug set to True will only process a single date - # or whatever is in the debug if satement below - debug = True - - # excluding burn and biogenics, and point fire sectors - # openburn - # beis4 - gkeys = """rwc - rail - onroad_gas - onroad_diesel - onroad_ca_adj_gas - onroad_ca_adj_diesel - np_solvents - np_oilgas - nonroad_gas - nonroad_diesel - nonpt - mexico_onroad - livestock - canada_ptdust_adj - canada_onroad - canada_og2D - canada_afdust_adj - airports - afdust_adj""".split() - - pkeys = """ptegu - canmex_point - cmv_c1c2_12 - cmv_c3_12 - pt_oilgas - ptnonipm""".split() + # parse arguments + parser = argparse.ArgumentParser(description="Converts SMOKE premerged or inln/stack groups files to GEOS-Chem style file") + parser.add_argument("--sector", type=str, help="Sector to be processed") + parser.add_argument("--filetype", type=str, help="Either 2D 'premerged' or 'inln'") + parser.add_argument("--startdate", type=str, help="Starting date in YYYY-MM-DD format") + parser.add_argument("--enddate", type=str, help="Ending date in YYYY-MM-DD format") + parser.add_argument("--outdir", type=str, help="Output directory of files. Should be the directory above the sector specific directory.") + parser.add_argument("--indir", type=str, help="Input directory of sector files. Should be the directory above the sector specific directory.") + + args = parser.parse_args() + + if args.filetype == 'premerged': + gkeys = [args.sector] + pkeys = list() + elif args.filetype == 'inln': + gkeys = list() + pkeys = [args.sector] + else: + raise ValueError(f"Filetype not recognized: {args.filetype}") + + + outdir = args.outdir + indir = args.indir #Generate speciation molecular weights geos_chem_mw() # Date range to process emissions - dates = pd.date_range('2022-01-01', '2022-12-31', freq='d') + dates = pd.date_range(args.startdate, args.enddate, freq='d') # Define grid by edges # Grid is latitude/longitude grid tenth degree grid elat = np.linspace(15, 65, 501) elon = np.linspace(-135, -50, 851) - if debug: - dates = pd.to_datetime(['2022-01-04']) - # pkeys = pkeys[:1] - # gkeys = gkeys[:1] - print('**WARNING: in debug mode; only processing') - print(dates) - print(pkeys) - print(gkeys) - for date in dates: for gkey in gkeys: outpath = ( - f'epa2022v2/{gkey}/{gkey}_{date:%Y-%m-%d}_geoschem14_6_3_2022he_geoschem.ncf' + f'{outdir}/{gkey}/{gkey}_{date:%Y-%m-%d}_epa2022v2_he_22m.ncf' ) if os.path.exists(outpath): continue @@ -123,19 +120,25 @@ def main(): os.makedirs(os.path.dirname(outpath), exist_ok=True) try: # open_gdemis could fail if date is not available (eg., only weekday) - gf = gd_file(open_file(date, tmpl='/home/mroark/git_repositories/cmaq2hemco/examples/rail/emis_mole_rail_%Y%m%d_12US1_geoschem14_6_3_2022he_geoschem.ncf')) + gf = gd_file(open_gdfile(date, tmpl=f'{indir}/premerged/{gkey}/emis_mole_{gkey}_%Y%m%d_12US1_geoschem14_6_3_2022he_geoschem.ncf', + engine='scipy')) except Exception as e: - print(f'**WARNING:: Skipping {date} {gkey}: {e}') - continue + try: + # open_gdemis could fail if date is not available (eg., only weekday) + gf = gd_file(open_gdfile(date, tmpl=f'{indir}/premerged/{gkey}/emis_mole_{gkey}_%Y%m%d_12US1_geoschem14_6_3_2022he_geoschem.ncf.gz', + engine='scipy')) + except Exception as e: + print(f'**WARNING:: Skipping {date} {gkey}: {e}') + continue # using bilinear interpolation of fluxes + gf = drop_zero_values(gf) rgf = gd2hemco_fast(outpath, gf, elat, elon, gc='geoschem14_6_3', nr='geoschem14_6_3') # use matrix interpolation for fractional area overlap (slow) # rgf = gd2hemco(outpath, gf, elat, elon, matrix=matrix) del rgf, gf - for pkey in pkeys: outpath = ( - f'epa2022v2/{pkey}/{pkey}_{date:%Y-%m-%d}_geoschem14_6_3_2022he_geoschem.ncf' + f'{outdir}/{pkey}/{pkey}_{date:%Y-%m-%d}_epa2022v2_he_22m.ncf' ) if os.path.exists(outpath): continue @@ -143,18 +146,28 @@ def main(): os.makedirs(os.path.dirname(outpath), exist_ok=True) try: # open_ptemis could fail if date is not available (eg., only weekday) - pf = gd_file(open_file(date, tmpl='/home/mroark/git_repositories/cmaq2hemco/examples/rail/emis_mole_rail_%Y%m%d_12US1_geoschem14_6_3_2022he_geoschem.ncf')) - #pf = gd_file(open_date(date, tmpl='sector/emis_mole_rail_%Y%m%d_12US1_geoschem14_6_3_2022he_geoschem.ncf', \ - # bucket='/home/mroark/git_repositories/cmaq2hemco/examples')) - except IOError as e: - print(f'**WARNING:: Skipping {date} {pkey}: {e}') - continue + pf = open_sefile(date, \ + stack_groups_tmpl=f'{indir}/smoke_out/{pkey}/stack_groups_{pkey}_{date:%Y%m%d}_12US1_2022he_geoschem.ncf', \ + inln_tmpl=f'{indir}/smoke_out/{pkey}/inln_mole_{pkey}_{date:%Y%m%d}_12US1_geoschem14_6_3_2022he_geoschem.ncf', + engine='scipy') + except Exception as e: + print(e) + try: + pf = open_sefile(date, \ + stack_groups_tmpl=f'{indir}/smoke_out/{pkey}/stack_groups_{pkey}_12US1_2022he_geoschem.ncf', \ + inln_tmpl=f'{indir}/smoke_out/{pkey}/inln_mole_{pkey}_{date:%Y%m%d}_12US1_geoschem14_6_3_2022he_geoschem.ncf', + engine='scipy') + except IOError as e: + print(f'**WARNING:: Skipping {date} {pkey}: {e}') + continue + pf = drop_zero_values(pf) rpf = pt2hemco(outpath, pf, elat, elon, gc='geoschem14_6_3', nr='geoschem14_6_3') # apply plume rise del rpf, pf + for sector in gkeys + pkeys: - hcpath = f'epa2022v2/{sector}/HEMCO_{sector}.rc' - secttmpl = f'epa2022v2/{sector}/{sector}_%Y-%m-%d_geoschem14_6_3_2022he_geoschem.ncf' + hcpath = f'{outdir}/{sector}/HEMCO_{sector}.rc' + secttmpl = f'{outdir}/{sector}/{sector}_%Y-%m-%d_epa2022v2_he_22m.ncf' writeconfig(hcpath, 2022, sector, secttmpl)