Inventory Analysis and Visualization

To use this notebook, the GRA2PES dataset must be downloaded and regridded using the process outlined in https://github.com/agmeyer4/atmos/tree/merge-atmos-column/workflows/gra2pes. There is also data included from the EPA and EDGAR datasets which must be downloaded and linked independently.

#Import libraries
import os
import xarray as xr
import numpy as np
import pandas as pd
import xesmf as xe
import sys 
from plotly import graph_objects as go
from plotly.subplots import make_subplots
from cartopy.io.img_tiles import OSM
import calendar
import matplotlib.pyplot as plt
import pyproj
import matplotlib
import datetime
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib.ticker import ScalarFormatter, FuncFormatter
import seaborn as sns
import cartopy.feature as cfeature
import cartopy.io.img_tiles as cimgt
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
sys.path.append('/uufs/chpc.utah.edu/common/home/u0890904/LAIR_1/atmos')
from configs.gra2pes import gra2pes_config
from workflows.gra2pes import gra2pes_regrid
from utils import gen_utils, datetime_utils, xr_utils, plot_utils, gra2pes_utils

#Autoreload changed local modules
%load_ext autoreload
%autoreload 2

#Define the paths
project_path = '/uufs/chpc.utah.edu/common/home/u0890904/LAIR_1/Projects/Ratio_paper_2025_v5'
project_figures_path = os.path.join(project_path,'Figures')
project_data_path = os.path.join(project_path,'Data')
project_pickle_path = os.path.join(project_path,'Data','Pickled')
inventory_figures_path = os.path.join(project_figures_path,'inventories')

wbb_loc = {'lat':40.768,'lon':-111.854}
uusyr_loc = {'lat':41.089,'lon':-112.119}

molar_masses = {'CO2':44.01, 'CO':28.01, 'CH4':16.04}

slv_extent = {'lon_min':-112.15,
                  'lon_max':-111.7,
                  'lat_min':40.351,
                  'lat_max':41.2} 

map_extent={'lon_min':-112.22,
            'lon_max':-111.65,
            'lat_min':40.35,
            'lat_max':41.3} 

specs = ['CO2','CO','CH4']

# Define point sources
point_sources = {

    'landfills' :   {   'locs':[(-112.043, 40.746),
                             (-112.055, 40.558),
                             (-111.917, 40.911)],
                        'plot_details':{'marker':'^','color':'red','edgecolor':'black','s':150}
                    },

    'refineries': {     'locs':[(-111.924, 40.825),
                                (-111.920, 40.838),
                                (-111.909, 40.794),
                                (-111.904, 40.887),
                                (-111.910, 40.868)],
                        'plot_details':{'marker':'o','s':150,'color':'red','edgecolor':'black'}
                    },
    'wastewater' : {    'locs':[(-112.074, 40.728),
                                (-111.900, 40.682),
                                (-111.919, 40.504),
                                (-111.931, 40.813),
                                (-111.924, 40.615),
                                (-111.942, 40.842),
                                (-111.933, 40.903),
                                (-111.945, 41.001)],
                        'plot_details':{'marker':'^','color':'green'}
                    },

}

Load Data

GRA2PES

# Load the GRA2PES data organized using atmos.inventories.gra2pes

config = gra2pes_config.Gra2pesConfig('/uufs/chpc.utah.edu/common/home/u0890904/LAIR_1/atmos/configs/gra2pes/gra2pes_config_1.yaml')
regrid_id = '0.025x0.025'
dtr = datetime_utils.DateTimeRange('2021-01-01','2021-12-31 23:59',tz = 'UTC')
rgh = gra2pes_utils.RegriddedGra2pesHandler(config,regrid_id)
regridded_ds = rgh.open_ds_inrange(dtr,sectors='all',chunks = {'utc_hour':1,'month':1,'lat':100,'lon':100})
grapes_ds = rgh.rework_ds_dt(regridded_ds[['CO2','CO','HC01']])
grapes_ds = grapes_ds.rename({'HC01':'CH4'})
grapes_gca = xr.open_dataset(os.path.join(rgh.regridded_path,'details','grid_out_area.nc'))['cell_area']

# Trim to the slv extent
grapes_slv_ds = xr_utils.trim_ds_to_extent(grapes_ds,slv_extent)
grapes_slv_gca = xr_utils.trim_ds_to_extent(grapes_gca,slv_extent)

# Convert Units
unit_converter = xr_utils.UnitConverter()
grapes_slv_ds = unit_converter.convert(grapes_slv_ds, 'km^-2', 'm^-2')
grapes_slv_ds = unit_converter.convert_mole_to_g(grapes_slv_ds, molar_masses)
grapes_slv_ds = unit_converter.convert(grapes_slv_ds, 'g', 'metric_Ton')
absolute_grapes_slv_ds = unit_converter.convert_m2_to_gridcell(grapes_slv_ds, grapes_slv_gca)

# Coarsen the data to 0.05x0.05 to match the resolution of GRA2PES in LCC
absolute_grapes_slv_ds_coarse = absolute_grapes_slv_ds.coarsen(lat = 2, lon = 2, boundary = 'trim').sum()

EDGAR

# Define a loading function for the EDGAR data

