<jupyter_start><jupyter_text>Some usefull functions and housekeepingInitalize a `sqlite` database<jupyter_code>sqdb = ut.sqlitedb('fcat_biomass')<jupyter_output><empty_output><jupyter_text>Functions for<jupyter_code>def distFromRange(total, maxr = 32, minr = 2):
    av = (maxr + minr)/2
    stdev = (float(maxr) - float(minr))/4 
    d_frac = (total-np.floor(total))*np.random.normal(av, stdev, 1).clip(min=0)[0]
    t_bdt = np.random.normal(av,stdev,np.floor(total)).clip(min=0)
    return np.append(t_bdt, d_frac)
def sumFromDist(total, maxr = 0.32, minr = 0.02):
    av = (maxr + minr)/2
    stdev = (float(maxr) - float(minr))/4 
    d_frac = (total-np.floor(total))*np.random.normal(av, stdev, 1).clip(min=0)[0]
    t_bdt = np.sum(np.random.normal(av,stdev,np.floor(total)).clip(min=0))
    return (d_frac+t_bdt)<jupyter_output><empty_output><jupyter_text>Black Carbon from pile burningThe  CARB criteria pollutant emissions inventory reports Particualte Matter (PM 2.5) emissions from anthropogenic buring of forest residuals. The following estimates elemental carbon (Black Carbon) based on empirically derived relationships between PM2.5 and EC from.1. Ward DE, Hardy CC. Organic and elemental profiles for smoke from prescribed fires. In: <NAME>, editor. International specialty conference of the Air and Waste Management Association [Internet]. San Francisco: Air and Waste Management Association; 1989. Available from: http://www.frames.gov/documents/smoke/serdp/ward_hardy_1989a.pdf<jupyter_code>ward = ut.gData('13UQtRfNBSJ81PXxbYSnB2LrjHePNcvhJhrsxRBjHpoY', 475419971)
#Units are ratio of EC to PM produced 
wardDF = ward[['source','pct_sm','pct_f','tc_f_est','tc_f_cv','tc_s_est','tc_s_cv','ec_f_est','ec_f_cv', 'ec_s_est','ec_s_cv','ecton-1_h_pm','ecton-1_l_pm']].transpose()
wardDF.columns = wardDF.iloc[0]
eFact = wardDF.to_dict()
d = dict(zip([eFact[i]['source'] for i in eFact.keys()],['ALL VEGETATION','WILDLAND FIRE USE (WFU)','FOREST MANAGEMENT']))
for k in eFact.keys():
    eFact[k]['arblink'] = d[k]
pd.DataFrame.from_dict(eFact).transpose().to_sql('bc_pm_ratio', sqdb['cx'], if_exists= 'replace')
#print tabulate(pd.DataFrame.from_dict(eFact).transpose(), headers = ['CARB CPE Cat.','BC/t PM 2.5 (high)','BC/t PM 2.5 (low)', 'Source'], tablefmt="pipe")
pd.DataFrame.from_dict(eFact).transpose()
pd.DataFrame.from_dict(eFact)
ward
ec_rat_table = pd.DataFrame([ward.source, ward.tc_f_est*ward.ec_f_est, ward.tc_f_cv, ward.ec_f_cv,ward.tc_s_est*ward.ec_s_est, ward.tc_s_cv, ward.ec_s_cv]).transpose()
ec_rat_table.to_sql('ec_ratios', sqdb['cx'], if_exists = 'replace')
ec_rat_table
print tabulate(ec_rat_table, headers = ['Source',
                                        'BC/t PM 2.5 (F, est.)',
                                        'TC/t PM 2.5 (F, CV)',
                                        'BC/t TC (F, CV)',
                                        'BC/t PM 2.5 (S, est.)',
                                        'TC/t PM 2.5 (S, CV)',
                                        'BC/t TC (S, CV)'], tablefmt="pipe")
