<jupyter_start><jupyter_text>TRI: Descriptive Analysis<NAME> 10.01.2020 **WHAT IS THE TOXIC RELEASE INVENTORY (TRI)** TRI tracks the management of certain toxic chemicals that may pose a threat to human health and the environment. U.S. facilities in different industry sectors must report annually how much of each chemical is released to the environment and/or managed through recycling, energy recovery and treatment. (A "release" of a chemical means that it is emitted to the air or water, or placed in some type of land disposal.) [1]---**WHICH CHEMICALS AND COMPANIES ARE REQUIRED TO REPORT?**- Have 10 or more full-time employee equivalents- Facility is included in the following NAICS codes:    - Manufacturing of: food, beverage, apparel, leather, wood, paper, petroleum, coal, chemical, plastics, rubber, minerals, metals, machinery, computers, electronics, applicance, transportation, furniture and more. - Facility uses chemical in quanitities greater than establihed threshold in the course of a calendar year. (changes year to year)---**WHICH CHEMICALS ARE TRACKED?**Those associated with:1. Cancer or other chronic human health effects2. Significant adverse acute human health effects3. Significant adverse environmental effectsThere are currently 755 individually listed chemicals and 33 chemical categories covered by the TRI Program. Facilities that manufacture, process or otherwise use these chemicals in amounts above established levels must submit annual reporting forms for each chemical. Note that the TRI chemical list doesn't include all toxic chemicals used in the U.S.--- **WHEN WAS THE PROGRAM STARTED**In 1986, Congress passed the Emergency Planning and Community Right-to-Know Act (EPCRA) to support and promote emergency planning and to provide the public with information about releases of toxic chemicals in their community. Section 313 of EPCRA established the Toxics Release Inventory. Data is collected on a yearly basis. [1] - https://www.epa.gov/toxics-release-inventory-tri-program/what-toxics-release-inventory[2] - https://www.epa.gov/toxics-release-inventory-tri-program/tri-covered-industry-sectors<jupyter_code>#Libraries

#Basic py: 
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import datetime
import os
import math

#Geo
import geopandas as gpd
import fiona
from shapely.geometry import Point
import descartes
import contextily as ctx #Basemaps<jupyter_output><empty_output><jupyter_text>Load the Data<jupyter_code>#Import all data and merge into one dataframe
path = "/home/boogie2/Hanson_Lab/TRI_STILT/data/raw/Toxic_Release_Inventory_raw" 
for index, files in enumerate(os.listdir(path)):
    temp_df = pd.read_csv(path + "/" + files)

    #Columns are standard, merge the dataset on that
    if index == 0:
        tri_df = temp_df
        print('check')
    else:
        tri_df = tri_df.append(temp_df)
        
#Clean up the data
tri_df.rename(columns=lambda x: ' '.join(x.split('.')[1:]), inplace=True)
tri_df.rename(columns=lambda x: x.replace(' ',''), inplace=True)

#Select the data from 1990-1999 which is air pollution
select_tri = tri_df.loc[(tri_df['YEAR']>=1990) & (tri_df["YEAR"] < 2000)]

