HWu = pd.read_sql('''SELECT *
                        FROM so4
                        WHERE harvestslash = "X"
                        AND processingresidues = "X" 
                        AND "post-usewoodproduct" = "X"
                        AND stumps is null''', sqdb['cx'],index_col = 'index')
HWo = pd.read_sql('''SELECT *
                        FROM so4
                        WHERE harvestslash is null
                        AND processingresidues = "X" 
                        AND stumps is null
                        AND "post-usewoodproduct" = "X"''', sqdb['cx'], index_col = 'index')


#HWo
print tabulate(HWo[['reference','df']], headers = ['index','reference','displacement factor'],tablefmt="pipe")

#HWu
print tabulate(HWu[['reference','df']], headers = ['index','reference','displacement factor'],tablefmt="pipe")

constants = {'me' : {'value':0.5,
                     'desc': 'Mill Efficiency'},
             'DFu' : {'value': np.average(HWu.df),
                      'desc': 'Displacement factor with logging residual utilization',
                      'source': '''\cite{Sathre2010}'''},
             'DFo' : {'value': np.average(HWo.df),
                      'desc': 'Displacement factor without logging residual utilization'},
             'wDens' : {'value': sum(wood_dens.pct/100 * wood_dens.density_lbscuft),
                        'units' : 'lbs/cuft',
                        'desc': 'average harvested wood density weighted by species harvested',
                        'source': '\cite{Mciver2012}'}
             }

constants['wDens']['value']

"""
Explanation: Wood products DF with and without logging slash utilization
Using studies from Sathre and O'Connor in the US.
End of explanation
"""


tpoData = ut.gData('1GDdquzrCoq2cxVN2fbCpP4gwi2yrMnONNrWbfhZKZu4', 872275354, hrow=1)
tpoData.to_sql('tpo', sqdb['cx'], if_exists = 'replace')
print tabulate(tpoData, headers = ['Ownership','Roundwood Products','Logging Residues', 'Year'],tablefmt="pipe")

"""
Explanation: Timber Products Output
The TPO estimates logging redisues produced from commercial timber harvesting operations. The follwoing is in million cubic feet (MCF)
End of explanation
"""


pd.read_csv('boe_hh.csv').to_sql('boe',sqdb['cx'], if_exists = 'replace')

pd.read_sql('select * from boe', sqdb['cx'])

"""
Explanation: Board of Equalization Data
The board of equalization
End of explanation
"""


mm_histHarvest = ut.gData('13UQtRfNBSJ81PXxbYSnB2LrjHePNcvhJhrsxRBjHpoY', 2081880100).fillna(value=0)

mm_histHarvest.to_sql('mm_hist', sqdb['cx'], if_exists = 'replace')

mm_histHarvest

"""
Explanation: McIver and Morgan annual in cubic
Figure 2 from morgan and mciver presents total roundwood harvest from 1947 through 2012 in MMBF. To convert MMBF to MCF we use a sawlog conversion of 5.44. This is an approximation as the actual sawlog conversion varies with the log size on average over time has changed.
End of explanation
"""


bioEnergy = ut.gData('138FWlGeW57MKdcz2UkWxtWV4o50SZO8sduB1R6JOFp8', 529610043)
bioEnergy.set_index('producttype').transpose().to_sql('mciver_bio', sqdb['cx'], if_exists = 'replace')
bio_pct = pd.read_sql('select "index" as year,"Bioenergy"/100 as biopct from mciver_bio where "Bioenergy" is not null', sqdb['cx'])
bio_dict = bio_pct.set_index('year').to_dict('index')
print tabulate(bio_pct, headers = ['year', 'bioenergy % of harvest'],tablefmt="pipe")

def bioPct(year):
#    if year < 1980:
#        return 0
    if year in bio_dict.keys():
        return bio_dict[year]['biopct']
    else:
        return np.average(bio_pct.biopct)

"""
Explanation: Bioenergy consumption
To apply the apropriate DF for harvested wood we need to know what fraction of the logging residues were utilized as bioenergy feedstock. McIver and Morgan (Table 6) reports bioenergy consumption from 2000 forward. For years previous, we use the average bioenergy consumption from 2000 -- 2012.
End of explanation
"""