def load_edgar_nc(edgar_path,spec,emi_or_flx,year):
    if spec in ['CH4','CO2']:
        nc_folder_name = f'v8.0_FT2022_GHG_{spec}_{year}_TOTALS_{emi_or_flx}_nc'
        nc_file_name = '_'.join(nc_folder_name.split('_')[:-1]) + '.nc'
    elif spec == 'CO':
        nc_folder_name = f'EDGARv6.1_{spec}_{year}_TOTALS.0.1x0.1'
        nc_file_name = (nc_folder_name) + '.nc'

    nc_fullpath = os.path.join(edgar_path,spec,emi_or_flx,nc_folder_name,nc_file_name)
    try:    
        ds = xr.open_dataset(nc_fullpath)
    except FileNotFoundError as f:
        print(f)
        raise FileNotFoundError(f'There doesnt seem to be data for this configuration:{spec},{emi_or_flx},{year}')
    if spec == 'CO':
        ds = coord_change_360_180(ds)
        ds = ds.rename({'emi_co':'fluxes'})
        ds['fluxes'].attrs['substance'] = 'CO'

    gca_fname = nc_file_name.split('.nc')[0] + '_gca' + '.nc'
    try:
        gca_fullpath = os.path.join(edgar_path,spec,emi_or_flx,nc_folder_name,gca_fname)
        gca_ds = xr.open_dataset(gca_fullpath)
        if spec == 'CO':
            gca_ds = coord_change_360_180(gca_ds)
        ds = add_gca_to_ds(ds,gca_ds)
    except FileNotFoundError:
        print(f'No grid cell area file {gca_fullpath} found. Create with cdo to use. Returning ds without gca')
        return ds

    return ds

def coord_change_360_180(ds):
    return ds.assign_coords(lon=(((ds.lon + 180) % 360) - 180)).sortby('lon')

def add_gca_to_ds(orig_ds,gca_ds):
    orig_ds['grid_cell_area'] = gca_ds['cell_area']
    return orig_ds

#Load Edgar inventories
edgar_path = '/uufs/chpc.utah.edu/common/home/lin-group9/agm/inventories/EDGAR/'
emi_or_flx = 'flx'

# For each species, load the dataset and process
edgar_dss = {
    'CO2': {'year':2021,'ds':None},
    'CH4': {'year':2021,'ds':None},
    'CO': {'year':2018,'ds':None}
}

for spec, info in edgar_dss.items():
    ds = load_edgar_nc(edgar_path,spec,emi_or_flx,info['year'])
    ds['emissions'] = (ds['grid_cell_area'] * ds['fluxes']) * 86400 * 365 * 1E-3  # Emissions in metric tons per year
    ds['emissions'].attrs = {'substance': spec, 'units': 'metric_Ton gridcell-1 year-1'}
    edgar_dss[spec]['ds'] = xr_utils.trim_ds_to_extent(ds,slv_extent)

EPA Massakers

Flux values in original nc files are $Flux_{ch4}$ in $[{molec_{CH4}}/{cm^2s}]$.

Convert flux values to emissions values by $Emission_{ch4} = Flux_{ch4} * grid cell area$
Where $Emission_{ch4}$ in $[molec_{ch4}/s]$

CONVERSION to g/yr

$Emissions_{ch4} [g/year] = Emissions_{ch4} [molec_{ch4}/s] * Conversion factor $

$Conversion factor = 1/6.022x10^23 [mole/molec_{ch4}] * 16.04 [g_{ch4}/mole] * 86400 [s/day] * y [days/year] $

The number of days per year (y) depends on whether the given year is a leap year.

# Define a loading function for the EPA data
def load_epa_nc(path,year):
    fname = f'Express_Extension_Gridded_GHGI_Methane_v2_{year}.nc'
    ds = xr.open_dataset(os.path.join(path,fname))
    name_map = {k: k.replace('emi','flx') for k in ds} #in the nc file, the data variables are names "emi..." which is confusing as they are fluxes. so rename to flx...
    ds = ds.rename(name_map)
    return ds

#Load the Maasakkers data
epa_path = '/uufs/chpc.utah.edu/common/home/lin-group9/agm/inventories/Maasakkers'
year = 2020

epa_full = load_epa_nc(epa_path,year).sum(dim='time',keep_attrs = True)
epa_slv = xr_utils.trim_ds_to_extent(epa_full,slv_extent)

#Sum variables
data_variables = list(epa_slv.keys()) #all of the vars in the nc file
emi_variables = [var for var in data_variables if var != 'grid_cell_area'] #drop the grid cell area from the variables list so we can sum
og_units = epa_slv[emi_variables[0]].attrs['units']
for emi_var in emi_variables: #check to make sure the units are the same on all of the emissions variables
    if epa_slv[emi_var].attrs['units'] != og_units:
        raise ValueError(f'The units are different in {emi_var}, so you shouldnt sum')
epa_slv['flx_ch4_sum'] = epa_slv[emi_variables].to_array().sum('variable') #sum along the emissions variables
epa_slv['flx_ch4_sum'].attrs = {'units':og_units} #add the units

#convert flux units to emissions units
# emissions[tonne/yr] = grid_cell_area[cm2] * flx_ch4[molec/cm2/s] * (1/6.022E23)[mol/molec] * 16.04[gch4/mol] * 86400[s/d] * 365[d/y] * 1E-6[tonne/g]
epa_slv['emissions'] = (epa_slv['grid_cell_area']*epa_slv['flx_ch4_sum'])*(1/6.022E23)*16.04*86400*365*1E-6  #multiply the sums by the grid cell area
epa_slv['emissions'].attrs['units'] = 'metric_Ton gridcell-1 year-1'

absolute_epa_slv_yrsum = epa_slv

Compute and Plot Yearly totals

# Aggregate to yearly emissions

absolute_grapes_slv_yrsum_coarse = absolute_grapes_slv_ds_coarse.sel(sector='total').sum(dim='datetime',keep_attrs = True)
absolute_grapes_slv_yrsum_coarse = unit_converter.update_all_units(absolute_grapes_slv_yrsum_coarse, 'hr^-1','yr^-1')

