# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from IPython.display import IFrame
from sklearn.model_selection import train_test_split
import numpy as np
from rasterstats import point_query
import rasterio
from mpl_toolkits.axes_grid1 import make_axes_locatable
import category_encoders as ce
from sklearn.ensemble import RandomForestClassifier
from sklearn.impute import SimpleImputer
from sklearn.pipeline import make_pipeline
from sklearn.metrics import classification_report


# +
DATA_PATH = 'https://raw.githubusercontent.com/LambdaSchool/DS-Unit-2-Kaggle-Challenge/master/data/'

train = pd.merge(pd.read_csv(DATA_PATH+'waterpumps/train_features.csv'), 
                 pd.read_csv(DATA_PATH+'waterpumps/train_labels.csv'))
test = pd.read_csv(DATA_PATH+'waterpumps/test_features.csv')


# +
def wrangle(X):
    """Wrangle train, validate, and test sets in the same way"""
    
    # Prevent SettingWithCopyWarning
    X = X.copy()
    
    # About 3% of the time, latitude has small values near zero,
    # outside Tanzania, so we'll treat these values like zero.
    X['latitude'] = X['latitude'].replace(-2e-08, 0)
    
    # When columns have zeros and shouldn't, they are like null values.
    # So we will replace the zeros with nulls, and impute missing values later.
    # Also create a "missing indicator" column, because the fact that
    # values are missing may be a predictive signal.
    cols_with_zeros = ['longitude', 'latitude', 'construction_year', 
                       'gps_height', 'population']
    for col in cols_with_zeros:
        X[col] = X[col].replace(0, np.nan)
        X[col+'_MISSING'] = X[col].isnull()
            
    # Drop duplicate columns
    duplicates = ['quantity_group', 'payment_type']
    X = X.drop(columns=duplicates)
    
    # Drop recorded_by (never varies) and id (always varies, random)
    unusable_variance = ['recorded_by']
    X = X.drop(columns=unusable_variance)
    
    # Convert date_recorded to datetime
    X['date_recorded'] = pd.to_datetime(X['date_recorded'], infer_datetime_format=True)
    
    # Extract components from date_recorded, then drop the original column
    X['year_recorded'] = X['date_recorded'].dt.year
    X['month_recorded'] = X['date_recorded'].dt.month
    X['day_recorded'] = X['date_recorded'].dt.day
    X = X.drop(columns='date_recorded')
    
    # Engineer feature: how many years from construction_year to date_recorded
    X['years'] = X['year_recorded'] - X['construction_year']
    X['years_MISSING'] = X['years'].isnull()
    
    
    # return the wrangled dataframe
    return X

def convert_to_gdf(df):
    gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.longitude, df.latitude))
    gdf.crs = "EPSG:4326"
    return gdf


# +
train = wrangle(train)
test = wrangle(test)

train_gdf = convert_to_gdf(train)
test_gdf = convert_to_gdf(test)

# +
fig, ax = plt.subplots(figsize=(10,10))

train_gdf.plot(ax=ax, color="green")
test_gdf.plot(ax=ax, color="orange")

plt.show()
# -

IFrame("https://data.humdata.org/dataset/2012-census-tanzania-wards-shapefiles", width="95%", height=400)

wards_data = "/home/alex/data/2012 Wards Shapefiles/TZwards.shp"
wards_gdf = gpd.read_file(wards_data)
wards_gdf.head()

fig, ax = plt.subplots(figsize=(10,10))
wards_gdf.plot("District_N", ax=ax)
plt.show()

wards_gdf_prj = wards_gdf.to_crs("EPSG:4326")

# +
fig, ax = plt.subplots(figsize=(10,10))

wards_gdf_prj.plot(ax=ax, color="black")
train_gdf.plot(ax=ax, color="blue", markersize=1)
test_gdf.plot(ax=ax, color="blue", markersize=1)

plt.show()

# +
train_join = gpd.sjoin(train_gdf, wards_gdf_prj, how="left")
test_join = gpd.sjoin(test_gdf, wards_gdf_prj, how="left")

train_ward_key = train_join[["id","Region_Nam","District_N","Ward_Name","Division"]]
test_ward_key = test_join[["id","Region_Nam","District_N","Ward_Name","Division"]]