def WPu (rw_harvest, lr, year, mill_efficiency = constants['me']['value'], wdens = constants['wDens']['value'], df = constants['DFu']['value']):
    '''
    Calculates the emissions reduction resulting from harvested wood with with utilization of loggin residuals for bioenergy
    '''
    # establish the aproporiate bioenergy consumption, if no data on bioenergy consumption exists, use average from 2000-2012
    if year in bio_dict.keys():
        bioe_pct = bio_dict[year]['biopct']
    else:
        bioe_pct = np.average(bio_pct.biopct)
    #Calculate total volume used in bioenergy
    bioevol = bioe_pct * rw_harvest
    #Establish utilization ratio for bioenergy
    lrUsed = bioevol/lr
    #Calcuate roundwood harvest volume fromwhich loggin residues were utilized
    HWu = lrUsed * rw_harvest
    #Calculate volume of final wood product produced using mill efficiency. S&O use volume of wood product not sawlogs for DF
    WPu = HWu * mill_efficiency * wdens * 1000000 / 2204.62 * 0.5 * df
    #Per comment from Roger Sathre,needs to be reduced by 50% before applying the DF as the DF is meant for tC not tWood..
    #This is in MT
    
    return WPu
    

def WPo (rw_harvest, lr, year, mill_efficiency = constants['me']['value'], wdens = constants['wDens']['value'], df = constants['DFo']['value']):
    '''
    Calculates the emissions reduction resulting from harvested wood without utilization of 
    logging residuals for bioenergy
    '''
    # establish the aproporiate bioenergy consumption lever for a given year, if no data on bioenergy consumption exists, use average from 2000-2012
    if year in bio_dict.keys():
        bioe_pct = bio_dict[year]['biopct']
    else:
        bioe_pct = np.average(bio_pct.biopct)
    #Calculate total volume used in bioenergy
    bioevol = bioe_pct * rw_harvest
    #Establish utilization ratio for bioenergy
    lrUsed = bioevol/lr
    #Calcuate roundwood harvest volume fromwhich loggin residues were utilized
    HWo = (1-lrUsed) * rw_harvest
    #Calculate volume of final wood product produced using mill efficiency. S&O use volume of wood product not sawlogs for DF
    WPo = HWo * mill_efficiency * wdens * 1000000 / 2204.62 * 0.5 * df
    #This is in MT
    
    return WPo

"""
Explanation: Logging residuals
The BOE data does not specifically estimate logging residuals, it simply reports harvested roundwood. To accurately ascribe fate to roundwood harvested, an estimation of logging residuals must be made
Calculating emissions reductions
The following functions calculate the displaced emissions resulting from wood harvested with and without logging residue utilization. They return estimates in metric tons of CO2 equivalents.
End of explanation
"""


erWPu = []
for row in tpoData.index:
    rw,lr,yr = tpoData.iloc[row][['roundwoodproducts','loggingresidues', 'year']].tolist()
    erWPu.append(WPu(rw,lr,yr))
tpoData['erWPu'] = erWPu

erWPo = []
for row in tpoData.index:
    rw,lr,yr = tpoData.iloc[row][['roundwoodproducts','loggingresidues', 'year']].tolist()
    erWPo.append(WPo(rw,lr,yr))
tpoData['erWPo'] = erWPo

"""
Explanation: Emissions reduction from harvested wood with LR utilized
Emissions reductions resulting from harvested roundwood with logging residue utilized in bioenergy
End of explanation
"""


tpoData['erTotal'] = tpoData.erWPo+tpoData.erWPu
tpoData.to_sql('tpo_emreduc', sqdb['cx'], if_exists='replace')
tpoData['bioe_pct'] = tpoData.year.apply(bioPct)
tpoData['bioe_t'] = tpoData.bioe_pct * tpoData.loggingresidues * 1e6* constants['wDens']['value']/2204.62

"""
Explanation: Emissions reduction from harvested wood without LR utilization
Emissions reductions resulting from harvested roundwood without logging residue utilized in bioenergy. Though wood with LR utilization rate has a higher displacement factor, the majority of loggin residues wer not utilized.
End of explanation
"""


erWPo = []
for row in mm_histHarvest.index:
    r = mm_histHarvest.iloc[row]
    yr = r['year'] ## year
    rw = (r.state+r.private+r.tribal+r.blm+r.nat_forest)/5.44 
    qry = 'select avg(loggingresidues/roundwoodproducts) lr from tpo where year = {}'.format(yr)
    if yr in tpoData.year.tolist():
        lr = pd.read_sql(qry, sqdb['cx'])*rw
    else:
        lr = pd.read_sql('select avg(loggingresidues/roundwoodproducts) lr from tpo', sqdb['cx'])*rw
    erWPo.append(WPo(rw,lr,yr).lr[0])