#CO2 inventory plot
fig_id = f'2021_gra2pes_co2_fullyear_inventory'
savefig = True
showfig = True

plot_ds = absolute_grapes_slv_yrsum_coarse
spec = 'CO2'
cmap = 'Blues'

fig = plt.figure(figsize = (15,15))
labsize = 44
proj = ccrs.PlateCarree()

ax = plt.axes(projection = proj)
ax.set_extent([map_extent['lon_min'],map_extent['lon_max'],map_extent['lat_min'],map_extent['lat_max']],crs=proj)
request = cimgt.GoogleTiles(style='satellite')
scale = 12.0 # prob have to adjust this
ax.add_image(request,int(scale))

max_value = plot_ds[spec].max().values
vmax_value = max_value*0.6
map = plot_ds[spec].plot.pcolormesh('lon','lat',ax = ax,alpha=0.6,cmap=cmap,add_colorbar=False,edgecolors = (0,0,0,0.1),linewidth = 0.2,vmax = vmax_value)

ax.scatter(wbb_loc['lon'],wbb_loc['lat'],c = 'white',marker='X',s=500,zorder = 10,label = "WBB",edgecolor = 'k',linewidth = 2)

cbformat = matplotlib.ticker.ScalarFormatter()   # create the formatter
cbformat.set_powerlimits((-2,2)) 
cbar = plt.colorbar(map,fraction=0.05,pad = 0.02,format=cbformat)
exp = gen_utils.calculate_exponent(max_value)
cbar.set_label(label = f"$CO_{2}$ emissions (Tonne$\\times 10^{{{exp}}}$)",size = labsize)
cbar.ax.tick_params(labelsize = labsize)


emissions_sum = float(plot_ds[spec].sum())
exp= gen_utils.calculate_exponent(emissions_sum)
sn_es = round(emissions_sum/10**(exp),2)
text = f"$CO_{2}$ sum ={sn_es}$\\times 10^{{{exp}}}$ Tonne"

t1 = fig.text(0.275,0.92,text,fontsize = labsize)
t1.set_bbox(dict(facecolor = (1,1,1,0.7)))
#plt.title("2021 GRA2PES $CO_{2}$ Inventory",fontsize = labsize*1.2)
fig.tight_layout()
ax.set_title("")


# Explicitly set longitude and latitude ticks
xticks = [-112.1,-111.9,-111.7]
yticks = [40.4,40.6,40.8,41.0,41.2]
ax.set_xticks(xticks, crs=proj)
ax.set_yticks(yticks, crs=proj)
lon_formatter = LongitudeFormatter()
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_yticklabels([f"{y:.1f}°" for y in yticks], fontsize=labsize/2)
ax.set_xticklabels([f"{x:.1f}°" for x in xticks], fontsize=labsize/2)
ax.tick_params(
    left=False, labelleft=True,   # hide left labels
    top=False, labeltop=False,     # hide top labels
    right=False, labelright=False,   # show right labels
    bottom=False, labelbottom=True, # show bottom labels
    length=0  # remove the tick marks themselves
)
ax.set_xlabel('')  # remove x-axis label
ax.set_ylabel('')  # remove y-axis label


if savefig:
    fig_name = f'{fig_id}.png'
    fig.savefig(os.path.join(inventory_figures_path,fig_name),dpi=500,bbox_inches = "tight")
if showfig:
    plt.show()
else:
    plt.close()

#CH4 inventory plot
fig_id = f'2021_gra2pes_ch4_fullyear_inventory'
savefig = True
showfig = True

plot_ds = absolute_grapes_slv_yrsum_coarse

spec = 'CH4'
cmap = 'Reds'

fig = plt.figure(figsize = (15,15))
labsize = 44
proj = ccrs.PlateCarree()

ax = plt.axes(projection = proj)
ax.set_extent([map_extent['lon_min'],map_extent['lon_max'],map_extent['lat_min'],map_extent['lat_max']],crs=proj)
request = cimgt.GoogleTiles(style='satellite')
scale = 12.0 # prob have to adjust this
ax.add_image(request,int(scale))


max_value = plot_ds[spec].max().values
vmax_value = max_value*0.9
map = plot_ds[spec].plot.pcolormesh('lon','lat',ax = ax,alpha=0.6,cmap=cmap,add_colorbar=False,edgecolors = (0,0,0,0.1),linewidth = 0.2,vmax = vmax_value)

ax.scatter(wbb_loc['lon'],wbb_loc['lat'],c = 'white',marker='X',s=500,zorder = 10,label = "WBB",edgecolor = 'k',linewidth = 2)

cbformat = matplotlib.ticker.ScalarFormatter()   # create the formatter
cbformat.set_powerlimits((-2,2)) 
cbar = plt.colorbar(map,fraction=0.05,pad = 0.02,format=cbformat)
exp = gen_utils.calculate_exponent(max_value)
cbar.set_label(label = f"$CH_{4}$ emissions (Tonne$\\times 10^{{{exp}}}$)",size = labsize)
cbar.ax.tick_params(labelsize = labsize)

emissions_sum = float(plot_ds[spec].sum())
exp= gen_utils.calculate_exponent(emissions_sum)
sn_es = round(emissions_sum/10**(exp),2)
text = f"$CH_{4}$ sum ={sn_es}$\\times 10^{{{exp}}}$ Tonne"
ax.set_title("")

t1 = fig.text(0.275,0.92,text,fontsize = labsize)
t1.set_bbox(dict(facecolor = (1,1,1,0.8)))
fig.tight_layout()

