CMAQ to GeoTIFF

cmaq2geotiff.py

"""
# Convert a CMAQ IOAPI File to a GeoTiff

---
    author: Barron H. Henderson
    last updated: 2025-09-05
---

Demonstrates converting a CMAQ file to a GeoTiff. GeoTiff is useful for ArcGIS or QGIS.

Requries numpy, xarray, and rasterio
"""
import numpy as np
import xarray as xr
import argparse
import os

prsr = argparse.ArgumentParser()
prsr.add_argument('-O', '--overwrite', default=False, action='store_true')
prsr.add_argument('cmaqpath')
prsr.add_argument('gtifpath')

args = prsr.parse_args()
if os.path.exists(args.gtifpath) and not args.overwrite:
    print(f'Keeping existing {args.gtifpath}')
    exit(1)

f = xr.open_dataset(args.cmaqpath)
f.coords['COL'] = f.XORIG + np.arange(0.5, f.NCOLS) * f.XCELL
f.coords['ROW'] = f.YORIG + np.arange(0.5, f.NROWS) * f.YCELL
proj4str = f'+proj=lcc +lat_1={f.P_ALP} +lat_2={f.P_BET} +lat_0={f.YCENT} +lon_0={f.XCENT} +R=6370000.0 +no_defs'
# Choose just the variable you want to output
# Sum all layers (this time just one)
# Average all times (this time just one)
# Rename the x/y dimensions
tmpf = f[['NOX']].mean('TSTEP').sum('LAY').rename(COL='x', ROW='y')
# Add the projection information in RIO conventions
tmpf.rio.write_crs(proj4str, inplace=True)
# Save out result.
tmpf.rio.to_raster(args.gtifpath)