# Import libraries import pandas as pd import numpy as np import random import os from tqdm.notebook import tqdm import geopandas as gpd from shapely.geometry import Point import folium import matplotlib.pyplot as plt import seaborn as sns from sklearn.ensemble import RandomForestRegressor from sklearn.model_selection import train_test_split from sklearn.metrics import mean_squared_error pd.options.display.float_format = '{:.5f}'.format pd.options.display.max_rows = None %matplotlib inline import warnings warnings.filterwarnings('ignore')

# Set seed for reproducability SEED = 2023 random.seed(SEED) np.random.seed(SEED)

# Load files train = pd.read_csv('/datasets/ey-carbon-dioxide-data/Train.csv') test = pd.read_csv('/datasets/ey-carbon-dioxide-data/Test.csv') samplesubmission = pd.read_csv('/datasets/ey-carbon-dioxide-data/SampleSubmission.csv') # Preview train dataset train.head()

test.head()

samplesubmission.head()

# Train to test sets ratio (test.shape[0]) / (train.shape[0] + test.shape[0])

# Train statistical summary train.describe(include = 'all')

# Check size and shape of datasets train.shape, test.shape, samplesubmission.shape

# Target variable distribution sns.set_style('darkgrid') plt.figure(figsize = (13, 7)) sns.histplot(train.emission, kde = True, bins = 15) plt.title('Target variable distribution', y = 1.02, fontsize = 15) display(plt.show(), train.emission.skew())

# Plotting boxplot for the CO2 emissions sns.set_style('darkgrid') plt.figure(figsize = (13, 7)) sns.boxplot(train.emission) plt.title('Boxplot showing CO2 emission outliers', y = 1.02, fontsize = 15) plt.show()

# Combine train and test for easy visualisation train_coords = train.drop_duplicates(subset = ['latitude', 'longitude']) test_coords = test.drop_duplicates(subset = ['latitude', 'longitude']) train_coords['set_type'], test_coords['set_type'] = 'train', 'test' all_data = pd.concat([train_coords, test_coords], ignore_index = True) # Create point geometries geometry = gpd.points_from_xy(all_data.longitude, all_data.latitude) geo_df = gpd.GeoDataFrame( all_data[["latitude", "longitude", "set_type"]], geometry=geometry ) # Preview the geopandas df geo_df.head()

# Create a canvas to plot your map on all_data_map = folium.Map(prefer_canvas=True) # Create a geometry list from the GeoDataFrame geo_df_list = [[point.xy[1][0], point.xy[0][0]] for point in geo_df.geometry] # Iterate through list and add a marker for each volcano, color-coded by its type. i = 0 for coordinates in geo_df_list: # assign a color marker for the type set if geo_df.set_type[i] == "train": type_color = "green" elif geo_df.set_type[i] == "test": type_color = "orange" # Place the markers all_data_map.add_child( folium.CircleMarker( location=coordinates, radius = 1, weight = 4, zoom =10, popup= "Set: " + str(geo_df.set_type[i]) + "<br>" "Coordinates: " + str([round(x, 2) for x in geo_df_list[i]]), color = type_color), ) i = i + 1 all_data_map.fit_bounds(all_data_map.get_bounds()) all_data_map

# Check for missing values train.isnull().sum().any(), test.isnull().sum().any()

# Plot missing values in train set ax = train.isna().sum().sort_values(ascending = False)[:15][::-1].plot(kind = 'barh', figsize = (9, 10)) plt.title('Percentage of Missing Values Per Column in Train Set', fontdict={'size':15}) for p in ax.patches: percentage ='{:,.0f}%'.format((p.get_width()/train.shape[0])*100) width, height =p.get_width(),p.get_height() x=p.get_x()+width+0.02 y=p.get_y()+height/2 ax.annotate(percentage,(x,y))

# Check for duplicates train.duplicated().any(), test.duplicated().any()

# Year countplot plt.figure(figsize = (14, 7)) sns.countplot(x = 'year', data = train) plt.title('Year count plot') plt.show()

# Week countplot plt.figure(figsize = (14, 7)) sns.countplot(x = 'week_no', data = train) plt.title('Week count plot') plt.show()

train.drop_duplicates(subset = ['year', 'week_no']).groupby(['year'])[['week_no']].count()

# Top 20 correlated features to the target top20_corrs = abs(train.corr()['emission']).sort_values(ascending = False).head(20) top20_corrs

# Quantify correlations between features corr = train[list(top20_corrs.index)].corr() plt.figure(figsize = (13, 8)) sns.heatmap(corr, cmap='RdYlGn', annot = True, center = 0) plt.title('Correlogram', fontsize = 15, color = 'darkgreen') plt.show()

# Sample a unique location and visualize its emissions across the years train.latitude, train.longitude = round(train.latitude, 2), round(train.longitude, 2) sample_loc = train[(train.latitude == -23.73) & (train.longitude == 28.77)] # Plot a line plot sns.set_style('darkgrid') fig, axes = plt.subplots(nrows = 4, ncols = 1, figsize = (13, 10)) fig.suptitle('Co2 emissions for location lat -23.75 lon 28.75', y=1.02, fontsize = 15) for ax, data, year, color, in zip(axes.flatten(), sample_loc, sample_loc.year.unique(), ['#882255','#332288', '#999933' , 'orangered']): df = sample_loc[sample_loc.year == year] sns.lineplot(df.emission, ax = ax, label = year, color = color) plt.legend() plt.tight_layout()