# Explicitly set longitude and latitude ticks
xticks = [-112.1,-111.9,-111.7]
yticks = [40.4,40.6,40.8,41.0,41.2]
ax.set_xticks(xticks, crs=proj)
ax.set_yticks(yticks, crs=proj)
lon_formatter = LongitudeFormatter()
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_yticklabels([f"{y:.1f}°" for y in yticks], fontsize=labsize/2)
ax.set_xticklabels([f"{x:.1f}°" for x in xticks], fontsize=labsize/2)
ax.tick_params(
    left=False, labelleft=True,   # hide left labels
    top=False, labeltop=False,     # hide top labels
    right=False, labelright=False,   # show right labels
    bottom=False, labelbottom=True, # show bottom labels
    length=0  # remove the tick marks themselves
)
ax.set_xlabel('')  # remove x-axis label
ax.set_ylabel('')  # remove y-axis label


if savefig:
    fig_name = f'{fig_id}.png'
    fig.savefig(os.path.join(inventory_figures_path,fig_name),dpi=500,bbox_inches = "tight")
if showfig:
    plt.show()
else:
    plt.close()

#CO inventory plot
fig_id = f'2021_gra2pes_co_fullyear_inventory'
savefig = True
showfig = True

plot_ds = absolute_grapes_slv_yrsum_coarse

spec = 'CO'
cmap = 'Oranges'

fig = plt.figure(figsize = (15,15))
labsize = 44
proj = ccrs.PlateCarree()

ax = plt.axes(projection = proj)
ax.set_extent([map_extent['lon_min'],map_extent['lon_max'],map_extent['lat_min'],map_extent['lat_max']],crs=proj)
request = cimgt.GoogleTiles(style='satellite')
scale = 12.0 # prob have to adjust this
ax.add_image(request,int(scale))


max_value = plot_ds[spec].max().values
vmax_value = max_value*0.9
map = plot_ds[spec].plot.pcolormesh('lon','lat',ax = ax,alpha=0.6,cmap=cmap,add_colorbar=False,edgecolors = (0,0,0,0.1),linewidth = 0.2,vmax = vmax_value)

ax.scatter(wbb_loc['lon'],wbb_loc['lat'],c = 'white',marker='X',s=500,zorder = 10,label = "WBB",edgecolor = 'k',linewidth = 2)

cbformat = matplotlib.ticker.ScalarFormatter()   # create the formatter
cbformat.set_powerlimits((-2,2)) 
cbar = plt.colorbar(map,fraction=0.05,pad = 0.02,format=cbformat)
exp = gen_utils.calculate_exponent(max_value)
cbar.set_label(label = f"$CO$ emissions (Tonne$\\times 10^{{{exp}}}$)",size = labsize)
cbar.ax.tick_params(labelsize = labsize)

# Explicitly set longitude and latitude ticks
xticks = [-112.1,-111.9,-111.7]
yticks = [40.4,40.6,40.8,41.0,41.2]
ax.set_xticks(xticks, crs=proj)
ax.set_yticks(yticks, crs=proj)
lon_formatter = LongitudeFormatter()
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_yticklabels([f"{y:.1f}°" for y in yticks], fontsize=labsize/2)
ax.set_xticklabels([f"{x:.1f}°" for x in xticks], fontsize=labsize/2)
ax.tick_params(
    left=False, labelleft=True,   # hide left labels
    top=False, labeltop=False,     # hide top labels
    right=False, labelright=False,   # show right labels
    bottom=False, labelbottom=True, # show bottom labels
    length=0  # remove the tick marks themselves
)
ax.set_xlabel('')  # remove x-axis label
ax.set_ylabel('')  # remove y-axis label

emissions_sum = float(plot_ds[spec].sum())
exp=gen_utils.calculate_exponent(emissions_sum)
sn_es = round(emissions_sum/10**(exp),2)
text = f"$CO$ sum ={sn_es}$\\times 10^{{{exp}}}$ Tonne"
ax.set_title("")

t1 = fig.text(0.255,0.92,text,fontsize = labsize)
t1.set_bbox(dict(facecolor = (1,1,1,0.8)))
fig.tight_layout()

if savefig:
    fig_name = f'{fig_id}.png'
    fig.savefig(os.path.join(inventory_figures_path,fig_name),dpi=500,bbox_inches = "tight")
if showfig:
    plt.show()
else:
    plt.close()

EDGAR

#CO2 inventory plot
fig_id = f'2021_edgar_co2_fullyear_inventory'
savefig = True
showfig = True

spec = 'CO2'
plot_ds = edgar_dss[spec]['ds']
var = 'emissions'
cmap = 'Blues'

fig = plt.figure(figsize = (15,15))
labsize = 44
proj = ccrs.PlateCarree()

ax = plt.axes(projection = proj)
ax.set_extent([map_extent['lon_min'],map_extent['lon_max'],map_extent['lat_min'],map_extent['lat_max']],crs=proj)
request = cimgt.GoogleTiles(style='satellite')
scale = 12.0 # prob have to adjust this
ax.add_image(request,int(scale))


max_value = plot_ds[var].max().values
vmax_value = max_value*0.6
map = plot_ds[var].plot.pcolormesh('lon','lat',ax = ax,alpha=0.6,cmap=cmap,add_colorbar=False,edgecolors = (0,0,0,0.1),linewidth = 0.2,vmax = vmax_value)

ax.scatter(wbb_loc['lon'],wbb_loc['lat'],c = 'white',marker='X',s=500,zorder = 10,label = "WBB",edgecolor = 'k',linewidth = 2)

cbformat = matplotlib.ticker.ScalarFormatter()   # create the formatter
cbformat.set_powerlimits((-2,2)) 
cbar = plt.colorbar(map,fraction=0.05,pad = 0.02,format=cbformat)
exp = gen_utils.calculate_exponent(max_value)
cbar.set_label(label = f"$CO_{2}$ emissions (Tonne$\\times 10^{{{exp}}}$)",size = labsize)
cbar.ax.tick_params(labelsize = labsize)