mm_histHarvest['erWPo'] = erWPo

erWPu = []
lrVect = []
tHarv = []
for row in mm_histHarvest.index:
    r = mm_histHarvest.iloc[row]
    yr = r['year'] ## year
    rw = (r.state+r.private+r.tribal+r.blm+r.nat_forest)/5.44 
    qry = 'select avg(loggingresidues/roundwoodproducts) lr from tpo where year = {}'.format(yr)
    if yr in tpoData.year.tolist():
        lr = pd.read_sql(qry, sqdb['cx'])*rw
    else:
        lr = pd.read_sql('select avg(loggingresidues/roundwoodproducts) lr from tpo', sqdb['cx'])*rw
    lrVect.append(lr.lr[0])
    tHarv.append(rw)
    erWPu.append(WPu(rw,lr,yr).lr[0])
mm_histHarvest['erWPu'] = erWPu
mm_histHarvest['loggingresidues'] = lrVect
mm_histHarvest['totalharvest'] = tHarv

"""
Explanation: Using M&M Historical data
End of explanation
"""


mm_histHarvest['erTotal'] = mm_histHarvest.erWPo+mm_histHarvest.erWPu
mm_histHarvest.to_sql('mm_emreduc', sqdb['cx'], if_exists='replace')
mm_histHarvest['bioe_pct'] = mm_histHarvest.year.apply(bioPct)
mm_histHarvest['bioe_t'] = mm_histHarvest.bioe_pct * mm_histHarvest.loggingresidues * 1e6* constants['wDens']['value']/2204.62
mm_histHarvest.to_sql()

mm_histHarvest

sns.set_style("whitegrid")
fig2, ax2 = plt.subplots(figsize=(12, 10))
ax2 = sns.barplot(x ='year', y='erTotal', data=mm_histHarvest.sort_values('year'))
ax2.set_ylabel('Emissions reduction (MT CO2e)')
ax2.set_title('Emissions reductions resulting \nfrom roundwood harvest in CA')
ax2.set_xticklabels(ax2.get_xticklabels(),rotation=90)

[fig2.savefig('graphics/ann_hh_em_reduc.{}'.format(i)) for i in ['pdf','png']]

sns.set_style("whitegrid")
fig, ax = plt.subplots()
ax = sns.barplot(x ='year', y='erTotal', hue="ownership", data=tpoData.sort_values('year'))
ax.set_ylabel('Emissions reduction (MT CO2e)')
ax.set_title('Emissions reductions resulting from roundwood harvest in CA')

[fig.savefig('graphics/harv_em_reductions.{}'.format(i)) for i in ['pdf','png']]

"""
Explanation: Total emissions reduction from harvested wood products
Sum of emissions reductions from harvested wood with and without LR utilization
End of explanation
"""


pd.read_sql('select sum("erTotal") from tpo_emreduc where year = "2012"', sqdb['cx'])

"""
Explanation: Total emissions reductions from roundwood harvesting in CA, 2012
End of explanation
"""


tName = 'cpe_allyears'
sqdb['crs'].executescript('drop table if exists {0};'.format(tName))
for y in [2000, 2005, 2010, 2012, 2015]:
    url = 'http://www.arb.ca.gov/app/emsinv/2013/emsbyeic.csv?F_YR={0}&F_DIV=0&F_SEASON=A&SP=2013&SPN=2013_Almanac&F_AREA=CA'
    df = pd.read_csv(url.format(y))
    df.to_sql(tName, sqdb['cx'], if_exists = 'append')

pmAnn = pd.read_sql('''
                        select year,
                                eicsoun,
                                "PM2_5"*365 an_pm25_av
                        from cpe_allyears
                        where eicsoun = 'FOREST MANAGEMENT';
                    ''', sqdb['cx'])
pmAnn

"""
Explanation: Emissions from un-utilized logging residuals
From logging residuals not used in bioenergy, emmisions are produced from combustion of the residual material or from decomposition of the material over time. To calculate the ratio of burned to decompsed logging residues I begin with the CARB estimate of PM2.5 produced from forest management:
End of explanation
"""