select_tri_air = select_tri[(select_tri['51-FUGITIVEAIR']>0) | (select_tri['52-STACKAIR']>0)]
print("From 1990-1999 there are {0} data entries and {1} features".format(select_tri.shape[0],select_tri.shape[1]))<jupyter_output>check
From 1990-1999 there are 5733 data entries and 114 features<jupyter_text>The loaded shape of the data is what we expect. The 114 features represent details about the location, industry and toxin. For more information on this, please visit the EPAs website (https://www.epa.gov/toxics-release-inventory-tri-program/descriptions-tri-data-termspopup_1)  Data Cleaning<jupyter_code>#Percentage of data missing values
pd.set_option('display.max_rows', 500)
print(select_tri_air.isna().sum().sort_values()/select_tri_air.shape[0])<jupyter_output>YEAR                       0.000000
62-M26                     0.000000
62-M24                     0.000000
62-M20                     0.000000
OFF-SITERELEASETOTAL       0.000000
62-M99                     0.000000
62-M94                     0.000000
62-M90                     0.000000
62-M79                     0.000000
62-M73                     0.000000
62-M65                     0.000000
62-M64                     0.000000
62-M67                     0.000000
62-M66                     0.000000
62-M63                     0.000000
62-M72                     0.000000
62-M82                     0.000000
62-M81                     0.000000
62-M71                     0.000000
62-M61METAL                0.000000
62-M40METAL                0.000000
62-M62                     0.000000
62-M41                     0.000000
62-M10                     0.000000
62-M28                     0.000000
62-M93                     0.000000
OFF-SITERECYCLEDTOTAL      0.000000
62-M56                     0[...]<jupyter_text>Examination of the missing values showcases the highest number of missing values are found in the NAIC and SIC codes. This makes sense, as not all companies had multiple codes. Additionally, approximately 23% of the data is missing parent company data and binary release data. We will only keep data for which there is less than 20% of missing data.<jupyter_code>clean_index = select_tri_air.isna().sum()<math.floor(select_tri_air.shape[0]*0.20)
select_tri_air_clean = select_tri_air.loc[:,clean_index.index[clean_index==True]]

#Let's collect and analyze a subset of the data which focuses on TRI air releases
select_tri_air_clean = select_tri_air_clean[['YEAR','TRIFD','FRSID','FACILITYNAME','CITY','COUNTY','ST','ZIP','LATITUDE','LONGITUDE',
                           'INDUSTRYSECTORCODE','INDUSTRYSECTOR','CHEMICAL','CAS#/COMPOUNDID','METAL','CARCINOGEN' ,
                           'UNITOFMEASURE','51-FUGITIVEAIR','52-STACKAIR','PRIMARYSIC','INDUSTRYSECTORCODE',
                                             '88-ONE-TIMERELEASE','PRODUCTIONWSTE(81-87)']]<jupyter_output><empty_output><jupyter_text>Selecting IARC Labelled "Known" CarcinogensThe international agency for research on cancer created a listing of carcinogens with varying confidence. While the TRI data contains a carcinogen column, we decided to use the IARC classifications for chemical filtering. We find this measure to be more conservative. Ratings: 1 - Known Carcinogen2A - Probably carcinogenic to humans2B - Possibly carcinogenic to humans3 - not classifiable as to its carcinogenicity to humans<jupyter_code>#Import the IARC Data
iarc_df = pd.read_csv('/home/boogie2/Hanson_Lab/TRI_STILT/data/raw/IARC_Class_Full_List.csv')

#Standardize the Column Names
iarc_df.rename(columns=lambda x: x.replace(' ','_'), inplace=True)
iarc_df.rename(columns=lambda x: x.replace('.',''), inplace=True)

#The CAS numbers are not standardized between groups. Let's change that!
#A CASRN contains a first number (2-7 digits) - second number (2 digits) - third number - 1 digit=
iarc_df.CAS_No = iarc_df.CAS_No.str.replace('-','').apply('{:0>9}'.format)
iarc_df.CAS_No = iarc_df.CAS_No.str.replace('000000nan',str(np.nan))

print("There are {0} chemicals with {1} features".format(iarc_df.shape[0],iarc_df.shape[1]))
#Select those chemicals with a 1,2A or 2B designation (KNOWN CARCINOGENS)
known_carc = iarc_df[iarc_df.Group != "3"]
#Merge the data
carc_sel_tri_air = known_carc.merge(select_tri_air_clean, 
                                    left_on='CAS_No', 
                                    right_on='CAS#/COMPOUNDID',
                                    how='inner')

print("Data are present from {0} to {1}".format(carc_sel_tri_air['YEAR'].min(),carc_sel_tri_air['YEAR'].max()))
print("There are {0} total observations in the data".format(carc_sel_tri_air.shape[0]))
#Fix those entries which have a "see xx (other chemical)"
carc_sel_tri_air.loc[carc_sel_tri_air.Agent== "Strong-inorganic-acid mists containing sulfuric acid (see Acid mists)",
       'Group'] = '1'

carc_sel_tri_air.loc[carc_sel_tri_air.Agent== "Bis(2-ethylhexyl) phthalate (see Di(2-ethylhexyl) phthalate)",
       'Group'] = '2B'

#Clean up unnessessary/unusable information
carc_sel_tri_air = carc_sel_tri_air.drop(columns=['CAS_No','Agent','Volume', 'Year', 'Additional_information'])
#Histogram: total number of chemicals per each IARC bin 
fig,ax = plt.subplots(1,1,figsize=(12,6))
ax = carc_sel_tri_air['Group'].value_counts().plot.bar()
ax.set_xlabel('IARC Classification')
ax.set_ylabel('Observations')
ax.set_title('Histogram of TRI Carcinogenic Status')

plt.close()<jupyter_output><empty_output><jupyter_text>Descriptive Analysis<jupyter_code>#Unique locations and chemicals
uniq_locs = carc_sel_tri_air.groupby(['LATITUDE', 'LONGITUDE']).size().reset_index(name='Freq')
uniq_chems = carc_sel_tri_air.groupby(['CAS#/COMPOUNDID']).size().reset_index(name='Freq')
print("There are a total of {0} unique locations".format(uniq_locs.shape[0]))
print("There are a total of {0} unique chemicals".format(uniq_chems.shape[0]))
#What are the most and least common chemicals
agg_chems =carc_sel_tri_air['CHEMICAL'].value_counts()
print("Most common chemical")
print(agg_chems.head(5))

print("\nLeast common chemicals")
print(agg_chems.tail(5))
val = carc_sel_tri_air[carc_sel_tri_air['51-FUGITIVEAIR']>0].shape[0]/carc_sel_tri_air.shape[0]
val2 = carc_sel_tri_air[carc_sel_tri_air['52-STACKAIR']>0].shape[0]/carc_sel_tri_air.shape[0]
val3 = carc_sel_tri_air[(carc_sel_tri_air['51-FUGITIVEAIR']>0) & (carc_sel_tri_air['52-STACKAIR']>0)].shape[0]/carc_sel_tri_air.shape[0]

print('Percentage of TRI Releases which have only Fugitive Release: {0:.2f}%'.format(val-val3))
print('Percentage of TRI Releases which have only Stack Release: {0:.2f}%'.format(val2-val3))
print('Percentage of TRI Releases with both Fugitive and Stack: {0:.2f}%'.format(val3))
val4 = carc_sel_tri_air[carc_sel_tri_air['88-ONE-TIMERELEASE']>0].shape[0]/carc_sel_tri_air.shape[0]
print('Percentage of TRI Releases which contain a one time release: {0:.2f}%'.format(val4*100.00))
print('Approximately 20% of the TRI values are missing this label though and remain unknown')
#Plotting the total releases per year
fig,ax = plt.subplots(1,1,figsize=(12, 7))
carc_sel_tri_air['YEAR'].value_counts().sort_index().plot.bar(ax=ax)
ax.set_xticklabels(sorted(carc_sel_tri_air['YEAR'].value_counts().index), rotation = 45, ha="right")
ax.set_xlabel('Year')
ax.set_ylabel('Observations')
ax.set_title('TRI Occurences per Year')


plt.savefig('/home/boogie2/Hanson_Lab/TRI_STILT//figures/descriptive_analysis/Average_TRI.png')
plt.close()<jupyter_output><empty_output><jupyter_text>Visualizing Spatial Distributions<jupyter_code>#Where are the TRI Releases?
gdf = gpd.GeoDataFrame(carc_sel_tri_air, geometry=gpd.points_from_xy(carc_sel_tri_air.LONGITUDE, carc_sel_tri_air.LATITUDE),crs='epsg:4326')
gdf = gdf.to_crs(3857)

fig, ax = plt.subplots(figsize = (10,10)) 
gdf.plot(ax = ax,markersize=20, marker='^', color='purple')

#Add Title
plt.title("Reporting TRI Facilities from 1990-1999")

#Add Legend
plt.legend(['TRI Reporting Facility Locations'], loc='lower right')

#Add Basemap
ctx.add_basemap(ax)

#Clean Up the Plot
ax.set_axis_off()
fig.patch.set_facecolor('xkcd:white')

#Save the figure
plt.savefig('/home/boogie2/Hanson_Lab/TRI_STILT//figures/descriptive_analysis/TRI_Locations.png')
plt.close()
#Load the County Shapefiles
counties =fiona.open('/home/boogie2/Hanson_Lab/TRI_STILT/data/raw/utah_counties.gdb')
counties_gdf = gpd.GeoDataFrame.from_features([feature for feature in counties], crs=counties.crs)
columns = list(counties.meta["schema"]["properties"]) + ["geometry"]
counties_gdf = counties_gdf[columns]
counties_gdf = counties_gdf.to_crs(epsg=3857) 

#Get the total number per country
temp=carc_sel_tri_air
temp['sumer'] = 1
temp = temp.groupby(['COUNTY']).agg({'sumer': 'sum'}).reset_index()

#Merge on County Name
merge = temp.merge(counties_gdf, 
                   left_on='COUNTY',
                   right_on='NAME',
                   how='inner')

#Place back into a geodataframe
merge = gpd.GeoDataFrame(merge)


#Check the shape match. The rows should not change
print(temp.shape[0] == merge.shape[0])
#Which Counties have the most observations? 
fig,ax = plt.subplots(1,2,figsize=(16, 8))
carc_sel_tri_air['COUNTY'].value_counts().plot.bar(ax=ax[0])
ax[0].set_xticklabels(carc_sel_tri_air['COUNTY'].value_counts().index, rotation = 45, ha="right")
ax[0].set_xlabel('County')
ax[0].set_ylabel('Observations')
ax[0].set_title('TRI Occurences by County')

#Geographic Plotting
merge.plot(cmap='YlOrRd',
           ax=ax[1],
           column='sumer',
           legend=True,
           alpha=0.8,
           legend_kwds={'label': "Occurences within County Limits"})
ax[1].set_title("Utah TRI Occurence Locations")

#Pulling off the axis
plt.rcParams.update({'font.size': 12})
plt.tick_params(
    axis='both',          # changes apply to the x-axis
    which='both',      # both major and minor ticks are affected
    bottom=False,      # ticks along the bottom edge are off
    top=False,# ticks along the top edge are off
    left=False,
    labelbottom=False,
    labelleft=False)

ctx.add_basemap(ax[1])

plt.savefig('/home/boogie2/Hanson_Lab/TRI_STILT//figures/descriptive_analysis/TRI_county_occurences.png')
plt.close()
#Per year, how do the pounds of released material change? 

fig,ax = plt.subplots(1,2,figsize=(16, 8))

data = []
for year in sorted(carc_sel_tri_air['YEAR'].drop_duplicates()):
    d=carc_sel_tri_air[carc_sel_tri_air['YEAR']==year]['51-FUGITIVEAIR']
    data.append(d)
    
ax[0].boxplot(data)
ax[0].set_xticklabels(sorted(carc_sel_tri_air['YEAR'].drop_duplicates()), rotation = 45, ha="right")
ax[0].set_xlabel('Year')
ax[0].set_ylabel('Lbs of Released Chemicals')
ax[0].set_title('Yearly Fugitive Releases')

data = []
for year in sorted(carc_sel_tri_air['YEAR'].drop_duplicates()):
    d=carc_sel_tri_air[carc_sel_tri_air['YEAR']==year]['52-STACKAIR']
    data.append(d)
    
ax[1].boxplot(data)
ax[1].set_xticklabels(sorted(carc_sel_tri_air['YEAR'].drop_duplicates()), rotation = 45, ha="right")
ax[1].set_xlabel('County')
ax[1].set_ylabel('Lbs of Released Chemicals')
ax[1].set_title('Yearly Stack Releases')

plt.savefig('/home/boogie2/Hanson_Lab/TRI_STILT//figures/descriptive_analysis/Stack_Fug_boxplot_yearly_releases.png')
plt.close()
#Per county are the releases centered in any sort of manner?
fig,ax = plt.subplots(1,2,figsize=(16, 8))

data = []
for year in sorted(carc_sel_tri_air['COUNTY'].drop_duplicates()):
    d=carc_sel_tri_air[carc_sel_tri_air['COUNTY']==year]['51-FUGITIVEAIR']
    data.append(d)
    
ax[0].boxplot(data)
ax[0].set_xticklabels(sorted(carc_sel_tri_air['COUNTY'].drop_duplicates()), rotation = 45, ha="right")
ax[0].set_xlabel('County')
ax[0].set_ylabel('Lbs of Released Chemicals')
ax[0].set_title('County Fugitive Releases')

data = []
for year in sorted(carc_sel_tri_air['COUNTY'].drop_duplicates()):
    d=carc_sel_tri_air[carc_sel_tri_air['COUNTY']==year]['52-STACKAIR']
    data.append(d)
    
ax[1].boxplot(data)
ax[1].set_xticklabels(sorted(carc_sel_tri_air['COUNTY'].drop_duplicates()), rotation = 45, ha="right")
ax[1].set_xlabel('County')
ax[1].set_ylabel('Lbs of Released Chemicals')
ax[1].set_title('County Stack Releases')


plt.savefig('/home/boogie2/Hanson_Lab/TRI_STILT//figures/descriptive_analysis/Stack_Fug_boxplot_county_releases.png')
plt.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12294764059469941, 'fraction_numerical': 0.038332255979314805, 'mean_word_length': 4.490063875088715, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 37, 'https://': 3, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '654869', 'license': 'MIT', 'n_tokens_mistral': 5801, 'n_tokens_neox': 5219, 'n_words': 1217, 'path': 'notebooks/descriptive_analysis.ipynb', 'repo_name': 'Hanson-Research-Lab/TRI_STILT', 'size': '15472'}