#ax.legend(fontsize = labsize-8,loc = 4)

emissions_sum = float(plot_ds[var].sum())
exp= gen_utils.calculate_exponent(emissions_sum)
sn_es = round(emissions_sum/10**(exp),2)
text = f"$CO_{2}$ sum ={sn_es}$\\times 10^{{{exp}}}$ Tonne"

t1 = fig.text(0.24,0.92,text,fontsize = labsize)
t1.set_bbox(dict(facecolor = (1,1,1,0.8)))
fig.tight_layout()
ax.set_title("")

if savefig:
    fig_name = f'{fig_id}.png'
    fig.savefig(os.path.join(inventory_figures_path,fig_name),dpi=500,bbox_inches = "tight")
if showfig:
    plt.show()
else:
    plt.close()

#CH4 inventory plot
fig_id = f'2021_edgar_ch4_fullyear_inventory'
savefig = True
showfig = True

spec = 'CH4'
plot_ds = edgar_dss[spec]['ds']
var = 'emissions'
cmap = 'Reds'

fig = plt.figure(figsize = (15,15))
labsize = 44
proj = ccrs.PlateCarree()

ax = plt.axes(projection = proj)
ax.set_extent([map_extent['lon_min'],map_extent['lon_max'],map_extent['lat_min'],map_extent['lat_max']],crs=proj)
request = cimgt.GoogleTiles(style='satellite')
scale = 12.0 # prob have to adjust this
ax.add_image(request,int(scale))


max_value = plot_ds[var].max().values
vmax_value = max_value*0.9
map = plot_ds[var].plot.pcolormesh('lon','lat',ax = ax,alpha=0.6,cmap=cmap,add_colorbar=False,edgecolors = (0,0,0,0.1),linewidth = 0.2,vmax = vmax_value)

ax.scatter(wbb_loc['lon'],wbb_loc['lat'],c = 'white',marker='X',s=500,zorder = 10,label = "WBB",edgecolor = 'k',linewidth = 2)

cbformat = matplotlib.ticker.ScalarFormatter()   # create the formatter
cbformat.set_powerlimits((-2,2)) 
cbar = plt.colorbar(map,fraction=0.05,pad = 0.02,format=cbformat)
exp = gen_utils.calculate_exponent(max_value)
cbar.set_label(label = f"$CH_{4}$ emissions (Tonne$\\times 10^{{{exp}}}$)",size = labsize)
cbar.ax.tick_params(labelsize = labsize)

#ax.legend(fontsize = labsize-8,loc = 4)

emissions_sum = float(plot_ds[var].sum())
exp= gen_utils.calculate_exponent(emissions_sum)
sn_es = round(emissions_sum/10**(exp),2)
text = f"$CH_{4}$ sum ={sn_es}$\\times 10^{{{exp}}}$ Tonne"
ax.set_title("")

t1 = fig.text(0.275,0.92,text,fontsize = labsize)
t1.set_bbox(dict(facecolor = (1,1,1,0.8)))
fig.tight_layout()

if savefig:
    fig_name = f'{fig_id}.png'
    fig.savefig(os.path.join(inventory_figures_path,fig_name),dpi=500,bbox_inches = "tight")
if showfig:
    plt.show()
else:
    plt.close()

#CO inventory plot
fig_id = f'2018_edgar_co_fullyear_inventory'
savefig = True
showfig = True

spec = 'CO'
plot_ds = edgar_dss[spec]['ds']
var = 'emissions'
cmap = 'Oranges'

fig = plt.figure(figsize = (15,15))
labsize = 44
proj = ccrs.PlateCarree()

ax = plt.axes(projection = proj)
ax.set_extent([map_extent['lon_min'],map_extent['lon_max'],map_extent['lat_min'],map_extent['lat_max']],crs=proj)
request = cimgt.GoogleTiles(style='satellite')
scale = 12.0 # prob have to adjust this
ax.add_image(request,int(scale))


max_value = plot_ds[var].max().values
vmax_value = max_value*0.9
map = plot_ds[var].plot.pcolormesh('lon','lat',ax = ax,alpha=0.6,cmap=cmap,add_colorbar=False,edgecolors = (0,0,0,0.1),linewidth = 0.2,vmax = vmax_value)

ax.scatter(wbb_loc['lon'],wbb_loc['lat'],c = 'white',marker='X',s=500,zorder = 10,label = "WBB",edgecolor = 'k',linewidth = 2)

cbformat = matplotlib.ticker.ScalarFormatter()   # create the formatter
cbformat.set_powerlimits((-2,2)) 
cbar = plt.colorbar(map,fraction=0.05,pad = 0.02,format=cbformat)
exp = gen_utils.calculate_exponent(max_value)
cbar.set_label(label = f"$CO$ emissions (Tonne$\\times 10^{{{exp}}}$)",size = labsize)
cbar.ax.tick_params(labelsize = labsize)

#ax.legend(fontsize = labsize-8,loc = 4)

emissions_sum = float(plot_ds[var].sum())
exp=gen_utils.calculate_exponent(emissions_sum)
sn_es = round(emissions_sum/10**(exp),2)
text = f"$CO$ sum ={sn_es}$\\times 10^{{{exp}}}$ Tonne"
ax.set_title("")

t1 = fig.text(0.285,0.92,text,fontsize = labsize)
t1.set_bbox(dict(facecolor = (1,1,1,0.8)))
fig.tight_layout()

if savefig:
    fig_name = f'{fig_id}.png'
    fig.savefig(os.path.join(inventory_figures_path,fig_name),dpi=500,bbox_inches = "tight")