ecPct = ward[['ec_f_est', 'pct_f','pct_sm','tc_f_est','tc_f_cv','tc_s_est','tc_s_cv','ec_f_cv','ec_s_est', 'ec_s_cv']].set_index(ward.source).to_dict('index')<jupyter_output><empty_output><jupyter_text>Black Carbon Global Warming Potential[Several estimates](https://docs.google.com/spreadsheets/d/13UQtRfNBSJ81PXxbYSnB2LrjHePNcvhJhrsxRBjHpoY/editgid=195715938) exist for the GWP of Black Carbon. We use GWP20 estimates for black carbon from the CARB Short-Lived Climate Pollutant Strategy and from Fuglestvedt et. al 2000. References1. California Air Resources Board. [Short Lived Climate Pollutant Reduction Strategy](https://docs.google.com/uc?id=0B9-9Vlx0SkkFWmxxQ2xtSkNxSlU&export=download). Sacramento, CA; 2015. 2. Fuglestvedt JS, Berntsen TK, Godal O, Skodvin T. Climate implications of GWP-based reductions in greenhouse gas emissions. Geophys Res Lett [Internet]. 2000 Feb 1 [cited 2015 Sep 5];27(3):409–12. Available from: http://doi.wiley.com/10.1029/1999GL010939<jupyter_code>bc_gwp = ut.gData('<KEY>', 195715938)
bc_gwp.to_sql('bc_gwp', sqdb['cx'], if_exists = 'replace')
#print tabulate(bc_gwp.drop('est_id', 1), headers = [i for i in bc_gwp.drop('est_id', 1).columns],tablefmt="pipe")
bc_gwp = bc_gwp.set_index(bc_gwp.est_id).drop('est_id',1)
bc_gwp<jupyter_output><empty_output><jupyter_text>CARB Criteria Air Pollutant (CAP) emissions data 2015Particulate Matter, Carbon Monoxide, Nitrogen Oxide, Sulfur Dioxide, Lead, and Reactive Organic Gases<jupyter_code>cpe_data = pd.read_csv('http://www.arb.ca.gov/app/emsinv/2013/emsbyeic.csv?F_YR=2015&F_DIV=0&F_SEASON=A&SP=2013&SPN=2013_Almanac&F_AREA=CA')
cpe_data.columns = [i.lower() for i in cpe_data.columns]
cpe_data.to_sql('cpe_2015', sqdb['cx'], if_exists = 'replace')
cpe_2015 = pd.read_sql('''SELECT
                                eicsoun,
                                pm.source type,
                                pm2_5 pm25_tpd,
                                pm2_5*365*"ecton-1_h_pm" as t_ec_high,
                                pm2_5*365*"ecton-1_l_pm" as t_ec_low,
                                pm2_5*365*"ecton-1_h_pm"*gwp_20 as co2e_high,
                                pm2_5*365*"ecton-1_l_pm"*gwp_20 as co2e_low,
                                pm2_5*365*"ecton-1_l_pm"*(gwp_20-gwp_20_std) as co2e20_low_m1std,
                                pm2_5*365*"ecton-1_l_pm"*(gwp_20+gwp_20_std) as co2e20_low_p1std,
                                pm2_5*365*"ecton-1_h_pm"*(gwp_20-gwp_20_std) as co2e20_hi_m1std,
                                pm2_5*365*"ecton-1_h_pm"*(gwp_20+gwp_20_std) as co2e20_hi_p1std,
                                pm2_5*365*"ecton-1_l_pm"*(gwp_100-gwp_100_std) as co2e100_low_m1std,
                                pm2_5*365*"ecton-1_l_pm"*(gwp_100+gwp_100_std) as co2e100_low_p1std,
                                pm2_5*365*"ecton-1_h_pm"*(gwp_100-gwp_100_std) as co2e100_hi_m1std,
                                pm2_5*365*"ecton-1_h_pm"*(gwp_100+gwp_100_std) as co2e100_hi_p1std,
                                gwp.source
                            FROM cpe_2015
                            JOIN bc_pm_ratio pm  on  (arblink = eicsoun)
                            CROSS JOIN bc_gwp gwp 
                            WHERE eicsoun in ('FOREST MANAGEMENT','WILDLAND FIRE USE (WFU)','ALL VEGETATION')''',
            con = sqdb['cx'])

foo = pd.melt(cpe_2015.drop('eicsoun',1), id_vars = ['type', 'source'])
bar = foo[foo['variable'].str.contains("co2e")].dropna()
bar['type'] = bar.type.astype('category')
bar.value = bar.value
bar<jupyter_output><empty_output><jupyter_text>Make the BC GHG emissions plot A secondary apporoach based on Ward (1989) and Jenkins (1996) uses a range of 2-32% of PM as BC. 1. <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, et al. Atmospheric Pollutant Emissions Factors From Open Burning of Agricultural and Forest Biomass By Wind Tunnel Simulations [Internet]. Sacramento, CA; 1996. Available from: http://www.arb.ca.gov/ei/speciate/r01t20/rf9doc/refnum9.htm<jupyter_code>pm2015=pd.read_sql('''select eicsoun,
                    pm2_5*365 as pm25
                from cpe_2015
                WHERE eicsoun in ('FOREST MANAGEMENT','WILDLAND FIRE USE (WFU)','ALL VEGETATION');''', con = sqdb['cx'])
pm2015=pm2015.set_index([['Wildfire','Pile Burn','Prescribed']])
pm2015
def ecDist(ecEst, ecCV, PM):
    '''
    PM is a mass measure, its intended to be PM2.5
    ecEst is the estimate of the percentage of elemental carbon comprising the PM
    ecCV is the coefficient of variation around the estimate of elemental carbon
    -----
    returns a random selection from a normal distribution of size `len(pm)`
    centered on `ecEst` with standard deviation of `ecCV`*`ecEst`
    '''
    ecStdev = ecCV * ecEst
#    ec_frac = (PM-np.floor(PM))*np.random.normal(ecEst, ecStdev, 1).clip(min=0)[0]
    t_ec = np.random.normal(ecEst,ecStdev,PM).clip(min=0)#+ec_frac
    return t_ec #+ ec_frac
ecPct
for k in ecPct.keys():
    pm = pm2015.loc[k]['pm25']
    ecPct[k]['t_pm'] = pm
    #PM2.5 smoldering 
    pm_sm = ecPct[k]['pct_sm']*pm
    ecPct[k]['pm_sm'] = pm_sm
    #TC smoldering
    tc_sm = ecPct[k]['tc_s_est']*pm_sm
    ecPct[k]['tc_sm'] = tc_sm
    #PM2.5 flaming
    pm_f = ecPct[k]['pct_f']*pm
    ecPct[k]['pm_f'] = pm_f
    #TC flaming
    tc_f = ecPct[k]['tc_f_est']*pm_f
    ecPct[k]['tc_f'] = tc_f
    res1k = []
    for t in range(1000):
        rnd = t
        # Total Carbon in PM
        tc_s = ecDist(ecPct[k]['tc_s_est'], ecPct[k]['tc_s_cv'], pm)
        tc_f = ecDist(ecPct[k]['tc_f_est'], ecPct[k]['tc_f_cv'], pm)
        #Elemental Carbon in Total Carbon
        ec_s = ecDist(ecPct[k]['ec_s_est'], ecPct[k]['ec_s_cv'], pm)
        ec_f = ecDist(ecPct[k]['ec_f_est'], ecPct[k]['ec_f_cv'], pm)
        ec1tpm = (tc_s*ec_s) + (tc_f*ec_f)
        ec_total_rnd = sum(ec1tpm)
        ec_total_gwp = ec_total_rnd * bc_gwp.loc['carb_slcp']['gwp_20']
        res1k.append([k,k+str(rnd),ec_total_rnd, ec_total_gwp])
    ecPct[k]['ec_tdist'] = pd.DataFrame(res1k)
    #ecPct[k]['ecSMDist'] = np.array([(ecDist(ecPct[k]['ec_s_est'], ecPct[k]['ec_s_cv'], pm)) for i in range(1000)])
    #ecPct[k]['ecFDist'] = np.array([(ecDist(ecPct[k]['ec_f_est'], ecPct[k]['ec_f_cv'], pm)) for i in range(1000)])
    #ecPct[k]['tECDist'] = ecPct[k]['ecSMDist']+ecPct[k]['ecFDist']
    #ecPct[k]['df'] = pd.DataFrame(np.column_stack((ecPct[k]['tECDist'], [k]*len(ecPct[k]['tECDist']))))
ecPct['Prescribed']['tc_f']
(pm_sm*0.2*3200)+(pm_f*0.2*3200)
df=pd.concat([ecPct[i]['ec_tdist'] for i in ecPct.keys()])
df.columns = ['source', 'sourcernd', 'mt_ec','co2e_ec']
t=sns.FacetGrid(df,col='source', sharey=False)
t.map(sns.boxplot, 'source','co2e_ec')
t.set_ylabels('t $CO_2e$')
[c.xaxis.set_visible(False) for c in t.axes[0]]
sns.plt.subplots_adjust(top = 0.8)
t.fig.suptitle('Black carbon emissions in CO2 equivalent units from burning in CA, 2015', fontsize= 15)
t.fig.text(0.1, -0.08,'''Sources: CARB Criteria Pollutant Emissions Inventory(2015), Ward and Hardy (1989)''',
      fontsize=10)
for f in ['.png','.pdf']:
    t.savefig('graphics/bc_prob_gwp{0}'.format(f))
122*0.66
pd.read_csv('fera_pile_cemissions.csv', header=1)<jupyter_output><empty_output><jupyter_text>Forested Lands and Wood Products Biodegradable Carbon Emissions & Sinks (MMTCO2)ARB (2007). Technical support document for Land Use, Land Use Change & Forestry - Biodegradable Carbon Emissions & Sinks. Table available as a pdf document at: http://www.arb.ca.gov/cc/inventory/archive/tables/net_co2_flux_2007-11-19.pdf<jupyter_code>t=foo[foo['variable'].str.contains("co2e20")].dropna()
co2 = ut.gData('1GDdquzrCoq2cxVN2fbCpP4gwi2yrMnONNrWbfhZKZu4', 1636249481)
co2.columns = co2.iloc[0]
co2plot=pd.melt(co2.reindex(co2.index.drop(0)), id_vars = ['Year'])
co2plot.columns = ['sc_cat','year','mmtco2e']
co2plot.to_sql('arb_co2', sqdb['cx'], if_exists = 'replace')
cdata = pd.read_sql('''select sc_cat, avg(mmtco2e) from arb_co2 where sc_cat in ('Forest and rangeland fires', 'Timber harvest slash')  group by sc_cat''', con = sqdb['cx'])
co2
cdata
#Output GitHub Markdown using the following:
#print tabulate(cdata, headers = ['Source Category','MMTCO2'],tablefmt="pipe")<jupyter_output><empty_output><jupyter_text>GHG equivalent emissions from management residualsTo arrive at an estimate of total annual emissions from buring forest management residuals in CO2 equivalent terms from published CARB estimates we can combine the CO2 emissions reported for 2004 in the LULUC Biodegradable Carbon Emissions and Sinks with black carbon emissions extrapolated from the CARB Criteria Air Pollutant Emissions inventory estimates. The dime discreppancy between the 2004 and 2015 is acknowledged as an irreconcilable source of uncertainty in this estimation among others. This does however reflect that a baseline of substantial emissions from forest management residuals has been reported in CARB emissions inventories and should be recognized as a baseline condition.<jupyter_code>t=foo[foo['variable'].str.contains("co2e20")].dropna()
tE = pd.DataFrame([cdata['avg(mmtco2e)'][1]*1.10231,
              avg(t[t['type'].str.contains('piles')]['value'])/1000000,
             (cdata['avg(mmtco2e)'][1]*1.10231)+(avg(t[t['type'].str.contains('piles')]['value'])/1000000)],columns = ['Mt CO2e'])
tE['Source']=['CO2 pile buring', 'CO2e BC pile burning', 'Total Mt CO2e']
tE.to_sql('pile_em', sqdb['cx'], if_exists = 'replace')
#print tabulate(tE, headers = tE.columns.tolist(), tablefmt ='pipe')
tE<jupyter_output><empty_output><jupyter_text>Timber Products OutputThe TPO estimates logging redisues produced from commercial timber harvesting operations.<jupyter_code>tpoData = ut.gData('1GDdquzrCoq2cxVN2fbCpP4gwi2yrMnONNrWbfhZKZu4', 872275354, hrow=1)
tpoData<jupyter_output><empty_output><jupyter_text>Biomass residuals from non-commercial management activitiesData from TPO does not account for forest management activities that do not result in commercial products (timber sales, biomass sales). To estimate the amount of residual material produced from non commercial management activities we use data from the US Forest Service (FACTS) and from CalFires timber harvest plan data.  Forest Service ACtivity Tracking System (FACTS)Data from TPO does not account for forest management activities that do not result in commercial products (timber sales, biomass sales). We use a range of 10-35 BDT/acre to convert acres reported in FACTS to volume.<jupyter_code>pd.read_excel('FACTS_Tabular_092115.xlsx', sheetname = 'CategoryCrosswalk').to_sql('facts_cat', sqdb['cx'], if_exists = 'replace')
pd.read_csv('pd/facts_notimber.csv').to_sql('facts_notimber', sqdb['cx'], if_exists='replace')<jupyter_output><empty_output><jupyter_text>Querying FACTSThe USFS [reports](http://data.fs.usda.gov/geodata/edw/datasets.php?dsetParent=Activities) Hazardous Fuels Treatment (HFT) activities as well as Timber Sales (TS) derived from the FACTS database. We use these two datasets to estimate the number of acres treated that did not produce commercial material (sawlogs or biomass) and where burning was not used. The first step is to elimina all treatments in the HFT dataset that included timber sales. We accomplish this by eliminating all rows in the HFT dataset that have identical `FACTS_ID` fields in the TS dataset. We further filter the HFT dataset by removing any planned but not executed treatements (`nbr_units1 >0` below -- `nbr_units1` references `NBR_UNITS_ACCOMPLISHED` in the USFS dataset, see metadata for HFT [here](http://data.fs.usda.gov/geodata/edw/edw_resources/meta/S_USA.Activity_HazFuelTrt_PL.xml)), and use text matching in the 'ACTIVITY' and 'METHOD' fields to remove any rows that contain reference to 'burning' or 'fire'. Finally, we remove all rows that that reference 'Biomass' in the method category as it is assumed that this means material was removed for bioenergy.<jupyter_code>usfs_acres = pd.read_sql('''select
                                sum(nbr_units1) acres,
                                method,
                                strftime('%Y',date_compl) year,
                                cat."ACTIVITY" activity,
                                cat."TENTATIVE_CATEGORY" r5_cat
                            from facts_notimber n 
                                join facts_cat cat
                                    on (n.activity = cat."ACTIVITY") 
                            where date_compl is not null
                                and nbr_units1 > 0
                                and cat."TENTATIVE_CATEGORY" != 'Burning'
                                and cat."ACTIVITY" not like '%ire%'
                                and method not like '%Burn%'
                                and method != 'Biomass'
                            group by cat."ACTIVITY",
                                year,
                                method,
                                cat."TENTATIVE_CATEGORY"
                            order by year;''', con = sqdb['cx'])<jupyter_output><empty_output><jupyter_text>Converting acres to cubic feetFACTS reports in acres. To estimate the production of biomass from acres treated we use a range of 10-35 BDT/acre. We assume that actual biomass residuals per acre are normally distributed with a mean of 22.5 and a standard deviation of (35-10)/4 = 6.25<jupyter_code>def sumBDT(ac, maxbdt = 35, minbdt = 10):
    av = (maxbdt + minbdt)/2
    stdev = (float(maxbdt) - float(minbdt))/4 
    d_frac = (ac-np.floor(ac))*np.random.normal(av, stdev, 1).clip(min=0)[0]
    t_bdt = np.sum(np.random.normal(av,stdev,np.floor(ac)).clip(min=0))
    return d_frac+t_bdt
usfs_acres['bdt'] = usfs_acres['acres'].apply(sumBDT)
usfs_an_bdt = usfs_acres.groupby(['year']).sum()<jupyter_output><empty_output><jupyter_text>Weighted average wood densityAverage wood density weighted by harvested species percent. Derived from McIver and Morgan,  [Table 4](https://docs.google.com/spreadsheets/d/138FWlGeW57MKdcz2UkWxtWV4o50SZO8sduB1R6JOFp8/editgid=393414465)<jupyter_code>wood_dens = ut.gData('138FWlGeW57MKdcz2UkWxtWV4o50SZO8sduB1R6JOFp8', 1297253755)
wavg_dens =sum(wood_dens.pct/100 * wood_dens.density_lbscuft)<jupyter_output><empty_output><jupyter_text>Annual unutilized management residuals- [x] Public lands non-commercial management residuals - [ ] Private land non-commercial management residuals- [x] Public lands logging residuals- [x] Private lands logging residuals<jupyter_code>cat_codes = {'nf_ncmr': 'Unburned, non-commercial management residuals from National Forest lands',
             'nf_lr': 'Logging residuals generated from timber sales on National Forest lands',
             'opriv_lr': 'Logging residuals generated from timber sales on non-industrial private forest lands',
             'fi_lr': 'Logging residuals generated from timber sales on industrial private lands',
             'opub_lr': 'Logging residuals generated from timber sales on industrial private lands'}
usfs_an_bdt['cuft']= usfs_an_bdt.bdt *wavg_dens
resid_stats=pd.DataFrame((usfs_an_bdt.iloc[6:,2]/1000000).describe())
resid_stats.columns = ['nf_ncmr']
resid_stats['nf_lr']=tpoData[tpoData.ownership.str.contains('National Forest')]['loggingresidues'].describe()
resid_stats['opriv_lr']=tpoData[tpoData.ownership.str.contains('Other Private')]['loggingresidues'].describe()
resid_stats['fi_lr']=tpoData[tpoData.ownership.str.contains('Forest Industry')]['loggingresidues'].describe()
resid_stats['opub_lr']=tpoData[tpoData.ownership.str.contains('Other Public')]['loggingresidues'].describe()
resid_stats
print tabulate(resid_stats, headers = resid_stats.columns.tolist(), tablefmt ='pipe')
import os
[os.path.splitext(i)[0] for i in os.listdir('lf/') if os.path.splitext(i)[1] =='.csv']
ureg = UnitRegistry()
ureg.define('cubic foot =  cubic_centimeter/ 3.53147e-5 = cubic_foot' )
ureg.define('million cubic foot =  cubic_foot*1000000 = MMCF' )
ureg.define('board foot sawlog = cubic_foot / 5.44 = BF_saw')
ureg.define('board foot veneer = cubic_foot / 5.0 = BF_vo')
ureg.define('board foot bioenergy = cubic_foot / 1.0 = BF_bio')
ureg.define('bone-dry unit = cubic_foot * 96 = BDU')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1292655014565127, 'fraction_numerical': 0.036933000416146486, 'mean_word_length': 3.7692880178615726, 'pattern_counts': {'":': 0, '<': 65, '<?xml version=': 0, '>': 67, 'https://': 3, 'lorem ipsum': 0, 'www.': 4, 'xml': 1}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '484146', 'license': 'MIT', 'n_tokens_mistral': 7375, 'n_tokens_neox': 6744, 'n_words': 1517, 'path': 'biomass_burning.ipynb', 'repo_name': 'WoodResourcesGroup/RoundwoodHarvestGHG', 'size': '19293'}