{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How Much Will This Field Yield?\n",
    "\n",
    "crop_yield.py\n",
    "\n",
    "Capstone project: predicting the yield of a farm plot (tonnes per hectare)\n",
    "from its crop, weather, soil and the fertiliser used.\n",
    "\n",
    "CBSE Class XII, Artificial Intelligence (843), Part C - Capstone Project.\n",
    "SDG 2 (Zero Hunger) and SDG 12 (Responsible Consumption and Production).\n",
    "\n",
    "The program follows the 10 steps of the Data Science Methodology:\n",
    "  1 Business Understanding   2 Analytic Approach   3 Data Requirements\n",
    "  4 Data Collection          5 Data Understanding  6 Data Preparation\n",
    "  7 Modelling                8 Evaluation          9 Deployment\n",
    "  10 Feedback\n",
    "\n",
    "Two regression models are compared - Linear Regression, as in the\n",
    "handbook, and a Random Forest - with a train-test split and with 5-fold\n",
    "cross validation, using MAE, MSE, RMSE and R squared.\n",
    "\n",
    "Libraries: pandas, numpy, matplotlib, scikit-learn, joblib."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1 - Business Understanding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"=\" * 64)\n",
    "print(\"HOW MUCH WILL THIS FIELD YIELD?\")\n",
    "print(\"=\" * 64)\n",
    "print(\"Problem : farmers and planners cannot tell before harvest how much\")\n",
    "print(\"          a plot will yield, so input use and storage are guesswork.\")\n",
    "print(\"Goal    : predict yield (t/ha) from crop, weather, soil and inputs.\")\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2 - Analytic Approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The answer is a number (tonnes per hectare), so this is PREDICTIVE\n",
    "# analytics using REGRESSION.\n",
    "print(\"Analytic approach : predictive - regression (the target is a number)\")\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3 - Data Requirements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.model_selection import train_test_split, KFold, cross_val_score\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score\n",
    "import joblib\n",
    "\n",
    "pd.set_option(\"display.width\", 130)\n",
    "pd.set_option(\"display.max_columns\", 20)\n",
    "\n",
    "required = [\"crop\", \"rainfall_mm\", \"temperature_c\", \"soil_ph\", \"nitrogen_kg_ha\",\n",
    "            \"fertiliser_kg_ha\", \"irrigated\", \"yield_t_ha\"]\n",
    "print(\"Columns required :\", required)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4 - Data Collection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv(\"crop_yield.csv\")\n",
    "print(\"Rows collected   :\", len(data))\n",
    "missing_columns = [c for c in required if c not in data.columns]\n",
    "print(\"Required columns missing from the file :\", missing_columns)\n",
    "print()\n",
    "print(data.head(8))\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5 - Data Understanding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"--- Step 5: Data Understanding ---\")\n",
    "print(\"Missing values in each column:\")\n",
    "print(data.isnull().sum().to_string())\n",
    "print()\n",
    "print(\"Summary statistics:\")\n",
    "print(data.describe().round(2))\n",
    "print()\n",
    "\n",
    "# Chart 1: the distribution of yield.\n",
    "plt.figure(figsize=(7, 4))\n",
    "plt.hist(data[\"yield_t_ha\"].clip(0, 10), bins=30, color=\"#59a14f\", edgecolor=\"white\")\n",
    "plt.title(\"Distribution of yield (values above 10 shown at 10)\")\n",
    "plt.xlabel(\"Yield (tonnes per hectare)\")\n",
    "plt.ylabel(\"Number of plots\")\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"chart1_yield_distribution.png\")\n",
    "plt.close()\n",
    "\n",
    "# Chart 2: fertiliser against yield.\n",
    "ok = data[(data[\"yield_t_ha\"] > 0) & (data[\"yield_t_ha\"] < 15)]\n",
    "plt.figure(figsize=(7, 5))\n",
    "plt.scatter(ok[\"fertiliser_kg_ha\"], ok[\"yield_t_ha\"], s=12, alpha=0.5, color=\"#4e79a7\")\n",
    "plt.title(\"Fertiliser used and yield\")\n",
    "plt.xlabel(\"Fertiliser (kg per hectare)\")\n",
    "plt.ylabel(\"Yield (t/ha)\")\n",
    "plt.grid(True, linestyle=\"--\", alpha=0.4)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"chart2_fertiliser_yield.png\")\n",
    "plt.close()\n",
    "\n",
    "# Chart 3: yield of each crop, as box plots.\n",
    "plt.figure(figsize=(7, 4.5))\n",
    "crops_seen = sorted(ok[\"crop\"].str.strip().str.title().unique())\n",
    "plt.boxplot([ok[ok[\"crop\"].str.strip().str.title() == c][\"yield_t_ha\"] for c in crops_seen], labels=crops_seen)\n",
    "plt.title(\"Yield of each crop\")\n",
    "plt.ylabel(\"Yield (t/ha)\")\n",
    "plt.grid(True, axis=\"y\", linestyle=\"--\", alpha=0.4)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"chart3_yield_by_crop.png\")\n",
    "plt.close()\n",
    "for c in crops_seen:\n",
    "    part = ok[ok[\"crop\"].str.strip().str.title() == c][\"yield_t_ha\"]\n",
    "    print(\"Median yield of \" + c.ljust(6), \":\", round(part.median(), 2), \"t/ha\")\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 6 - Data Preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"--- Step 6: Data Preparation ---\")\n",
    "before = len(data)\n",
    "data = data.drop_duplicates()\n",
    "print(\"Duplicate rows removed          :\", before - len(data))\n",
    "\n",
    "data[\"crop\"] = data[\"crop\"].str.strip().str.title()\n",
    "print(\"Crops after tidying the names   :\", sorted(data[\"crop\"].unique()))\n",
    "\n",
    "impossible = ((data[\"yield_t_ha\"] <= 0) | (data[\"yield_t_ha\"] > 15)).sum()\n",
    "data = data[(data[\"yield_t_ha\"] > 0) & (data[\"yield_t_ha\"] <= 15)]\n",
    "print(\"Impossible yields removed       :\", impossible)\n",
    "\n",
    "# Missing rainfall is filled with the median rainfall of the same crop:\n",
    "# rice fields get far more rain than wheat fields, so one overall median\n",
    "# would be wrong for both.\n",
    "blank_rain = data[\"rainfall_mm\"].isnull().sum()\n",
    "data[\"rainfall_mm\"] = data.groupby(\"crop\")[\"rainfall_mm\"].transform(lambda s: s.fillna(s.median()))\n",
    "print(\"Blank rainfall values filled    :\", blank_rain, \"(median of the same crop)\")\n",
    "print(\"Rows ready for modelling        :\", len(data))\n",
    "\n",
    "# Text columns become numbers. irrigated: Yes = 1, No = 0. crop: one\n",
    "# column per crop (one-hot encoding), so no crop is treated as \"bigger\".\n",
    "data[\"irrigated\"] = (data[\"irrigated\"] == \"Yes\").astype(int)\n",
    "data = pd.get_dummies(data, columns=[\"crop\"], dtype=int)\n",
    "features = [\"rainfall_mm\", \"temperature_c\", \"soil_ph\", \"nitrogen_kg_ha\", \"fertiliser_kg_ha\",\n",
    "            \"irrigated\", \"crop_Maize\", \"crop_Rice\", \"crop_Wheat\"]\n",
    "print(\"Features used                   :\", len(features))\n",
    "print()\n",
    "\n",
    "# Chart 4: correlation of every feature with yield.\n",
    "corr = data[features + [\"yield_t_ha\"]].corr()[\"yield_t_ha\"].drop(\"yield_t_ha\").sort_values()\n",
    "print(\"Correlation of each feature with yield:\")\n",
    "print(corr.round(2).to_string())\n",
    "print()\n",
    "plt.figure(figsize=(7, 4.5))\n",
    "colours = [\"#e15759\" if v < 0 else \"#59a14f\" for v in corr.values]\n",
    "plt.barh(corr.index, corr.values, color=colours)\n",
    "plt.axvline(0, color=\"#999999\")\n",
    "plt.title(\"Correlation with yield\")\n",
    "plt.xlabel(\"Correlation coefficient\")\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"chart4_correlation.png\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 7 - Modelling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"--- Step 7: Modelling ---\")\n",
    "X = data[features]\n",
    "y = data[\"yield_t_ha\"]\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, shuffle=True, random_state=7)\n",
    "print(\"Training rows :\", len(X_train))\n",
    "print(\"Testing rows  :\", len(X_test))\n",
    "\n",
    "linear = LinearRegression()\n",
    "linear.fit(X_train, y_train)\n",
    "forest = RandomForestRegressor(n_estimators=200, min_samples_leaf=3, random_state=0)\n",
    "forest.fit(X_train, y_train)\n",
    "\n",
    "print(\"Linear Regression - intercept :\", round(linear.intercept_, 3))\n",
    "print(\"Linear Regression - fertiliser coefficient :\",\n",
    "      round(linear.coef_[features.index(\"fertiliser_kg_ha\")], 4), \"t/ha per kg\")\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 8 - Evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"--- Step 8: Evaluation ---\")\n",
    "results = []\n",
    "for name, model in [[\"Linear Regression\", linear], [\"Random Forest\", forest]]:\n",
    "    pred = model.predict(X_test)\n",
    "    mae = mean_absolute_error(y_test, pred)\n",
    "    mse = mean_squared_error(y_test, pred)\n",
    "    rmse = np.sqrt(mse)\n",
    "    r2 = r2_score(y_test, pred)\n",
    "    results.append([name, mae, mse, rmse, r2])\n",
    "    print(name)\n",
    "    print(\"  \" + name + \" MAE  :\", round(mae, 3), \"t/ha\")\n",
    "    print(\"  \" + name + \" MSE  :\", round(mse, 3))\n",
    "    print(\"  \" + name + \" RMSE :\", round(rmse, 3), \"t/ha\")\n",
    "    print(\"  \" + name + \" R2   :\", round(r2, 3))\n",
    "print()\n",
    "\n",
    "# 5-fold cross validation: every row is in the test fold exactly once.\n",
    "folds = KFold(n_splits=5, shuffle=True, random_state=1)\n",
    "cv_lin = np.sqrt(-cross_val_score(LinearRegression(), X, y, cv=folds, scoring=\"neg_mean_squared_error\"))\n",
    "cv_rf = np.sqrt(-cross_val_score(RandomForestRegressor(n_estimators=200, min_samples_leaf=3, random_state=0),\n",
    "                                 X, y, cv=folds, scoring=\"neg_mean_squared_error\"))\n",
    "print(\"RMSE in each of the 5 folds (t/ha):\")\n",
    "print(\"  Fold  Linear  Forest\")\n",
    "for k in range(5):\n",
    "    print(\"  \", k + 1, \"  \", round(cv_lin[k], 3), \" \", round(cv_rf[k], 3))\n",
    "print(\"Mean cross-validated RMSE, Linear Regression :\", round(cv_lin.mean(), 3))\n",
    "print(\"Mean cross-validated RMSE, Random Forest     :\", round(cv_rf.mean(), 3))\n",
    "print()\n",
    "\n",
    "best_pred = forest.predict(X_test)\n",
    "\n",
    "# Chart 5: predicted against actual, for both models side by side.\n",
    "lin_pred = linear.predict(X_test)\n",
    "plt.figure(figsize=(10, 4.8))\n",
    "for k in range(2):\n",
    "    plt.subplot(1, 2, k + 1)\n",
    "    pred = [lin_pred, best_pred][k]\n",
    "    plt.scatter(y_test, pred, s=14, alpha=0.6, color=[\"#9c9c9c\", \"#59a14f\"][k])\n",
    "    plt.plot([1, 9], [1, 9], color=\"#e15759\", linestyle=\"--\", label=\"Perfect prediction\")\n",
    "    plt.title([\"Linear Regression\", \"Random Forest\"][k])\n",
    "    plt.xlabel(\"Actual yield (t/ha)\")\n",
    "    plt.ylabel(\"Predicted yield (t/ha)\")\n",
    "    plt.legend()\n",
    "    plt.grid(True, linestyle=\"--\", alpha=0.4)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"chart5_predicted_actual.png\")\n",
    "plt.close()\n",
    "\n",
    "# Chart 6: the residuals (actual minus predicted) of the forest.\n",
    "residuals = y_test - best_pred\n",
    "print(\"Mean residual (Random Forest)        :\", round(residuals.mean(), 3), \"t/ha\")\n",
    "print(\"Residuals within +/- 0.5 t/ha        :\", round((residuals.abs() <= 0.5).mean() * 100, 1), \"%\")\n",
    "plt.figure(figsize=(7, 4))\n",
    "plt.hist(residuals, bins=25, color=\"#59a14f\", edgecolor=\"white\")\n",
    "plt.axvline(0, color=\"#e15759\", linestyle=\"--\")\n",
    "plt.title(\"Random Forest: residuals on the test plots\")\n",
    "plt.xlabel(\"Actual minus predicted yield (t/ha)\")\n",
    "plt.ylabel(\"Number of plots\")\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"chart6_residuals.png\")\n",
    "plt.close()\n",
    "\n",
    "# Chart 7: RMSE in each fold for both models.\n",
    "x = np.arange(5)\n",
    "plt.figure(figsize=(7, 4))\n",
    "plt.bar(x - 0.2, cv_lin, width=0.4, label=\"Linear Regression\", color=\"#9c9c9c\")\n",
    "plt.bar(x + 0.2, cv_rf, width=0.4, label=\"Random Forest\", color=\"#59a14f\")\n",
    "plt.xticks(x, [\"Fold 1\", \"Fold 2\", \"Fold 3\", \"Fold 4\", \"Fold 5\"])\n",
    "plt.ylabel(\"RMSE (t/ha) - lower is better\")\n",
    "plt.title(\"5-fold cross validation\")\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"chart7_kfold.png\")\n",
    "plt.close()\n",
    "\n",
    "# Chart 8: which features the forest relied on.\n",
    "importance = pd.Series(forest.feature_importances_, index=features).sort_values()\n",
    "print(\"Feature importance (Random Forest):\")\n",
    "print(importance.round(3).to_string())\n",
    "print()\n",
    "plt.figure(figsize=(7, 4.5))\n",
    "plt.barh(importance.index, importance.values, color=\"#f28e2b\")\n",
    "plt.title(\"What the Random Forest relied on\")\n",
    "plt.xlabel(\"Importance\")\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"chart8_importance.png\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 9 - Deployment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"--- Step 9: Deployment ---\")\n",
    "# The trained model is saved to a file, so an app or website can load it\n",
    "# and predict without training again.\n",
    "joblib.dump(forest, \"yield_model.pkl\")\n",
    "loaded = joblib.load(\"yield_model.pkl\")\n",
    "print(\"Model saved to yield_model.pkl and loaded back\")\n",
    "\n",
    "new_plots = pd.DataFrame({\n",
    "    \"rainfall_mm\": [1180, 560, 640], \"temperature_c\": [26.5, 18.5, 27.0], \"soil_ph\": [6.7, 7.1, 7.9],\n",
    "    \"nitrogen_kg_ha\": [85, 90, 60], \"fertiliser_kg_ha\": [130, 150, 40], \"irrigated\": [0, 1, 0],\n",
    "    \"crop_Maize\": [0, 0, 1], \"crop_Rice\": [1, 0, 0], \"crop_Wheat\": [0, 1, 0],\n",
    "})\n",
    "names = [\"Rice, rain-fed, good rain\", \"Wheat, irrigated\", \"Maize, rain-fed, little fertiliser\"]\n",
    "predicted = loaded.predict(new_plots)\n",
    "for i in range(len(names)):\n",
    "    print(names[i].ljust(36), \"->\", round(predicted[i], 2), \"t/ha\")\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 10 - Feedback"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"--- Step 10: Feedback ---\")\n",
    "# After the next season's harvest, the real yields come back. Comparing\n",
    "# them with what the model predicted shows whether it still works.\n",
    "nxt = pd.read_csv(\"next_season.csv\")\n",
    "nxt[\"irrigated\"] = (nxt[\"irrigated\"] == \"Yes\").astype(int)\n",
    "nxt = pd.get_dummies(nxt, columns=[\"crop\"], dtype=int)\n",
    "next_pred = loaded.predict(nxt[features])\n",
    "next_rmse = np.sqrt(mean_squared_error(nxt[\"yield_t_ha\"], next_pred))\n",
    "print(\"Plots in the next season's feedback :\", len(nxt))\n",
    "print(\"RMSE on last season's test data     :\", round(results[1][3], 3), \"t/ha\")\n",
    "print(\"RMSE on the next season             :\", round(next_rmse, 3), \"t/ha\")\n",
    "print(\"Average rainfall, training season   :\", round(data[\"rainfall_mm\"].mean()), \"mm\")\n",
    "print(\"Average rainfall, next season       :\", round(nxt[\"rainfall_mm\"].mean()), \"mm\")\n",
    "if next_rmse > 1.2 * results[1][3]:\n",
    "    print(\"Feedback: error has grown by more than 20% - retrain with the new season's data.\")\n",
    "else:\n",
    "    print(\"Feedback: error is within 20% of last season - the model can stay in use.\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}