if showfig:
    plt.show()
else:
    plt.close()

EPA

#CH4 inventory plot
fig_id = f'2020_epa_ch4_fullyear_inventory'
savefig = True
showfig = True

spec = 'CH4'
plot_ds = absolute_epa_slv_yrsum
var = 'emissions'
cmap = 'Reds'

fig = plt.figure(figsize = (15,15))
labsize = 44
proj = ccrs.PlateCarree()

ax = plt.axes(projection = proj)
ax.set_extent([map_extent['lon_min'],map_extent['lon_max'],map_extent['lat_min'],map_extent['lat_max']],crs=proj)
request = cimgt.GoogleTiles(style='satellite')
scale = 12.0 # prob have to adjust this
ax.add_image(request,int(scale))


max_value = plot_ds[var].max().values
vmax_value = max_value*0.9
map = plot_ds[var].plot.pcolormesh('lon','lat',ax = ax,alpha=0.6,cmap=cmap,add_colorbar=False,edgecolors = (0,0,0,0.1),linewidth = 0.2,vmax = vmax_value)

ax.scatter(wbb_loc['lon'],wbb_loc['lat'],c = 'white',marker='X',s=500,zorder = 10,label = "WBB",edgecolor = 'k',linewidth = 2)

cbformat = matplotlib.ticker.ScalarFormatter()   # create the formatter
cbformat.set_powerlimits((-2,2)) 
cbar = plt.colorbar(map,fraction=0.05,pad = 0.02,format=cbformat)
exp = gen_utils.calculate_exponent(max_value)
cbar.set_label(label = f"$CH_{4}$ emissions (Tonne$\\times 10^{{{exp}}}$)",size = labsize)
cbar.ax.tick_params(labelsize = labsize)

#ax.legend(fontsize = labsize-8,loc = 4)

emissions_sum = float(plot_ds[var].sum())
exp= gen_utils.calculate_exponent(emissions_sum)
sn_es = round(emissions_sum/10**(exp),2)
text = f"$CH_{4}$ sum ={sn_es}$\\times 10^{{{exp}}}$ Tonne"
ax.set_title("")

t1 = fig.text(0.275,0.92,text,fontsize = labsize)
t1.set_bbox(dict(facecolor = (1,1,1,0.8)))
fig.tight_layout()

if savefig:
    fig_name = f'{fig_id}.png'
    fig.savefig(os.path.join(inventory_figures_path,fig_name),dpi=500,bbox_inches = "tight")
if showfig:
    plt.show()
else:
    plt.close()

Ratios

# Comparing ratios with EM27 means we only want to look at the daylight hours

daylight_hours = [14,15,16,17,18,19,20,21,22,23,0,1]
grapes_daylight_ds = absolute_grapes_slv_ds_coarse.where(absolute_grapes_slv_ds_coarse['datetime'].dt.hour.isin(daylight_hours),drop = True)
grapes_daylight_ds_total = grapes_daylight_ds.sel(sector='total')
grapes_daylight_ds = grapes_daylight_ds.sel(sector=grapes_daylight_ds.sector[grapes_daylight_ds.sector != "total"])

Yearly, daylight

# Aggregate to yearly emissions during daylight hours
absolute_grapes_slv_yrsum_daylight = grapes_daylight_ds_total.sum(dim='datetime',keep_attrs = True)
absolute_grapes_slv_yrsum_daylight = unit_converter.update_all_units(absolute_grapes_slv_yrsum_daylight, 'hr^-1','yr^-1')

# Print the ratios as reported in the manuscript
ir = xr_utils.InventoryRatios(absolute_grapes_slv_yrsum_daylight,molar_masses)
print(f"CO:CO2 ratio = {ir.get_total_ratio('CO','CO2')*1000} permil")
print(f"CH4:CO2 ratio = {ir.get_total_ratio('CH4','CO2')*1000} permil")

CO:CO2 ratio = 6.106649237470725 permil
CH4:CO2 ratio = 3.973079270326021 permil

Seasonal, daylight

# Define the seasons
seasons = {'DJF':{'months':[12,1,2]},
           'MAM':{'months':[3,4,5]},
           'JJA':{'months':[6,7,8]},
           'SON':{'months':[9,10,11]}
}

# Get and print the seasonal ratios

grapes_season_ratios = {}
for season in seasons.keys():
    months = seasons[season]['months']
    season_sum_ds = grapes_daylight_ds_total.where(grapes_daylight_ds_total['datetime'].dt.month.isin(months),drop = True).sum(dim='datetime',keep_attrs = True)
    season_sum_ds = unit_converter.update_all_units(season_sum_ds, 'hr^-1','season^-1')

    ir = xr_utils.InventoryRatios(season_sum_ds,molar_masses)
    ch4_co2 = ir.get_total_ratio('CH4','CO2')
    co_co2 = ir.get_total_ratio('CO','CO2')
    grapes_season_ratios[season] = {'ch4_co2':ch4_co2,'co_co2':co_co2}

    print(f"{season}")
    print(f"CH4:CO2 ratio = {ir.get_total_ratio('CH4','CO2')*1000} permil")
    print(f"CO:CO2 ratio = {ir.get_total_ratio('CO','CO2')*1000} permil")
    print("")

DJF
CH4:CO2 ratio = 3.3845618480081976 permil
CO:CO2 ratio = 5.506621140658151 permil

MAM
CH4:CO2 ratio = 3.9530358845863596 permil
CO:CO2 ratio = 6.105654891904161 permil

JJA
CH4:CO2 ratio = 4.570129302373411 permil
CO:CO2 ratio = 6.672402895752505 permil

