Document ID: EPA-R09-OAR-2015-0432-0203
Agency: epa
Document Type: Supporting & Related Material
Title: 
Posted Date: 2016-02-09T05:00Z

""" Python script to construct Carrying Capacity Diagrams for the Scaled-RRF
    method for the 2015 Annual PM2.5 SIP update.

    This works for old 2008 Annual Averge worksheets only.

    Usage:  python "this file".py
"""

import sys
import win32com.client
import os
import glob
import shutil
import csv

#
def scale(ws2, OC_scale, EC_scale, OTR_scale, NOx_scale, NH4_scale, SOx_scale):
    index = [ 4, 5, 6, 7 ]
    for k in index:
        ws2.Range("I" + str(k)).Value = (1. - ws2.Range("B" + str(k)).Value) * NH4_scale
        ws2.Range("J" + str(k)).Value = (1. - ws2.Range("C" + str(k)).Value) * NOx_scale
        ws2.Range("K" + str(k)).Value = (1. - ws2.Range("D" + str(k)).Value) * SOx_scale
        ws2.Range("L" + str(k)).Value = (1. - ws2.Range("E" + str(k)).Value) * EC_scale
        ws2.Range("M" + str(k)).Value = (1. - ws2.Range("F" + str(k)).Value) * OC_scale
        ws2.Range("N" + str(k)).Value = (1. - ws2.Range("G" + str(k)).Value) * OTR_scale
# Substact from unity and copy back.
        ws2.Range("B" + str(k)).Value = 1. - ws2.Range("I" + str(k)).Value
        ws2.Range("C" + str(k)).Value = 1. - ws2.Range("J" + str(k)).Value
        ws2.Range("D" + str(k)).Value = 1. - ws2.Range("K" + str(k)).Value
        ws2.Range("E" + str(k)).Value = 1. - ws2.Range("L" + str(k)).Value
        ws2.Range("F" + str(k)).Value = 1. - ws2.Range("M" + str(k)).Value
        ws2.Range("G" + str(k)).Value = 1. - ws2.Range("N" + str(k)).Value
#
def pbw(ws):
    Row_List = [14, 15, 16, 17, 24, 25, 26, 27]
    mystring_1 = "temp 1.0 1.0 1 0 rh 0. ammonium 0. sulphate nitrate 0. 0. 0. 0. 3 3 3 3 3 \
30 1 3 2 3 3 3 4 3 5 3 6 3 7 3 8 3 9 3 10 3 \
11 3 12 3 13 3 14 3 15 3 16 3 17 3 18 3 19 3 20 3 \
21 3 22 3 23 3 24 3 25 3 26 3 27 3 28 3 29 3 30 3 "
# temp,rh are fixed.
    temp,rh = 294.15, 0.35
    mystring_2 = mystring_1.replace("temp", str(temp))
    mystring_3 = mystring_2.replace("rh", str(rh))
#
    Values = []
    pbw = []
#
    for i in Row_List:
        cellnum = "%s" % (i)
        nitrate = ws.Range("C" + cellnum).Value/62.0
        sulphate = ws.Range("D" + cellnum).Value/96.0
        ammonium = (nitrate + 2.0*sulphate)
        mystring_4 = mystring_3.replace("nitrate", str(nitrate))
        mystring_5 = mystring_4.replace("sulphate", str(sulphate))
        mystring_6 = mystring_5.replace("ammonium", str(ammonium))
        Values.append([mystring_6])
#
    shutil.copy('eaim.dat.tmp', 'eaim.dat')
    f = open("eaim.dat", 'a')
    f.write("\n".join((item[0] for item in Values)))
    f.write("\n")
    f.close()
#
    os.system("eaim.exe")
#
# Search the output for (PBW in grams) and save it to PBW.
    fout = open("eaim.rs1", 'r')
    for line in fout:
        if line.find("H2O(aq)") > -1:
            pbw.append(float(line.split()[2]))
    fout.close()
#
    k = 0
    for j in Row_List:
        cellnum = "%s" % (j)
        ws.Range("H" + cellnum).Value = pbw[k]
        k = k + 1

def main():

# Open each book, loop over the diagram(copy the appropriate sheet to a
# scratch sheet, scale RRFs, calculate PBW, copy the DV to another book,
# delete the sheet), and close the book.  Make sure the book does not need any
# intervension when you open the book or copy sheets.

# Initialize Excel application.

        xl = win32com.client.gencache.EnsureDispatch("Excel.Application")

# Visible make the workbook visible.  This allows updates and saves later.

        xl.Visible = 1

# win32.com looks for the file in "My Documents".  The following few lines
# tells win32.com where the file is.  

        cwd = os.getcwd()
        print "%s" % cwd

# 2012 & 2020 emissions - With external adjustments - Baseline.
# From Steve Zelinka on Feb 13, 2015 (V1.01) + rules/incentives from Patricia
# on January 15, 2015.

        PM25_base = 65.956	
        PM25_future = 60.181