# Examples of feature engineering - Aggregations, cumulative differences, moving averages ... # Lets explore the rolling mean # First we create a unique location from lat lon train['location'] = [str(x) + '_' + str(y) for x, y in zip(train.latitude, train.longitude)] # Filter based on one location example_loc = train[train.location == '-23.73_28.77'] # Calculate rolling mean for SulphurDioxide_SO2_column_number_density_amf with a window of 2 weeks rolling_mean = example_loc['SulphurDioxide_SO2_column_number_density_amf'].rolling(window = 2).mean() # Visualise rolling mean plt.figure(figsize = (15, 7)) rolling_mean.plot() plt.title('Rolling mean with a window of 2 weeks for SulphurDioxide_SO2_column_number_density_amf', y = 1.02, fontsize = 15) plt.xlabel('week', y = 1.05, fontsize = 13) plt.ylabel('SulphurDioxide_SO2_column_number_density_amf', x = 1.05, fontsize = 13) plt.show()

# Generate the above feature - rolling mean for all locations for both the train and test # Feature engineering train train_roll_mean = train.sort_values(by = ['location', 'year', 'week_no']).groupby(['location'])[train.columns[5:].tolist()].rolling(window = 2).mean().reset_index() train_roll_mean.drop(['level_1', 'emission', 'location'], axis = 1, inplace = True) train_roll_mean.columns = [col + '_roll_mean' for col in train_roll_mean.columns] # Feature engineering test test.latitude, test.longitude = round(test.latitude, 2), round(test.longitude, 2) test['location'] = [str(x) + '_' + str(y) for x, y in zip(test.latitude, test.longitude)] test_roll_mean = test.sort_values(by = ['location', 'year', 'week_no']).groupby(['location'])[test.columns[5:].tolist()].rolling(window = 2).mean().reset_index() test_roll_mean.drop(['level_1', 'location'], axis = 1, inplace = True) test_roll_mean.columns = [col + '_roll_mean' for col in test_roll_mean.columns] test_roll_mean.head()

# Merge engineered features with train and test set #Train train_eng = train.sort_values(by = ['location', 'year', 'week_no'], ignore_index = True).merge(train_roll_mean, how = 'left', left_index=True, right_index=True) # Test test_eng = test.sort_values(by = ['location', 'year', 'week_no'], ignore_index = True).merge(test_roll_mean, how = 'left', left_index=True, right_index=True) # Preview engineered test set test_eng.head()

train[train_eng.isnull()].columns

!pip install xgboost

from xgboost import XGBRegressor # Selecting the independent variables and the target variable X = train_eng.drop(['ID_LAT_LON_YEAR_WEEK', 'location', 'emission'], axis = 1).fillna(0) y = train_eng.emission # Splitting the data into training and testing sets X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.275, random_state = SEED) # Instantiating the model clf = XGBRegressor(random_state = SEED, n_jobs=-1, learning_rate=0.1) clf.fit(X_train, y_train) # Making predictions y_pred = clf.predict(X_test) # Measuring the accuracy of the model print(f'RMSE Score: {mean_squared_error(y_test, y_pred, squared=False)}')

# Selecting the independent variables and the target variable X = train_eng.drop(['ID_LAT_LON_YEAR_WEEK', 'location', 'emission'], axis = 1).fillna(0) y = train_eng.emission # Splitting the data into training and testing sets X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.3, random_state = SEED) # Instantiating the model clf = RandomForestRegressor(random_state = SEED, n_jobs=-1) clf.fit(X_train, y_train) # Making predictions y_pred = clf.predict(X_test) # Measuring the accuracy of the model print(f'RMSE Score: {mean_squared_error(y_test, y_pred, squared=False)}') # 23432.342352754695

X_test.head()

# Analyse predictions pred_errors = X_test.copy() pred_errors['emission'] = y_test pred_errors['prediction'] = y_pred pred_errors['error'] = abs(pred_errors.prediction - pred_errors.emission) pred_errors = pred_errors[['latitude', 'longitude','year','week_no', 'emission', 'prediction', 'error']] pred_errors.sort_values(by = 'error', ascending = False, inplace = True) pred_errors.head()

pred_errors.tail()

train.emission.describe()

# Feature importance impo_df = pd.DataFrame({'feature': X.columns, 'importance': clf.feature_importances_}).set_index('feature').sort_values(by = 'importance', ascending = False) impo_df = impo_df[:12].sort_values(by = 'importance', ascending = True) impo_df.plot(kind = 'barh', figsize = (10, 10)) plt.legend(loc = 'center right') plt.title('Bar chart showing feature importance', fontsize = 14) plt.xlabel('Features', fontsize = 12) plt.show()

# Make prediction on the test set test_df = test_eng.drop(['ID_LAT_LON_YEAR_WEEK', 'location'], axis = 1).fillna(0) predictions = clf.predict(test_df) # # Create a submission file sub_file = pd.DataFrame({'ID_LAT_LON_YEAR_WEEK': test_eng.ID_LAT_LON_YEAR_WEEK, 'emission': predictions}) sub_file.head()

# Create file sub_file.to_csv('BaselineSubmission.csv', index = False) # Download subfile and submit to zindi for scoring