SON
CH4:CO2 ratio = 4.0778679621754 permil
CO:CO2 ratio = 6.235658154424306 permil

Seasonal Sector

# Get the seasonal data separated by sector 

flat_data = []
for season, config in seasons.items():
    months = config['months']
    season_sum_ds = grapes_daylight_ds.where(
        grapes_daylight_ds['datetime'].dt.month.isin(months),
        drop=True
    ).sum(dim='datetime', keep_attrs=True)
    
    # Convert to per-season units
    season_sum_ds = unit_converter.update_all_units(season_sum_ds, 'hr^-1', 'season^-1')
    season_sum = season_sum_ds.sum(dim=['lat', 'lon'], keep_attrs=True)
    season_sum = unit_converter.update_all_units(season_sum, 'gridcell^-1', '')

    season_df = season_sum.to_dataframe().reset_index()

    for var in season_sum.data_vars:
        for sector in season_df['sector'].unique():
            value = season_df[season_df['sector'] == sector][var].values[0]
            flat_data.append([season, var, sector, value])

season_df = pd.DataFrame(flat_data, columns=['Season', 'Variable', 'Sector', 'Sum'])

# Group sectors
sectors_to_keep = ['OFFROAD', 'ONROAD_GAS', 'RES']
season_df['Sector_Grouped'] = season_df['Sector'].apply(lambda x: x if x in sectors_to_keep else 'Other')

# Summarize data
season_df_grouped = (
    season_df
    .groupby(['Season', 'Variable', 'Sector_Grouped'], as_index=False, observed=False)
    .agg({'Sum': 'sum'})
    .rename(columns={'Sector_Grouped': 'Sector'})
)

# Rename sectors for clarity
sector_rename_map = {
    'RES': 'Residential',
    'ONROAD_GAS': 'On Road Gas',
    'OFFROAD': 'Off Road Gas'
}
season_df_grouped['Sector'] = season_df_grouped['Sector'].replace(sector_rename_map)

palette = ['#ffa600','#003f5c','#7a5195','#ef5675']
sector_order = ['Residential', 'On Road Gas', 'Off Road Gas', 'Other']
var = 'CO2'
plot_utils.plot_stacked_seasonal_bars(
    df_grouped=season_df_grouped,
    var=var,
    fig_path=inventory_figures_path,
    fig_id=f'{var}_seasonal_inventory_sectors',
    palette=palette,
    sector_order=sector_order,
    savefig=True,
    showfig=True,
    legfontsize=40,
)

/uufs/chpc.utah.edu/common/home/u0890904/LAIR_1/atmos/utils/plot_utils.py:228: UserWarning: Tight layout not applied. The bottom and top margins cannot be made large enough to accommodate all Axes decorations.
  fig_leg.tight_layout()

Gridded Ratios

These cells allow for visualization of the grid cell level ratios.

# Define the inventory regridder
ir = xr_utils.InventoryRatios(absolute_grapes_slv_yrsum_daylight,molar_masses)

ratio_ds = ir.get_gc_ratio('CH4','CO2',num_quantile_filter = 0.25, denom_quantile_filter = 0.5)
ratio_ds['CO2_eq_sum'] = ratio_ds['CO2_eq_sum'].fillna(0)

ratio_id = 'ch4_co2'
fig_id = f'ch4co2_gridcell'
savefig = True
showfig = True

da = ratio_ds[ratio_id]*1000
alpha_var = 'CO2_eq_sum'
min_alpha = 0.3
alphas = ratio_ds[alpha_var].values
alphas = (alphas - alphas.min())/(alphas.max()-alphas.min())*(1-min_alpha)+min_alpha

labsize = 40
proj = ccrs.PlateCarree()

fig = plt.figure(figsize = (12,12))
ax = plt.axes(projection = proj)
ax.set_extent([map_extent['lon_min'],map_extent['lon_max'],map_extent['lat_min'],map_extent['lat_max']],crs=proj)
request = cimgt.GoogleTiles(style='satellite')
scale = 12.0 # prob have to adjust this
ax.add_image(request,int(scale))

map = da.plot.pcolormesh('lon','lat',ax = ax,alpha=alphas,cmap='viridis',add_colorbar=False,edgecolors = (0,0,0,1),linewidth = 2,vmin = 0,vmax=10)#,vmin=0,vmax = 10)
da.plot.pcolormesh('lon','lat', ax = ax,facecolors='none', edgecolors=(0, 0, 0, 1), linewidth=2,add_colorbar=False)

ax.scatter(wbb_loc['lon'],wbb_loc['lat'],c = 'white',marker='X',s=500,zorder = 10,label = "WBB",edgecolor = 'k',linewidth = 2)

for lat in da.lat.values:
    for lon in da.lon.values:
        value = da.sel(lat=lat,lon=lon).values
        if value != value:
            continue
        #plt.text(lon,lat,f'{value:.1f}',color='k',ha='center',va='center',size=15)

source_types_to_plot = ['landfills','refineries']
for source_type,sources in point_sources.items():
    if source_type not in source_types_to_plot:
        continue
    locs = sources['locs']
    plot_details = sources['plot_details']
    for loc in locs:
        ax.scatter(loc[0],loc[1],**plot_details)
    

ax.get_figure().gca().set_title("")

cbar = plt.colorbar(map,fraction=0.045,pad = 0.02,location = 'right')
#cbar.set_label(label = r"${\alpha_{CH4}}^o/_{oo}$",size = labsize)
cbar.ax.tick_params(labelsize = labsize)
cbar.set_ticks([0,2,4,6,8,10])
cbar.ax.set_yticklabels([0,2,4,6,8,">10"])