# PM is now split in to EC, OC, and OTR (=PM-EC-OC).

        EC_base = 9.711	
        EC_future = 6.498

        OC_base = 12.108	
        OC_future = 11.013

        OTR_base = 44.137	
        OTR_future = 43.306

        NOx_base = 332.191  # With Madera update
        NOx_future = 206.432
        
        SOx_base = 8.085
        SOx_future = 7.849
        
# 2014-2005 Baseline from the old Plan.

        PM25_2008_base = 86.0
        PM25_2008 = 22.7

        EC_2008_base = 11.4
        EC_2008 = 3.9

        OC_2008_base = 26.1
        OC_2008 = 7.5

        OTR_2008_base = 48.5
        OTR_2008 = 11.2

        NOx_2008_base = 575.4
        NOx_2008 = 284.2
        
        SOx_2008_base = 26.4
        SOx_2008 = 1.8

# SOx is not changed in the Carrying Capacity Diagram.  So, we can
# calculate that scaling factor here.

        SOx_scale = ((SOx_base - SOx_future)/SOx_base)*(SOx_2008_base/SOx_2008)
        print "SOx_scale %s" % SOx_scale
        
# Sheets in the Carrying Capacity Diagram.
        Scale_Sheets = ['BAC-CTRL_2013', \
                      'BEP-CTRL_2013', \
                      'FSF-CTRL_2013', 'MAD-CTRL_2013', 'MAD-CTRL_2014', \
                      'FSH-CTRL_2013', 'TRA-CTRL_2013', \
                      'CLO-CTRL_2013', \
                      'M14-CTRL_2013', 'MAN-CTRL_2013', 'TUR-CTRL_2013', \
                      'MRM-CTRL_2013', \
                      'MRC-CTRL_2013', \
                      'SOH-CTRL_2013', \
                      'VCS-CTRL_2013', 'HAN-CTRL_2013'  ]
        
# Design Values(DVs)are written into a csv file.
        dv_file = "CC_Diagrams_percent_2020.csv"
        dvWriter = csv.writer(open(dv_file, 'wb'), delimiter=',')
#
        files = glob.glob('4K-SENS_2013\*.xlsm')
        print "%s" % files
        for f in files:
            print "In the files loop"
            inname = f
            print "%s" % inname
# Open the workbooks and select the worksheets.
            wb = xl.Workbooks.Open(cwd + "/" + inname)
            number_sheets = wb.Sheets.Count
            print "Number of sheets %s" % number_sheets
            for ns in range(number_sheets):
                print "Sheet number %s" % ns
                wsname = wb.Sheets(ns+1).Name
                if wsname in Scale_Sheets:
                    print "In the worksheet %s" % wsname
                    dvWriter.writerow(wsname)
                    ws = wb.Sheets(ns+1)
                    for pm in range(30, 101, 5):
#                    for pm in range(100, 101, 5):
# DV array is used to write DVs into the csv file.
                        DV = []
# We are scaling from the 2013 future value.
#                       OC_scale = EC_scale = OTR_scale = \
#                       ((PM25_base - PM25_future * pm/100.)/PM25_base)*(PM25_2008_base/PM25_2008)
                        OC_scale = \
                        ((OC_base - OC_future * pm/100.)/OC_base)*(OC_2008_base/OC_2008)
                        EC_scale = \
                        ((EC_base - EC_future * pm/100.)/EC_base)*(EC_2008_base/EC_2008)
                        OTR_scale = \
                        ((OTR_base - OTR_future * pm/100.)/OTR_base)*(OTR_2008_base/OTR_2008)
                        print "OC_scale %s" % OC_scale
                        print "EC_scale %s" % EC_scale
                        print "OTR_scale %s" % OTR_scale
                        for nox in range(30, 101, 5):
#                        for nox in range(100, 101, 5):
                            NOx_scale = NH4_scale = \
                                ((NOx_base - NOx_future * nox/100.)/NOx_base)*(NOx_2008_base/NOx_2008)
                            print "NOx_scale %s" % NOx_scale
                            wb.Sheets(ns+1).Copy(After=wb.Sheets(number_sheets))
                            ws2 = wb.Sheets(number_sheets+1)
                            ws2.Name = wsname + "_CC"
# Scale the RRFs.  We could just pass two scaling factors, but we pass all
# just in case they were made different in the future.
                            scale(ws2, OC_scale, EC_scale, OTR_scale, NOx_scale, NH4_scale, \
                                  SOx_scale)
# Calculate particle bound water (PBW).
                            pbw(ws2)
# Write the future DV onto a different sheet.
                            DV_value = ws2.Range("O28").Value
                            DV.append(DV_value)
                            xl.DisplayAlerts=False
                            ws2.Delete()
                            xl.DisplayAlerts=True
                        dvWriter.writerow(DV)

            wb.Close(SaveChanges=1)
#
if __name__ == "__main__":
    sys.exit(main())