{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd \nimport matplotlib.pyplot as plt \nimport numpy as np \nimport seaborn as sns\n\nimport plotly.express as px\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\n\nfrom tqdm.notebook import tqdm \nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.model_selection import GridSearchCV\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.preprocessing import StandardScaler\n\nimport matplotlib\nfrom matplotlib import rc\n\nrc(\"font\", **{\"family\": \"sans-serif\", \"sans-serif\": \"DejaVu Sans\"})\nrc(\"figure\", **{\"dpi\": 200})\nrc(\n    \"axes\",\n    **{\"spines.right\": False, \"spines.top\": False, \"xmargin\": 0.0, \"ymargin\": 0.05}\n)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-05-11T08:21:36.075019Z","iopub.execute_input":"2022-05-11T08:21:36.075408Z","iopub.status.idle":"2022-05-11T08:21:39.000507Z","shell.execute_reply.started":"2022-05-11T08:21:36.075302Z","shell.execute_reply":"2022-05-11T08:21:38.999775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.read_csv('../input/eurecom-aml-2022-challenge-1/public/train.csv')\ndf_test = pd.read_csv('../input/eurecom-aml-2022-challenge-1/public/test_feat.csv')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T08:21:39.002174Z","iopub.execute_input":"2022-05-11T08:21:39.002888Z","iopub.status.idle":"2022-05-11T08:23:11.265143Z","shell.execute_reply.started":"2022-05-11T08:21:39.002852Z","shell.execute_reply":"2022-05-11T08:23:11.263918Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data Exploration\n\nFirst off we take a look at some of the different available samples using some different plots.\n\n","metadata":{}},{"cell_type":"code","source":"# we only plot a subset because with too many values the plot takes too long\ngeo_plot = df[['fact_latitude', 'fact_longitude', 'fact_temperature']].head(10000)\nfig = px.scatter_geo(geo_plot, lat='fact_latitude', lon='fact_longitude', color='fact_temperature')\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T08:23:11.267132Z","iopub.execute_input":"2022-05-11T08:23:11.267404Z","iopub.status.idle":"2022-05-11T08:23:12.915993Z","shell.execute_reply.started":"2022-05-11T08:23:11.267365Z","shell.execute_reply":"2022-05-11T08:23:12.913443Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We start by manually looking at the data and building hypotheses and assumptions from that. \nWe then compare our hypotheses to a covariance analysis in comparison to the `fact_temperature` values. \n\nWe state the following points: \n* `cmc_0_0_7` probably highly correlated -> build average\n* `gfs_cloudness` many 0 values (buckets?)\n* use `fact_latitude` to identify \"heat\" region e.g. bucket into nothern/middle/southern\n* some `fact_latitude` intervals very common in data set -> reduce before training\n* all data on numeric scale (transform some to ordinal)\n* no \"errors\" \n\nThe following columns are possible candidates to be cut, looking at it manually: \n`cmc_0_1_65` - `cmc_0_1-68`, `cmc_available`, `gfs_available`, `gfs_clouds_sea`, `gfs_precipitations`, `gfs_soil_temperature_available`, `gfs_timedelta_s` (all timedeltas in seconds so rather unimportant), `wrf_available`\n\n\n","metadata":{}},{"cell_type":"code","source":"def co_variance_analysis(df, top=30):\n    cor = df[df.columns[1:]].corr()[\"fact_temperature\"]\n    fig = make_subplots(rows=1, cols=2)\n    df_plot = pd.DataFrame({'feature': cor.index.values, 'correlation': cor.values})\n    fig.add_trace(\n        go.Scatter(x=df_plot['feature'], y=df_plot['correlation'], text=f'Top {top} Correlation Values'),\n        row=1, col=1\n    )\n\n    # remove first entry which is fact_temperature itself\n    df_plot_top = df_plot.sort_values(by='correlation', ascending=False)[1:top]\n    fig.add_trace(\n        go.Scatter(x=df_plot_top['feature'], y=df_plot_top['correlation'], text=f'Top {top} Correlation Values'),\n        row=1, col=2\n    )\n    fig.show()\n    return df_plot_top\n\ncorr_features = co_variance_analysis(df)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T08:23:12.918547Z","iopub.execute_input":"2022-05-11T08:23:12.919165Z","iopub.status.idle":"2022-05-11T08:24:22.892985Z","shell.execute_reply.started":"2022-05-11T08:23:12.919091Z","shell.execute_reply":"2022-05-11T08:24:22.892165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data Preprocessing\n\nNow that we looked at the data and identified the steps we want to take we continue with that.\n\nWe bucket the geographical zones into the following: \n- The North Frigid Zone (90° N - 66°33′48.7\" N ) : 0,\n- The North Temperate Zone (66°33′48.7\" N (== 66.5) - 23°26'11.3\" N ) : 2\n- The Torrid Zone (23°26'11.3\" N - 23°26'11.3\" S) : 4\n- The South Temperate Zone (23°26'11.3\" S - 66°33'48.7\" S) : 3\n- The South Frigid Zone (66°33'48.7\" S - 90° S) : 1\n\nWe average the different `cmc` values, but give more weight to `cmc_0_0_7_1000`, because it contains the biggest differences in datapoints.\n\nWe also bucket the `gfs_cloudness` values.","metadata":{}},{"cell_type":"code","source":"def drop_features(df, selection):\n    X = df[selection].dropna().iloc[:, 1:-1].values\n    y = df[selection].dropna().iloc[:, -1].values\n    return X, y\n\n\ndef preprocess_columns(df):\n    df = df[(df['cmc_available'] == 1.0) & (df['gfs_available'] == 1.0) & (df['wrf_available'] == 1.0)]\n    \n    df['geographical_zone'] = [0 if x > 66.5 else 1 if x < -66.5 else 2 if x > 23.5 else 3 if x < 23.5 else 4 for x in\n                               df['fact_latitude']]\n    \n    df['cmc_0_0_7_1000'] = df['cmc_0_0_7_1000'] * 2\n    df['cmc_0_0_7'] = df[\n        ['cmc_0_0_7_1000', 'cmc_0_0_7_2', 'cmc_0_0_7_500', 'cmc_0_0_7_700', 'cmc_0_0_7_850', 'cmc_0_0_7_925']].mean(\n        axis=1)\n    df['gfs_cloudness_buckets'] = [0 if x == 0 else 1 if x <= 1 else 2 if x <= 2 else 3 for x in df['gfs_cloudness']]\n    return df\n","metadata":{"execution":{"iopub.status.busy":"2022-05-11T08:24:22.894272Z","iopub.execute_input":"2022-05-11T08:24:22.894509Z","iopub.status.idle":"2022-05-11T08:24:22.902734Z","shell.execute_reply.started":"2022-05-11T08:24:22.894478Z","shell.execute_reply":"2022-05-11T08:24:22.901755Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"selection = [\n    'index', 'fact_time', 'fact_latitude', 'fact_longitude',\n    'sun_elevation', 'cmc_precipitations', 'wrf_t2_interpolated',\n    'gfs_temperature_95000', 'gfs_temperature_90000', 'gfs_temperature_65000',\n    'gfs_temperature_55000', 'gfs_temperature_45000', 'gfs_temperature_35000',\n    'gfs_a_vorticity', 'gfs_humidity', 'climate_temperature', 'gfs_precipitable_water',\n    'cmc_0_1_0_0', 'cmc_0_0_6_2', 'cmc_0_0_7', 'gfs_cloudness_buckets', 'geographical_zone',\n    'fact_temperature',\n]\ndf = preprocess_columns(df)\nX, y = drop_features(df, selection)\nX_tr, X_val, y_tr, y_val = train_test_split(X, y, random_state=1, test_size=0.40)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T08:24:22.904005Z","iopub.execute_input":"2022-05-11T08:24:22.904363Z","iopub.status.idle":"2022-05-11T08:24:28.244958Z","shell.execute_reply.started":"2022-05-11T08:24:22.904326Z","shell.execute_reply":"2022-05-11T08:24:28.244231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Model\n\nAfter some initial testing we chose a basic Random Forest Regressor for our regression task.\nWe use `sklearn`'s `make_pipeline` function, to integrate a `StandardScaler` to scale our datapoints into a normalized interval and then add a `RandomForestRegressor` to find the best model. ","metadata":{}},{"cell_type":"code","source":"def try_out_model(model):\n    pipeline = make_pipeline(StandardScaler(), model)\n    pipeline.fit(X_tr, y_tr)\n    ypred_tr = pipeline.predict(X_tr)\n    ypred_val = pipeline.predict(X_val)\n    \n    valid_rmse = mean_squared_error(y_val, ypred_val, squared=False)\n    \n    print()\n    print(f'{model.__class__} | Train RMSE: {mean_squared_error(y_tr, ypred_tr, squared=False):.3f}')\n    print(f'{model.__class__} | Valid RMSE: {valid_rmse:.3f}')\n\n    return pipeline, valid_rmse\n\n\nlr, lr_rmse = try_out_model(LinearRegression())\nrfr, rfr_rmse = try_out_model(RandomForestRegressor(n_jobs=8, n_estimators=24, verbose=True))","metadata":{"execution":{"iopub.status.busy":"2022-05-11T08:24:28.24623Z","iopub.execute_input":"2022-05-11T08:24:28.246502Z","iopub.status.idle":"2022-05-11T08:30:13.183393Z","shell.execute_reply.started":"2022-05-11T08:24:28.24647Z","shell.execute_reply":"2022-05-11T08:30:13.182327Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model_plot= pd.DataFrame({\n    'model': ('Baseline', 'Random Forest Base', 'Random Forest Optimized'), \n    'rmse': (lr_rmse, rfr_rmse, 1.95253) # optimized version ran only locally, not on kaggle\n})\n\nfig = px.bar(model_plot, x='model', y='rmse')\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T08:44:29.523077Z","iopub.execute_input":"2022-05-11T08:44:29.524515Z","iopub.status.idle":"2022-05-11T08:44:29.583543Z","shell.execute_reply.started":"2022-05-11T08:44:29.524431Z","shell.execute_reply":"2022-05-11T08:44:29.582396Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Fairness Concerns: Is the amount of sensor measurements equally spread around all the different climate zones?\n\nLooking at the plots of the data we noticed, that most of the sensor measurements were taken in the \"global north\" where most of the richer countries in the world are. At first glance it seemed like countries in the \"global south\" seemed to have lesser funds for thorough measurements and thus also are underrepresented in the samples. \n\nSo we wanted to test this hypothesis. To do this we first mapped the `lat` and `lon` values to the corresponding country codes and then took a look at this distribution. \n\n","metadata":{}},{"cell_type":"code","source":"!pip install reverse_geocoder","metadata":{"execution":{"iopub.status.busy":"2022-05-11T08:30:13.265016Z","iopub.execute_input":"2022-05-11T08:30:13.26533Z","iopub.status.idle":"2022-05-11T08:30:28.92346Z","shell.execute_reply.started":"2022-05-11T08:30:13.265297Z","shell.execute_reply":"2022-05-11T08:30:28.922094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We use to library `reverse_geocoder` to map the coordinates to country codes. If the coordinates point to an unknown country, we assume thats its international territory and save it as that. ","metadata":{}},{"cell_type":"code","source":"import reverse_geocoder as rg\n\ndef to_country(rows):\n    coords = list(zip(rows['fact_latitude'], rows['fact_longitude']))\n    geos = rg.search(coords)\n    return [g['cc'] for g in geos]\n\ndf['country'] = to_country(df)\nprint(df.head(20)['country'])","metadata":{"execution":{"iopub.status.busy":"2022-05-11T08:30:28.9271Z","iopub.execute_input":"2022-05-11T08:30:28.927467Z","iopub.status.idle":"2022-05-11T08:30:34.661911Z","shell.execute_reply.started":"2022-05-11T08:30:28.9274Z","shell.execute_reply":"2022-05-11T08:30:34.660736Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now lets look at a simple box plot of the distribution. ","metadata":{}},{"cell_type":"code","source":"country_counts = df['country'].value_counts()\n\ndf_plot = pd.DataFrame({'country': country_counts.index, 'count': country_counts.values})\n\n\nfig_bar = px.bar(df_plot, x='country', y='count', title='Amount of sensor datapoints per country')\nfig_bar.show()\n\nfig_pie = px.pie(df_plot, names='country', values='count', labels=None, title='Amount of sensor datapoints per country')\nfig_pie.update_layout(showlegend=False)\nfig_pie.update_traces(textposition='inside')\nfig_pie.update_layout(uniformtext_minsize=12, uniformtext_mode='hide')\nfig_pie.show()\n\nprint('Top 10 Countries with most datapoints')\nprint(country_counts.index.values[:10])\nprint()\nprint('Bottom 10 Countries with least datapoints')\nprint(country_counts.index.values[-10:])","metadata":{"execution":{"iopub.status.busy":"2022-05-11T08:38:14.687653Z","iopub.execute_input":"2022-05-11T08:38:14.688422Z","iopub.status.idle":"2022-05-11T08:38:14.946122Z","shell.execute_reply.started":"2022-05-11T08:38:14.68837Z","shell.execute_reply":"2022-05-11T08:38:14.945169Z"},"trusted":true},"execution_count":null,"outputs":[]}]}