pfbec = pd.read_csv('fera_pile_cemissions.csv', header=1)
ward = ut.gData('13UQtRfNBSJ81PXxbYSnB2LrjHePNcvhJhrsxRBjHpoY', 475419971)
def sp2bio(pm, species = 'PM2.5 (tons)'):
    return pm * (pfbec[species]/pfbec['Pile Biomass (tons)'])

def bioPm(pm):
    return pm * (pfbec['Pile Biomass (tons)']/pfbec['PM2.5 (tons)'])
co2t = lambda x: sp2bio(x,'CO2 (tons)')
ch4t = lambda x: sp2bio(x,'CH4 (tons)')

pmAnn['biomass_t']=pmAnn.an_pm25_av.apply(bioPm)
pmAnn['co2_t'] = pmAnn.biomass_t.apply(co2t)
pmAnn['ch4_t'] = pmAnn.biomass_t.apply(ch4t)
pmAnn['ch4_co2e'] = pmAnn.ch4_t * 56
pmAnn['bc_co2e']= pmAnn.an_pm25_av.apply(ut.pm2bcgwpPiles)
pmAnn['t_co2e']= pmAnn.co2_t + pmAnn.ch4_co2e + pmAnn.bc_co2e

print tabulate(pmAnn[['YEAR','EICSOUN','co2_t','ch4_co2e','bc_co2e','t_co2e']], headers = ['Year','Emissions source','CO2 (t)', 'CH4 (tCO2e)', 'BC (tCO2e)', 'Pile Burn Total (tCO2e)'],tablefmt="pipe")

"""
Explanation: Estimate biomass, CO2, CH4 and BC from PM2.5
To estimate total biomass from PM2.5 I assume 90% consumption of biomass in piles and use the relationship of pile tonnage to PM emissions calculated using the Piled Fuels Biomass and Emissions Calculator provided by the Washington State Department of Natural resources. This calculator is based on the Consume fire behavior model published by the US Forest Service.
End of explanation
"""


annLrAvg = pd.read_sql('''with ann as (select sum(loggingresidues) lr
                                                        from tpo
                                                        group by year)
                                                    select avg(lr) foo
                                                    from ann;''', sqdb['cx'])['foo'][0]
pctLR_bio = (np.average(pmAnn.biomass_t)/1e6)/annLrAvg

annLrAvg

pmAnn

lr_t = 1e6*tpoData.loggingresidues*constants['wDens']['value']/2204.62
tpoData['unused_lr'] = 1e6*(lr_t-(pctLR_bio*lr_t))
tpoData['burned_lr'] = 1e6*lr_t*(np.average(pmAnn.biomass_t)/(annLrAvg*1e6))
tpoData['unburned_lr'] = (lr_t*1e6) - tpoData.bioe_t - tpoData.burned_lr
tpoData['unburned_lr_co2e'] = tpoData.unburned_lr.apply(ut.co2eDecomp)
tpoData

"""
Explanation: Estimating GHG emissions from decomposition of unitilized logging slash
To provide a full picture of the emissions from residual material produced from commercial timber harvesting in California, decomposition of unutilized logging residuals left on-site that are not burned must be accounted for. To establish the fraction of logging residue that is left to decompose, residues burned and used in bioenergy are subtracted from the total reported by the TPO:

To calculate the GHG emissions from decomposition of piles we use the following equation:
End of explanation
"""


pd.read_excel('lf/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')

"""
Explanation: Biomass residuals from non-commercial management activities
Data 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.
End of explanation
"""


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'])

"""
Explanation: Querying FACTS
The USFS reports 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 &gt;0 below -- nbr_units1 references NBR_UNITS_ACCOMPLISHED in the USFS dataset, see metadata for HFT here), 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.
End of explanation
"""


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()

"""
Explanation: Converting acres to cubic feet
FACTS 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
End of explanation
"""


wood_dens = ut.gData('138FWlGeW57MKdcz2UkWxtWV4o50SZO8sduB1R6JOFp8', 1297253755)
wavg_dens =sum(wood_dens.pct/100 * wood_dens.density_lbscuft)

"""
Explanation: Weighted average wood density
Average wood density weighted by harvested species percent. Derived from McIver and Morgan,  Table 4
End of explanation
"""


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')

"""
Explanation: 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
End of explanation
"""


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')


"""
Explanation: Estimating combined GHG and SLCP emissions from unutilized residues
Only a fraction of the
End of explanation
"""