# total_ratio = ir.get_total_ratio('CH4','CO2')*1000
# text = r"${\alpha}_{{CH_{4}},{EI},{tot}}$"
# text += f" = {total_ratio:.2f}"
# text = text+ "$^o/_{oo}$"
# t1 = fig.text(0.405,0.832,text,fontsize = labsize-5)
# t1.set_bbox(dict(facecolor = 'white'))

# Explicitly set longitude and latitude ticks
xticks = [-112.1,-111.9,-111.7]
yticks = [40.4,40.6,40.8,41.0,41.2]
ax.set_xticks(xticks, crs=proj)
ax.set_yticks(yticks, crs=proj)
lon_formatter = LongitudeFormatter()
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_yticklabels([f"{y:.1f}°" for y in yticks], fontsize=labsize/2)
ax.set_xticklabels([f"{x:.1f}°" for x in xticks], fontsize=labsize/2)
ax.tick_params(
    left=False, labelleft=True,   # hide left labels
    top=False, labeltop=False,     # hide top labels
    right=False, labelright=False,   # show right labels
    bottom=False, labelbottom=True, # show bottom labels
    length=0  # remove the tick marks themselves
)
ax.set_xlabel('')  # remove x-axis label
ax.set_ylabel('')  # remove y-axis label

if savefig:
    fig_name = f'{fig_id}.png'
    fig.savefig(os.path.join(inventory_figures_path,fig_name),dpi=500,bbox_inches = "tight")
if showfig:
    plt.show()
else:
    plt.close()

ratio_ds = ir.get_gc_ratio('CO','CO2',num_quantile_filter = 0.25, denom_quantile_filter = 0.5)
ratio_ds['CO2_eq_sum'] = ratio_ds['CO2_eq_sum'].fillna(0)

ratio_id = 'co_co2'
fig_id = f'coco2_gridcell'
savefig = True
showfig = True

da = ratio_ds[ratio_id]*1000
alpha_var = 'CO2_eq_sum'
min_alpha = 0.3
alphas = ratio_ds[alpha_var].values
alphas = (alphas - alphas.min())/(alphas.max()-alphas.min())*(1-min_alpha)+min_alpha

labsize = 40
proj = ccrs.PlateCarree()

fig = plt.figure(figsize = (12,12))
ax = plt.axes(projection = proj)
ax.set_extent([map_extent['lon_min'],map_extent['lon_max'],map_extent['lat_min'],map_extent['lat_max']],crs=proj)
request = cimgt.GoogleTiles(style='satellite')
scale = 12.0 # prob have to adjust this
ax.add_image(request,int(scale))

map = da.plot.pcolormesh('lon','lat',ax = ax,alpha=alphas,cmap='viridis',add_colorbar=False,edgecolors = (0,0,0,1),linewidth = 2,vmin = 0,vmax=10)#,vmin=0,vmax = 10)
da.plot.pcolormesh('lon','lat', ax = ax,facecolors='none', edgecolors=(0, 0, 0, 1), linewidth=2,add_colorbar=False)

ax.scatter(wbb_loc['lon'],wbb_loc['lat'],c = 'white',marker='X',s=500,zorder = 10,label = "WBB",edgecolor = 'k',linewidth = 2)

for lat in da.lat.values:
    for lon in da.lon.values:
        value = da.sel(lat=lat,lon=lon).values
        if value != value:
            continue
        #plt.text(lon,lat,f'{value:.1f}',color='k',ha='center',va='center',size=15)

# source_types_to_plot = ['landfills','refineries']
# for source_type,sources in point_sources.items():
#     if source_type not in source_types_to_plot:
#         continue
#     locs = sources['locs']
#     plot_details = sources['plot_details']
#     for loc in locs:
#         ax.scatter(loc[0],loc[1],**plot_details)
    

ax.get_figure().gca().set_title("")

cbar = plt.colorbar(map,fraction=0.045,pad = 0.02,location = 'right')
#cbar.set_label(label = r"${\alpha_{CH4}}^o/_{oo}$",size = labsize)
cbar.ax.tick_params(labelsize = labsize)
cbar.set_ticks([0,2,4,6,8,10])
cbar.ax.set_yticklabels([0,2,4,6,8,">10"])

# total_ratio = ir.get_total_ratio('CH4','CO2')*1000
# text = r"${\alpha}_{{CH_{4}},{EI},{tot}}$"
# text += f" = {total_ratio:.2f}"
# text = text+ "$^o/_{oo}$"
# t1 = fig.text(0.405,0.832,text,fontsize = labsize-5)
# t1.set_bbox(dict(facecolor = 'white'))

# Explicitly set longitude and latitude ticks
xticks = [-112.1,-111.9,-111.7]
yticks = [40.4,40.6,40.8,41.0,41.2]
ax.set_xticks(xticks, crs=proj)
ax.set_yticks(yticks, crs=proj)
lon_formatter = LongitudeFormatter()
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_yticklabels([f"{y:.1f}°" for y in yticks], fontsize=labsize/2)
ax.set_xticklabels([f"{x:.1f}°" for x in xticks], fontsize=labsize/2)
ax.tick_params(
    left=False, labelleft=True,   # hide left labels
    top=False, labeltop=False,     # hide top labels
    right=False, labelright=False,   # show right labels
    bottom=False, labelbottom=True, # show bottom labels
    length=0  # remove the tick marks themselves
)
ax.set_xlabel('')  # remove x-axis label
ax.set_ylabel('')  # remove y-axis label

if savefig:
    fig_name = f'{fig_id}.png'
    fig.savefig(os.path.join(inventory_figures_path,fig_name),dpi=500,bbox_inches = "tight")
if showfig:
    plt.show()
else:
    plt.close()