# -

IFrame("http://opendata.rcmrd.org/datasets/tanzania-srtm-dem-30-meters", width="95%", height=500)

tif = "/home/alex/data/Tanzania_SRTM30meters/Tanzania_SRTM30meters.tif"
rasterio.open(tif).crs

train_for_pq = train_gdf[train_gdf["latitude"].notna()].copy()
test_for_pq = test_gdf[test_gdf["latitude"].notna()].copy()

train_elev = point_query(train_for_pq, tif, interpolate="nearest")
test_elev = point_query(test_for_pq, tif, interpolate="nearest")

train_for_pq["SRTM_elevation"] = train_elev
test_for_pq["SRTM_elevation"] = test_elev

fig, ax = plt.subplots(figsize=(8,6))
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
train_for_pq.plot("SRTM_elevation", ax=ax, cax=cax, cmap="gist_earth", legend=True)
plt.suptitle("Water Pump Elevations in Tanzania", size="xx-large", y=0.95)
plt.show()

# +
train_elev_key = train_for_pq[["id","SRTM_elevation"]].copy()
test_elev_key = test_for_pq[["id","SRTM_elevation"]].copy()

train_elev_key.head()
# -

train = train.merge(train_ward_key, on="id", how="left").merge(train_elev_key, on="id", how="left").drop(["geometry"], axis=1)
test = test.merge(test_ward_key, on="id", how="left").merge(test_elev_key, on="id", how="left").drop(["geometry"], axis=1)

train.head()

IFrame("https://data.humdata.org/dataset/413174ce-fc61-42eb-99e5-41a96c9fda2f", width="95%", height=400)

ward_population = pd.read_excel('/home/alex/Downloads/tza_pop_popn_nbs_baselinedata.xlsx')
ward_population.head()

ward_population_clean = ward_population[["Ward_Name","total_both","total_male","total_fema","rural_both","rural_male","rural_fema", "urban_both","urban_male","urban_fema","ward_type","Pop_Den"]]
ward_population_clean

train_prepped = train.merge(ward_population_clean, on="Ward_Name", how="left").drop_duplicates(["id"]).copy()
test_prepped = test.merge(ward_population_clean, on="Ward_Name", how="left").drop_duplicates(["id"]).copy()

train_prepped, validate_prepped = train_test_split(
    train_prepped, 
    train_size=0.80, 
    test_size=0.20, 
    stratify=train['status_group'], 
    random_state=42)

train_prepped.shape, validate_prepped.shape, test_prepped.shape

# +
# The status_group column is the target
target = 'status_group'

# Get a dataframe with all train columns except the target
train_features = train_prepped.drop(columns=[target, "id"])

# Get a list of the numeric features
numeric_features = train_features.select_dtypes(include='number').columns.tolist()

# Get a series with the cardinality of the nonnumeric features
cardinality = train_features.select_dtypes(exclude='number').nunique()

# Get a list of all categorical features
categorical_features = cardinality.index.tolist()

# Combine the lists 
features = numeric_features + categorical_features

# +
X_train = train_prepped[features]
X_validate = validate_prepped[features]
X_test = test_prepped[features]

y_train = train_prepped[target]
y_validate = validate_prepped[target]

# +
pipeline = make_pipeline(
    ce.OrdinalEncoder(), 
    SimpleImputer(), 
    RandomForestClassifier(random_state=8, min_samples_leaf=2, max_depth=32, n_estimators=370, n_jobs=-1))
    
pipeline.fit(X_train, y_train)

print(pipeline.score(X_train, y_train))
print(pipeline.score(X_validate, y_validate))
# -

y_pred = pipeline.predict(X_validate)
print(classification_report(y_validate, y_pred))

# +
model = pipeline.named_steps["randomforestclassifier"]

encoder = pipeline.named_steps["ordinalencoder"]

encoded_columns = encoder.transform(X_validate).columns

importances = pd.Series(model.feature_importances_, encoded_columns)

plt.figure(figsize=(10,15))
importances.sort_values().plot(kind="barh", color="green")
plt.title("Random Forest Classifier Feature Importances", size="xx-large")
plt.show()

StarcoderdataJupyterScriptsDedupFiltered