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/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) 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) diff --git a/cmaq2hemco/utils.py b/cmaq2hemco/utils.py index b5c95f7..cd5fcf3 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' + 'to_ioapi', 'getmw', 'se_file', 'gd2matrix', 'gd2hemco_fast', 'gd2hemco_fast_3D', + 'unitconvert', 'hemco_area', 'symlinks', 'gd_file', 'open_gdfile', 'open_sefile', ] import numpy as np @@ -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 @@ -227,9 +227,74 @@ def plumerise_briggs( return dzbriggs +def open_gdfile( + date, tmpl, engine='scipy' +): + """ + 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=engine) + else: + 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 @@ -275,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: @@ -288,13 +353,13 @@ 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 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 @@ -413,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) @@ -492,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): +def gd2hemco_fast(path, gf, elat, elon, verbose=0, gc=None, nr=None): """ Bilinear interpolation of fluxes (w/ MSFX2 factor) @@ -508,6 +573,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 --------- @@ -522,11 +591,13 @@ def gd2hemco_fast(path, gf, elat, elon, verbose=0): ) 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 ) @@ -547,14 +618,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=None, nr=None): """ Uses a fractional aera overlap interoplation. @@ -572,6 +643,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 --------- @@ -662,14 +737,125 @@ 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 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. + + 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) -def unitconvert(key, val, unit, area=None, inplace=True): + # 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=None, nr=None): """ Arguments --------- @@ -688,6 +874,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 ------- @@ -709,7 +900,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: @@ -1187,7 +1378,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. @@ -1196,6 +1387,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 @@ -1203,9 +1400,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) diff --git a/examples/smoke_geos_chem.py b/examples/smoke_geos_chem.py new file mode 100755 index 0000000..d22d51b --- /dev/null +++ b/examples/smoke_geos_chem.py @@ -0,0 +1,176 @@ +#!/home/mroark/miniforge3/envs/cmaq2hemco/bin/python +import sys +import pandas as pd +import numpy as np +import os +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, \ + "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) + +# 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(): + # 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(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) + + for date in dates: + for gkey in gkeys: + outpath = ( + f'{outdir}/{gkey}/{gkey}_{date:%Y-%m-%d}_epa2022v2_he_22m.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_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: + 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'{outdir}/{pkey}/{pkey}_{date:%Y-%m-%d}_epa2022v2_he_22m.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 = 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'{outdir}/{sector}/HEMCO_{sector}.rc' + secttmpl = f'{outdir}/{sector}/{sector}_%Y-%m-%d_epa2022v2_he_22m.ncf' + writeconfig(hcpath, 2022, sector, secttmpl) + + + +if __name__ == '__main__': + main() 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': [