trained model using sentinel 1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9027603a-f5f7-492b-b991-01e6b58b09d2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: tifffile in c:\\users\\rishi\\anaconda3\\lib\\site-packages (2023.4.12)Note: you may need to restart the kernel to use updated packages.\n",
      "\n",
      "Requirement already satisfied: numpy in c:\\users\\rishi\\anaconda3\\lib\\site-packages (from tifffile) (1.26.4)\n"
     ]
    }
   ],
   "source": [
    "pip install tifffile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "308b5e94-99f6-41d0-9048-b9d8f030d84d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shapefile CRS: EPSG:4326\n",
      "Number of features collected: 1408\n",
      "Number of targets collected: 1408\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'X_resampled' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[1], line 100\u001b[0m\n\u001b[0;32m     97\u001b[0m X_new \u001b[38;5;241m=\u001b[39m selector\u001b[38;5;241m.\u001b[39mfit_transform(X_scaled, y)\n\u001b[0;32m     99\u001b[0m \u001b[38;5;66;03m# Split the data\u001b[39;00m\n\u001b[1;32m--> 100\u001b[0m X_train, X_test, y_train, y_test \u001b[38;5;241m=\u001b[39m train_test_split(X_resampled, y_resampled, test_size\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0.2\u001b[39m, random_state\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m42\u001b[39m)\n\u001b[0;32m    102\u001b[0m \u001b[38;5;66;03m# Model tuning with GridSearchCV\u001b[39;00m\n\u001b[0;32m    103\u001b[0m param_grid \u001b[38;5;241m=\u001b[39m {\n\u001b[0;32m    104\u001b[0m     \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mn_estimators\u001b[39m\u001b[38;5;124m'\u001b[39m: [\u001b[38;5;241m100\u001b[39m, \u001b[38;5;241m200\u001b[39m, \u001b[38;5;241m300\u001b[39m],\n\u001b[0;32m    105\u001b[0m     \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mmax_depth\u001b[39m\u001b[38;5;124m'\u001b[39m: [\u001b[38;5;241m10\u001b[39m, \u001b[38;5;241m20\u001b[39m, \u001b[38;5;241m30\u001b[39m, \u001b[38;5;28;01mNone\u001b[39;00m],\n\u001b[0;32m    106\u001b[0m     \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mmax_features\u001b[39m\u001b[38;5;124m'\u001b[39m: [\u001b[38;5;124m'\u001b[39m\u001b[38;5;124msqrt\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlog2\u001b[39m\u001b[38;5;124m'\u001b[39m]  \u001b[38;5;66;03m# Avoid using 'auto'\u001b[39;00m\n\u001b[0;32m    107\u001b[0m }\n",
      "\u001b[1;31mNameError\u001b[0m: name 'X_resampled' is not defined"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import geopandas as gpd\n",
    "import rasterio\n",
    "import rasterio.mask\n",
    "from shapely.geometry import box\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import train_test_split, GridSearchCV\n",
    "from sklearn.metrics import mean_squared_error, accuracy_score\n",
    "from sklearn.preprocessing import StandardScaler, LabelEncoder\n",
    "from sklearn.feature_selection import SelectKBest, f_classif\n",
    "from skimage.transform import resize\n",
    "from joblib import Parallel, delayed\n",
    "from imblearn.over_sampling import SMOTE\n",
    "import matplotlib.pyplot as plt\n",
    "import joblib\n",
    "\n",
    "# Define the paths\n",
    "raster_directory = r\"C:\\Users\\rishi\\OneDrive\\Desktop\\sentinel_1\"\n",
    "shapefile_path = r\"C:\\Users\\rishi\\OneDrive\\Desktop\\merged for training\\crop_data_merged.shp\"\n",
    "output_directory = r\"C:\\Users\\rishi\\OneDrive\\Desktop\\sentinel_1\"\n",
    "\n",
    "# Create the output directory if it does not exist\n",
    "os.makedirs(output_directory, exist_ok=True)\n",
    "\n",
    "# Load the shapefile\n",
    "gdf = gpd.read_file(shapefile_path)\n",
    "\n",
    "# Print the CRS of the shapefile\n",
    "print(f\"Shapefile CRS: {gdf.crs}\")\n",
    "\n",
    "# Check if shapefile has the expected 'layer' column\n",
    "assert 'layer' in gdf.columns, \"Shapefile does not contain 'layer' column\"\n",
    "\n",
    "# Define the fixed shape for the features (e.g., 128x128)\n",
    "fixed_shape = (128, 128)\n",
    "\n",
    "def process_raster(filename):\n",
    "    raster_path = os.path.join(raster_directory, filename)\n",
    "    features = []\n",
    "    targets = []\n",
    "\n",
    "    with rasterio.open(raster_path) as src:\n",
    "        print(f\"Raster CRS: {src.crs}\")\n",
    "        if src.crs != gdf.crs:\n",
    "            print(f\"Reprojecting shapefile to match raster CRS\")\n",
    "            gdf_reprojected = gdf.to_crs(src.crs)\n",
    "        else:\n",
    "            gdf_reprojected = gdf\n",
    "\n",
    "        raster_bounds = box(*src.bounds)\n",
    "        print(f\"Raster bounds: {raster_bounds.bounds}\")\n",
    "\n",
    "        intersecting_features = gdf_reprojected[gdf_reprojected.intersects(raster_bounds)]\n",
    "        print(f\"Number of intersecting features: {len(intersecting_features)}\")\n",
    "\n",
    "        if not intersecting_features.empty:\n",
    "            for feature in intersecting_features.itertuples():\n",
    "                geom = feature.geometry\n",
    "                Layer = feature.layer\n",
    "                out_image, out_transform = rasterio.mask.mask(src, [geom], crop=True)\n",
    "                resized_image = resize(out_image[0], fixed_shape, mode='reflect', anti_aliasing=True)\n",
    "                features.append(resized_image.flatten())\n",
    "                targets.append(Layer)\n",
    "    \n",
    "    return features, targets\n",
    "\n",
    "# Use parallel processing to handle multiple files concurrently\n",
    "results = Parallel(n_jobs=-1)(delayed(process_raster)(filename) for filename in os.listdir(raster_directory) if filename.endswith('.tif'))\n",
    "\n",
    "# Combine results from all parallel processes\n",
    "all_features = []\n",
    "all_targets = []\n",
    "for f, t in results:\n",
    "    all_features.extend(f)\n",
    "    all_targets.extend(t)\n",
    "\n",
    "# Debugging: Print the number of features and targets collected\n",
    "print(f\"Number of features collected: {len(all_features)}\")\n",
    "print(f\"Number of targets collected: {len(all_targets)}\")\n",
    "\n",
    "assert len(all_features) == len(all_targets), \"Mismatch between number of images and target values\"\n",
    "\n",
    "X = np.array(all_features)\n",
    "label_encoder = LabelEncoder()\n",
    "y = label_encoder.fit_transform(np.array(all_targets))\n",
    "\n",
    "if len(X) == 0 or len(y) == 0:\n",
    "    raise ValueError(\"No samples collected. Ensure that the raster data and shapefile have corresponding overlaps.\")\n",
    "\n",
    "# Data preprocessing\n",
    "scaler = StandardScaler()\n",
    "X_scaled = scaler.fit_transform(X)\n",
    "\n",
    "# Feature selection\n",
    "selector = SelectKBest(score_func=f_classif, k='all')\n",
    "X_new = selector.fit_transform(X_scaled, y)\n",
    "\n",
    "# Handling imbalanced data\n",
    "smote = SMOTE(random_state=42)\n",
    "X_resampled, y_resampled = smote.fit_resample(X_new, y)\n",
    "\n",
    "# Split the data\n",
    "X_train, X_test, y_train, y_test = train_test_split(X_resampled, y_resampled, test_size=0.2, random_state=42)\n",
    "\n",
    "# Model tuning with GridSearchCV\n",
    "param_grid = {\n",
    "    'n_estimators': [100, 200, 300],\n",
    "    'max_depth': [10, 20, 30, None],\n",
    "    'max_features': ['sqrt', 'log2']  # Avoid using 'auto'\n",
    "}\n",
    "\n",
    "model = RandomForestClassifier(random_state=42, class_weight='balanced')\n",
    "grid_search = GridSearchCV(estimator=model, param_grid=param_grid, cv=5, scoring='accuracy')\n",
    "grid_search.fit(X_train, y_train)\n",
    "\n",
    "best_model = grid_search.best_estimator_\n",
    "print(f\"Best parameters found: {grid_search.best_params_}\")\n",
    "\n",
    "# Train the best model\n",
    "best_model.fit(X_train, y_train)\n",
    "y_pred = best_model.predict(X_test)\n",
    "\n",
    "# Calculate Mean Squared Error\n",
    "mse = mean_squared_error(y_test, y_pred)\n",
    "print(f\"Mean Squared Error: {mse}\")\n",
    "\n",
    "# Calculate Accuracy\n",
    "accuracy = accuracy_score(y_test, y_pred)\n",
    "print(f\"Accuracy: {accuracy}\")\n",
    "\n",
    "# Visualize the true vs predicted values\n",
    "plt.scatter(y_test, y_pred)\n",
    "plt.xlabel(\"True Values\")\n",
    "plt.ylabel(\"Predictions\")\n",
    "plt.title(\"True vs Predicted Values\")\n",
    "plt.show()\n",
    "\n",
    "label_mapping = dict(zip(label_encoder.classes_, label_encoder.transform(label_encoder.classes_)))\n",
    "print(\"Label Mapping:\", label_mapping)\n",
    "\n",
    "# Save the model\n",
    "model_path = r\"C:\\Users\\rishi\\OneDrive\\Desktop\\best_model.joblib\"\n",
    "joblib.dump(best_model, model_path)\n",
    "\n",
    "print(f\"Model saved to {model_path}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c0b82ee0-ea33-4339-aa7c-3eb9ca93cff1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shapefile CRS: EPSG:4326\n",
      "Number of features collected: 1408\n",
      "Number of targets collected: 1408\n",
      "Best parameters found: {'max_depth': 10, 'max_features': 'sqrt', 'n_estimators': 100}\n",
      "Mean Squared Error: 3.5\n",
      "Accuracy: 0.3191489361702128\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHFCAYAAAAOmtghAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA94ElEQVR4nO3de1iUdf7/8dco54TxCKKiWBmeOphmonnEMCnTn+2V20lMLbHUVdZV0bK+bUVtbtlJ7aCQ+TWtRctNK6kENXUTxSxBsxI1gzwllMrR+/eHX2abODgzDgzcPh/XdV9X85n3fc/7vrlrXs19z2cshmEYAgAAMIkGnm4AAADAnQg3AADAVAg3AADAVAg3AADAVAg3AADAVAg3AADAVAg3AADAVAg3AADAVAg3AADAVAg3wEWwWCwOLWlpaZ5utUaFh4fb7W+jRo104403aunSpbXy+snJybJYLMrJybGNDRgwQAMGDHB6W08//bTef/99t/VWLicnRxaLRcnJyVXWTJs2TRaLRXv37q2yZs6cObJYLNq5c6fDrx0eHq4xY8Y40S1QvxFugIuwdetWuyUmJkb+/v4Vxq+//npPt1rj+vTpY9vf8rARGxurhQsXeqSfBQsWaMGCBU6vV1PhxhHjxo2TJC1ZsqTS58+dO6elS5fquuuuuyTOKcBVXp5uAKjPevXqZfe4RYsWatCgQYXxPzpz5owCAgJqsrVa17hxY7v9Hjx4sNq1a6fnn39eEydOrHSdsrIylZaWytfX1+39dO7c2e3brGldu3ZVz5499fbbb+vpp5+Wl5f9f6LXr1+vH3/8UTNnzvRQh0D9wCc3QA0bMGCAunbtqo0bN6p3794KCAjQ2LFjJZ2/rPX4449XWKeyywh5eXmaMGGC2rRpIx8fH7Vv317/8z//o9LS0mpff8SIEWrXrp3OnTtX4bkbb7zR7hOA9957TzfeeKOsVqsCAgJ0+eWX23p1VuPGjRUREaGDBw9K+u9lmX/84x968skn1b59e/n6+mrDhg2SpIyMDN1+++1q2rSp/Pz81K1bN7377rsVtrtt2zb16dNHfn5+atWqlRISElRSUlKhrrLLUkVFRXriiSfUqVMn+fn5qVmzZho4cKC2bNki6fzf4/Tp03rrrbdsl9h+vw1H/wY//fST7rzzTgUGBspqtWrUqFHKy8tz6LiNGzdOeXl5+uijjyo8l5SUJF9fX91zzz0qLCzUX//6V1133XWyWq1q2rSpIiMj9cEHH1zwNSq7jCdJaWlplV5G/fTTTxUVFaWgoCAFBASoT58++uyzz+xqjh07pgcffFBhYWHy9fVVixYt1KdPH3366acO7TfgTnxyA9SC3Nxc3XvvvZoxY4aefvppNWjg3P9X5OXlqWfPnmrQoIHmzp2rK664Qlu3btWTTz6pnJwcJSUlVbnu2LFjNXz4cH3++ecaPHiwbXzv3r368ssv9dJLL0k6f4lt1KhRGjVqlB5//HH5+fnp4MGD+vzzz13a55KSEh08eFAtWrSwG3/ppZd01VVXad68eQoKClKHDh20YcMG3XLLLbrxxhu1aNEiWa1WrVixQqNGjdKZM2dsQS8rK0tRUVEKDw9XcnKyAgICtGDBAi1fvvyC/ZSWlmro0KHatGmTpk6dqkGDBqm0tFTbtm3ToUOH1Lt3b23dulWDBg3SwIED9eijj0qSgoKCJDn+Nzh79qwGDx6sn376SYmJibrqqqu0du1ajRo1yqHjdtddd2natGlasmSJhg0bZhv/5Zdf9MEHH+j//b//pyZNmig/P18nT57U9OnT1bp1axUXF+vTTz/VyJEjlZSUpNGjRzv0eheybNkyjR49WsOHD9dbb70lb29vvfbaaxoyZIg++eQTRUVFSZLuu+8+7dy5U0899ZSuuuoqnTp1Sjt37tSJEyfc0gfgFAOA28TGxhqXXXaZ3Vj//v0NScZnn31WoV6S8dhjj1UYb9eunREbG2t7PGHCBKNRo0bGwYMH7ermzZtnSDL27NlTZU8lJSVGSEiIcffdd9uNz5gxw/Dx8TGOHz9ut61Tp05daDcr7TcmJsYoKSkxSkpKjAMHDhixsbGGJONvf/ubYRiGceDAAUOSccUVVxjFxcV263fs2NHo1q2bUVJSYjd+2223GaGhoUZZWZlhGIYxatQow9/f38jLy7PVlJaWGh07djQkGQcOHLCN9+/f3+jfv7/t8dKlSw1JxhtvvFHtvlx22WV2x76co3+DhQsXGpKMDz74wK7ugQceMCQZSUlJ1b6+YZw/j7y9vY2ff/7ZNvbyyy8bkozU1NRK1yktLTVKSkqMcePGGd26dbN77o/nU1JSUoXjZRiGsWHDBkOSsWHDBsMwDOP06dNG06ZNjWHDhtnVlZWVGddee63Rs2dP21ijRo2MqVOnXnDfgNrAZSmgFjRp0kSDBg1yef0PP/xQAwcOVKtWrVRaWmpbhg4dKklKT0+vcl0vLy/de++9WrVqlfLz8yWdv9fl7bff1vDhw9WsWTNJ0g033CBJuvPOO/Xuu+/qyJEjTvW4bt06eXt7y9vbW+3bt9e7776ryZMn68knn7Sru/322+Xt7W17/N1332nv3r265557JMlu/2JiYpSbm6t9+/ZJkjZs2KCoqCiFhITY1m/YsKFDn4p89NFH8vPzc/kym6N/gw0bNigwMFC333673fp33323w681btw4lZSU6O2337aNJSUlqV27drZPSqTzlxH79OmjRo0aycvLS97e3lq8eLGys7Nd2sc/2rJli06ePKnY2Fi7fT537pxuueUWbd++XadPn5Yk9ezZU8nJyXryySe1bdu2Si8VArWFcAPUgtDQ0Ita/+eff9a///1vW3goX7p06SJJOn78eLXrjx07VoWFhVqxYoUk6ZNPPlFubq7uv/9+W02/fv30/vvvq7S0VKNHj1abNm3UtWtXvfPOOw71eNNNN2n79u3KyMhQVlaWTp06pZdeekk+Pj52dX88Fj///LMkafr06RX276GHHrLbvxMnTqhly5YVXruysT86duyYWrVq5fQlwd/36cjf4MSJE3bhy5key/Xt21dXXXWV7VLX7t27tXPnTt1///2yWCySpFWrVunOO+9U69attWzZMm3dulXbt2+3/a3dofxv86c//anCfj/77LMyDEMnT56UJK1cuVKxsbF68803FRkZqaZNm2r06NEO32sEuBP33AC1oPwN6Y98fX1VVFRUYfyP9yk0b95c11xzjZ566qlKt9OqVatqX79z587q2bOnkpKSNGHCBCUlJalVq1aKjo62qxs+fLiGDx+uoqIibdu2TYmJibr77rsVHh6uyMjIal/DarWqR48e1dZIFY9F8+bNJUkJCQkaOXJkpetERERIkpo1a1bpm6Ujb6AtWrTQ5s2bde7cOZcCjqN/g2bNmunLL790qcffGzt2rGbNmqUvv/xSy5cvV4MGDexuMl+2bJnat2+vlStX2h3Tys6nP/Lz86u09o8hufxv8/LLL1f5DcDyINe8eXPNnz9f8+fP16FDh7RmzRrNmjVLR48e1ccff3zhHQbciHADeFB4eLh2795tN/b555/rt99+sxu77bbbtG7dOl1xxRVq0qSJS691//33a+LEidq8ebP+/e9/Kz4+Xg0bNqy01tfXV/3791fjxo31ySefKDMz84LhxlURERHq0KGDvvrqKz399NPV1g4cOFBr1qzRzz//bHtTLSsr08qVKy/4OkOHDtU777yj5OTkai9N+fr66uzZsxXGHf0bDBw4UO+++67WrFljd2nKkZuefy82NlaPPPKIXnvtNa1Zs0ZRUVFq166d7XmLxSIfHx+7YJOXl+fQt6XCw8Mlnf9EqDw4StKaNWvs6vr06aPGjRsrKytLkyZNcrj3tm3batKkSfrss8/0xRdfOLwe4C6EG8CD7rvvPj366KOaO3eu+vfvr6ysLL3yyiuyWq12dU888YRSU1PVu3dvTZkyRRERESosLFROTo7WrVunRYsWqU2bNtW+1l133aX4+HjdddddKioqqvBV87lz5+rHH39UVFSU2rRpo1OnTunFF1+Ut7e3+vfv7+5dt/Paa69p6NChGjJkiMaMGaPWrVvr5MmTys7O1s6dO/Xee+9Jkh555BGtWbNGgwYN0ty5cxUQEKBXX33Vdt9Hde666y4lJSUpLi5O+/bt08CBA3Xu3Dn95z//UadOnfTnP/9ZknT11VcrLS1N//73vxUaGqrAwEBFREQ4/DcYPXq0XnjhBY0ePVpPPfWUOnTooHXr1umTTz5x6pi0bNlSMTExSkpKkmEYtgn+yt12221atWqVHnroIf3pT3/S4cOH9fe//12hoaHav39/tdu+4YYbFBERoenTp6u0tFRNmjTR6tWrtXnzZru6Ro0a6eWXX1ZsbKxOnjypP/3pTwoODtaxY8f01Vdf6dixY1q4cKHy8/M1cOBA3X333erYsaMCAwO1fft2ffzxx1V+GgfUKE/f0QyYSVXflurSpUul9UVFRcaMGTOMsLAww9/f3+jfv7+xa9euCt9uMQzDOHbsmDFlyhSjffv2hre3t9G0aVOje/fuxpw5c4zffvvNof7uvvtuQ5LRp0+fCs99+OGHxtChQ43WrVsbPj4+RnBwsBETE2Ns2rTpgttt166dceutt1ZbU/5tqeeee67S57/66ivjzjvvNIKDgw1vb2+jZcuWxqBBg4xFixbZ1X3xxRdGr169DF9fX6Nly5bG3/72N+P111+/4LelDMMwzp49a8ydO9fo0KGD4ePjYzRr1swYNGiQsWXLFlvNrl27jD59+hgBAQGGJLttOPo3+PHHH4077rjDaNSokREYGGjccccdxpYtWxz+tlS5Dz74wJBkNG3a1CgsLKzw/DPPPGOEh4cbvr6+RqdOnYw33njDeOyxx4w//qe9svPp22+/NaKjo42goCCjRYsWxuTJk421a9fafVuqXHp6unHrrbcaTZs2Nby9vY3WrVsbt956q/Hee+8ZhmEYhYWFRlxcnHHNNdcYQUFBhr+/vxEREWE89thjxunTpx3eX8BdLIZhGJ6LVgAAAO7Ft6UAAICpEG4AAICpEG4AAICpEG4AAICpEG4AAICpEG4AAICpXHKT+J07d04//fSTAgMDq5wSHwAA1C2GYejXX3916DfiLrlw89NPPyksLMzTbQAAABccPnz4gjOyX3LhJjAwUNL5gxMUFOThbgAAgCMKCgoUFhZmex+vziUXbsovRQUFBRFuAACoZxy5pYQbigEAgKkQbgAAgKkQbgAAgKkQbgAAgKkQbgAAgKkQbgAAgKkQbgAAgKkQbgAAgKkQbgAAgKlccjMU15TfCks1bWWmDv1yVm2b+OuFUd3UyI/Di4uTd6pQt728UQWFpQry89KHk/upZWM/T7cFEzhbXKan12Up58QZhTcL0OyYzvL3aejptlDPLVj/jf7x+UHb4xmD2umh6K613ofFMAyj1l/1/yxcuFALFy5UTk6OJKlLly6aO3euhg4dWuU66enpio+P1549e9SqVSvNmDFDcXFxDr9mQUGBrFar8vPz3fbzC7e/skm7fyyoMH5NmyCtmdTXLa+BS0+nRz/S2ZJzFcb9vRso++9V/zsCXMgDS7crNetohfGbOwfrjdE3eKAjmEH4rLVVPpfzzK0XvX1n3r89elmqTZs2euaZZ5SRkaGMjAwNGjRIw4cP1549eyqtP3DggGJiYtS3b19lZmZq9uzZmjJlilJSUmq58/+qKthI0u4fC3T7K5tquSOYQVXBRpLOlpxTp0c/quWOYBZVBRtJSs06qgeWbq/ljmAG1QUbR553N4+Gm2HDhikmJkZXXXWVrrrqKj311FNq1KiRtm3bVmn9okWL1LZtW82fP1+dOnXS+PHjNXbsWM2bN6+WOz/vt8LSKoNNud0/Fui3wtJa6ghmkHeqsMpgU+5syTnlnSqspY5gFmeLy6oMNuVSs47qbHFZLXUEM1iw/hu31rlDnbmhuKysTCtWrNDp06cVGRlZac3WrVsVHR1tNzZkyBBlZGSopKSk0nWKiopUUFBgt7jLtJWZbq0DJOm2lze6tQ4o9/S6LLfWAZLs7rFxR507eDzcfP3112rUqJF8fX0VFxen1atXq3PnzpXW5uXlKSQkxG4sJCREpaWlOn78eKXrJCYmymq12pawsDC39X7ol7NurQMkqcDBT/ocrQPK5Zw449Y6oK7yeLiJiIjQrl27tG3bNk2cOFGxsbHKyqr6/xosFovd4/L7of84Xi4hIUH5+fm25fDhw27rvW0Tf7fWAZIU5OC37BytA8qFNwtwax1QV3k83Pj4+OjKK69Ujx49lJiYqGuvvVYvvvhipbUtW7ZUXl6e3djRo0fl5eWlZs2aVbqOr6+vgoKC7BZ3eWFUN7fWAZL04eR+bq0Dys2OqfxTcVfrAOn8173dWecOHg83f2QYhoqKiip9LjIyUqmpqXZj69evV48ePeTt7V0b7dlp5Oela9pUH5auaRPEfDdwSsvGfvL3rv5fTX/vBsx3A6f5+zTUzZ2Dq625uXMw893AKY7OY1Ob8914NNzMnj1bmzZtUk5Ojr7++mvNmTNHaWlpuueeeySdv6Q0evRoW31cXJwOHjyo+Ph4ZWdna8mSJVq8eLGmT5/uqV3Qmkl9qww4zHMDV2X/fWiVAYd5bnAx3hh9Q5UBh3lu4KoLzWPjjnlunOHRSfzGjRunzz77TLm5ubJarbrmmms0c+ZM3XzzzZKkMWPGKCcnR2lpabZ10tPTNW3aNNskfjNnzvT4JH4SMxSjZjBDMWoKMxSjJtTkDMXOvH97NNx4Qk2FGwAAUHPqzQzFAAAA7ka4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApuLRcJOYmKgbbrhBgYGBCg4O1ogRI7Rv375q10lLS5PFYqmw7N27t5a6BgAAdZlHw016eroefvhhbdu2TampqSotLVV0dLROnz59wXX37dun3Nxc29KhQ4da6BgAANR1Xp588Y8//tjucVJSkoKDg7Vjxw7169ev2nWDg4PVuHHjGuwOAADUR3Xqnpv8/HxJUtOmTS9Y261bN4WGhioqKkobNmyosq6oqEgFBQV2CwAAMK86E24Mw1B8fLxuuukmde3atcq60NBQvf7660pJSdGqVasUERGhqKgobdy4sdL6xMREWa1W2xIWFlZTuwAAAOoAi2EYhqebkKSHH35Ya9eu1ebNm9WmTRun1h02bJgsFovWrFlT4bmioiIVFRXZHhcUFCgsLEz5+fkKCgq66L4BAEDNKygokNVqdej9u058cjN58mStWbNGGzZscDrYSFKvXr20f//+Sp/z9fVVUFCQ3QIAAMzLozcUG4ahyZMna/Xq1UpLS1P79u1d2k5mZqZCQ0Pd3B0AAKiPPBpuHn74YS1fvlwffPCBAgMDlZeXJ0myWq3y9/eXJCUkJOjIkSNaunSpJGn+/PkKDw9Xly5dVFxcrGXLliklJUUpKSke2w8AAFB3eDTcLFy4UJI0YMAAu/GkpCSNGTNGkpSbm6tDhw7ZnisuLtb06dN15MgR+fv7q0uXLlq7dq1iYmJqq20AAFCH1ZkbimuLMzckAQCAuqHe3VAMAADgLoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKoQbAABgKl6ebgBA1crOGfrywEkd/bVQwYF+6tm+qRo2sHi6LZjA2eIyPb0uSzknzii8WYBmx3SWv09DT7cFuIVHw01iYqJWrVqlvXv3yt/fX71799azzz6riIiIatdLT09XfHy89uzZo1atWmnGjBmKi4urpa6B2vHxN7n6n39nKTe/0DYWavXTY8M665auoR7sDPXdA0u3KzXrqO3xpv3S29sO6ebOwXpj9A0e7AxwD49elkpPT9fDDz+sbdu2KTU1VaWlpYqOjtbp06erXOfAgQOKiYlR3759lZmZqdmzZ2vKlClKSUmpxc6BmvXxN7mauGynXbCRpLz8Qk1ctlMff5Proc5Q3/0x2PxeatZRPbB0ey13BLifxTAMw9NNlDt27JiCg4OVnp6ufv36VVozc+ZMrVmzRtnZ2baxuLg4ffXVV9q6desFX6OgoEBWq1X5+fkKCgpyW++Au5SdM3TTs59XCDblLJJaWv20eeYgLlHBKWeLy9Rp7scXrMt+4hYuUaHOceb9u07dUJyfny9Jatq0aZU1W7duVXR0tN3YkCFDlJGRoZKSkgr1RUVFKigosFuAuuzLAyerDDaSZEjKzS/UlwdO1l5TMIWn12W5tQ6oq+pMuDEMQ/Hx8brpppvUtWvXKuvy8vIUEhJiNxYSEqLS0lIdP368Qn1iYqKsVqttCQsLc3vvgDsd/bXqYONKHVAu58QZt9YBdVWdCTeTJk3S7t279c4771yw1mKx/yi+/MraH8clKSEhQfn5+bbl8OHD7mkYqCHBgX5urQPKhTcLcGsdUFfViXAzefJkrVmzRhs2bFCbNm2qrW3ZsqXy8vLsxo4ePSovLy81a9asQr2vr6+CgoLsFqAu69m+qUKtfqrqbhqLzn9rqmf7qi/fApWZHdPZrXVAXeXRcGMYhiZNmqRVq1bp888/V/v27S+4TmRkpFJTU+3G1q9frx49esjb27umWgVqTcMGFj027Pybyx8DTvnjx4Z15mZiOM3fp6Fu7hxcbc3NnYO5mRj1nkfDzcMPP6xly5Zp+fLlCgwMVF5envLy8nT27FlbTUJCgkaPHm17HBcXp4MHDyo+Pl7Z2dlasmSJFi9erOnTp3tiF4AacUvXUC2893q1tNpfempp9dPCe69nnhu47I3RN1QZcJjnBmbh0a+CV3aPjCQlJSVpzJgxkqQxY8YoJydHaWlptufT09M1bdo02yR+M2fOdHgSP74KjvqEGYpRU5ihGPWNM+/fdWqem9pAuAEAoP6pt/PcAAAAXCzCDQAAMBXCDQAAMBXCDQAAMBXCDQAAMBXCDQAAMBXCDQAAMBXCDQAAMBXCDQAAMBXCDQAAMBXCDQAAMBXCDQAAMBXCDQAAMBXCDQAAMBWXws3OnTv19ddf2x5/8MEHGjFihGbPnq3i4mK3NQcAAOAsl8LNhAkT9O2330qSfvjhB/35z39WQECA3nvvPc2YMcOtDQIAADjDpXDz7bff6rrrrpMkvffee+rXr5+WL1+u5ORkpaSkuLM/AAAAp7gUbgzD0Llz5yRJn376qWJiYiRJYWFhOn78uPu6AwAAcJJL4aZHjx568skn9fbbbys9PV233nqrJOnAgQMKCQlxa4MAAADOcCnczJ8/Xzt37tSkSZM0Z84cXXnllZKkf/3rX+rdu7dbGwQAAHCGxTAMw10bKywsVMOGDeXt7e2uTbpdQUGBrFar8vPzFRQU5Ol2AACAA5x5//a6mBcqLi7W0aNHbffflGvbtu3FbBYAAMBlLoWbb7/9VuPGjdOWLVvsxg3DkMViUVlZmVuaAwAAcJZL4eb++++Xl5eXPvzwQ4WGhspisbi7LwAAAJe4FG527dqlHTt2qGPHju7uBwAA4KK49G2pzp07M58NAACok1wKN88++6xmzJihtLQ0nThxQgUFBXYLAACAp7j0VfAGDc5noj/ea1Mfbijmq+AAANQ/Nf5V8A0bNrjUGAAAQE1zKdz079/f3X0AAAC4hcuT+J06dUqLFy9Wdna2LBaLOnfurLFjx8pqtbqzPwAAAKe4dENxRkaGrrjiCr3wwgs6efKkjh8/rueff15XXHGFdu7c6e4eAQAAHObSDcV9+/bVlVdeqTfeeENeXuc//CktLdX48eP1ww8/aOPGjW5v1F24oRgAgPrHmfdvl8KNv7+/MjMzK0zil5WVpR49eujMmTPObrLWEG4AAKh/nHn/dumyVFBQkA4dOlRh/PDhwwoMDHRlkwAAAG7hUrgZNWqUxo0bp5UrV+rw4cP68ccftWLFCo0fP1533XWXu3sEAABwmEvflpo3b54sFotGjx6t0tJSSZK3t7cmTpyoZ555xq0NAgAAOMOle27KnTlzRt9//70Mw9CVV16pgIAAd/ZWI7jnBgCA+qfGZyguFxAQoKuvvvpiNgEAAOBWDoebkSNHKjk5WUFBQRo5cmS1tatWrbroxgAAAFzhcLixWq22H8oMCgqq8KOZAAAAdcFF3XNTH3HPDQAA9U+Nz3MzaNAgnTp1qtIXHjRokCubBAAAcAuXwk1aWpqKi4srjBcWFmrTpk0X3RQAAICrnPq21O7du23/nJWVpby8PNvjsrIyffzxx2rdurX7ugMAAHCSU+Hmuuuuk8VikcViqfTyk7+/v15++WWHt7dx40Y999xz2rFjh3Jzc7V69WqNGDGiyvq0tDQNHDiwwnh2dnaF37kCAACXJqfCzYEDB2QYhi6//HJ9+eWXatGihe05Hx8fBQcHq2HDhg5v7/Tp07r22mt1//3364477nB4vX379tndTPT7PgAAwKXNqXDTrl07SdK5c+fc8uJDhw7V0KFDnV4vODhYjRs3dksPAADAXFy6oTgxMVFLliypML5kyRI9++yzF93UhXTr1k2hoaGKiorShg0bqq0tKipSQUGB3QIAAMzLpXDz2muvVXqPS5cuXbRo0aKLbqoqoaGhev3115WSkqJVq1YpIiJCUVFR2rhxY5XrJCYmymq12pawsLAa6w8AAHieS5P4+fn5KTs7W+3bt7cb/+GHH9S5c2cVFhY634jFcsEbiiszbNgwWSwWrVmzptLni4qKVFRUZHtcUFCgsLAwJvEDAKAeqfFJ/MLCwvTFF19UGP/iiy/UqlUrVzbpsl69emn//v1VPu/r66ugoCC7BQAAmJdLvwo+fvx4TZ06VSUlJbavhH/22WeaMWOG/vrXv7q1wQvJzMxUaGhorb4mAACou1wKNzNmzNDJkyf10EMP2WYq9vPz08yZM5WQkODwdn777Td99913tscHDhzQrl271LRpU7Vt21YJCQk6cuSIli5dKkmaP3++wsPD1aVLFxUXF2vZsmVKSUlRSkqKK7sBAABMyKVwY7FY9Oyzz+rRRx9Vdna2/P391aFDB/n6+jq1nYyMDLtJ+eLj4yVJsbGxSk5OVm5urg4dOmR7vri4WNOnT9eRI0fk7++vLl26aO3atYqJiXFlNwAAgAnxq+AAAKDOc+b92+FPbkaOHKnk5GQFBQVp5MiR1dauWrXK0c0CAAC4lcPhxmq1ymKx2P4ZAACgLuKyFAAAqPNqfJ4bAACAusrhy1LdunWzXZa6kJ07d7rcEAAAwMVwONz8/mcRCgsLtWDBAnXu3FmRkZGSpG3btmnPnj166KGH3N4kAACAoxwON4899pjtn8ePH68pU6bo73//e4Waw4cPu687AAAAJ7l0Q7HValVGRoY6dOhgN75//3716NFD+fn5bmvQ3bihGACA+qfGbyj29/fX5s2bK4xv3rxZfn5+rmwSAADALVz6+YWpU6dq4sSJ2rFjh3r16iXp/D03S5Ys0dy5c93aIAAAgDNcCjezZs3S5ZdfrhdffFHLly+XJHXq1EnJycm688473dogAACAM5jEDwAA1Hm1MonfqVOn9Oabb2r27Nk6efKkpPPz2xw5csTVTQIAAFw0ly5L7d69W4MHD5bValVOTo7Gjx+vpk2bavXq1Tp48KCWLl3q7j4BAAAc4tInN/Hx8RozZoz2799v9+2ooUOHauPGjW5rDgAAwFkuhZvt27drwoQJFcZbt26tvLy8i24KAADAVS6FGz8/PxUUFFQY37dvn1q0aHHRTQEAALjKpXAzfPhwPfHEEyopKZEkWSwWHTp0SLNmzdIdd9zh1gYBAACc4VK4mTdvno4dO6bg4GCdPXtW/fv315VXXqnAwEA99dRT7u4RAADAYS59WyooKEibN2/W559/rp07d+rcuXO6/vrrNXjwYHf3BwAA4BSnw01paan8/Py0a9cuDRo0SIMGDaqJvgAAAFzi9GUpLy8vtWvXTmVlZTXRDwAAwEVx6Z6bRx55RAkJCbaZiQEAAOoKl+65eemll/Tdd9+pVatWateunS677DK753fu3OmW5gAAAJzlUrgZMWKELBaLLrHf3AQAAPWAU+HmzJkz+tvf/qb3339fJSUlioqK0ssvv6zmzZvXVH8AAABOceqem8cee0zJycm69dZbddddd+nTTz/VxIkTa6o3AAAApzn1yc2qVau0ePFi/fnPf5Yk3XPPPerTp4/KysrUsGHDGmkQAADAGU59cnP48GH17dvX9rhnz57y8vLSTz/95PbGAAAAXOFUuCkrK5OPj4/dmJeXl0pLS93aFAAAgKucuixlGIbGjBkjX19f21hhYaHi4uLsvg6+atUq93UIAADgBKfCTWxsbIWxe++9123NAAAAXCynwk1SUlJN9QEAAOAWLv38AgAAQF1FuAEAAKZCuAEAAKZCuAEAAKZCuAEAAKZCuAEAAKZCuAEAAKZCuAEAAKZCuAEAAKZCuAEAAKZCuAEAAKbi1G9LoWonfyvWn1/foqO/Fis40EcrHuytpo18PN0W6rl1X/6oh1Z9ZXu8YOS1iunZxoMdAUDVzhaX6el1Wco5cUbhzQI0O6az/H0a1nofFsMwjFp/1f+zceNGPffcc9qxY4dyc3O1evVqjRgxotp10tPTFR8frz179qhVq1aaMWOG4uLiHH7NgoICWa1W5efnKygo6CL34LwbnkzVsd+KK4y3aOSj7Y/c7JbXwKUnfNbaKp/LeebWWuwEAC7sgaXblZp1tML4zZ2D9cboGy56+868f3v0stTp06d17bXX6pVXXnGo/sCBA4qJiVHfvn2VmZmp2bNna8qUKUpJSanhTqtWVbCRpGO/FeuGJ1NruSOYQXXBxpHnAaA2VRVsJCk166geWLq9Vvvx6GWpoUOHaujQoQ7XL1q0SG3bttX8+fMlSZ06dVJGRobmzZunO+64o4a6rNrJ34qrDDbljv1WrJO/FXOJCg5b9+WPDtdxiQqAp50tLqsy2JRLzTqqs8VltXaJql7dULx161ZFR0fbjQ0ZMkQZGRkqKSmpdJ2ioiIVFBTYLe7y59e3uLUOkGR3j4076gCgJj29Lsutde5Qr8JNXl6eQkJC7MZCQkJUWlqq48ePV7pOYmKirFarbQkLC3NbP0d/rf5TG2frAACob3JOnHFrnTvUq3AjSRaLxe5x+f3Qfxwvl5CQoPz8fNty+PBht/USHOjYpSZH6wAAqG/CmwW4tc4d6lW4admypfLy8uzGjh49Ki8vLzVr1qzSdXx9fRUUFGS3uMuKB3u7tQ6Qzn/d2511AFCTZsd0dmudO9SrcBMZGanUVPtvH61fv149evSQt7d3rffTtJGPWlzgRuEWjXy4mRhOcfQmYW4mBlAX+Ps01M2dg6utublzcK3Od+PRcPPbb79p165d2rVrl6TzX/XetWuXDh06JOn8JaXRo0fb6uPi4nTw4EHFx8crOztbS5Ys0eLFizV9+nRPtC9J2v7IzVUGHOa5gasuNI8N89wAqEveGH1DlQHHXfPcOMOjk/ilpaVp4MCBFcZjY2OVnJysMWPGKCcnR2lpabbn0tPTNW3aNNskfjNnzvT4JH4SMxSjZjBDMYD6pCZnKHbm/duj4cYTaircAACAmlNvZigGAABwN8INAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFcINAAAwFY+HmwULFqh9+/by8/NT9+7dtWnTpipr09LSZLFYKix79+6txY4BAEBd5tFws3LlSk2dOlVz5sxRZmam+vbtq6FDh+rQoUPVrrdv3z7l5ubalg4dOtRSxwAAoK7zaLh5/vnnNW7cOI0fP16dOnXS/PnzFRYWpoULF1a7XnBwsFq2bGlbGjZsWEsdAwCAus5j4aa4uFg7duxQdHS03Xh0dLS2bNlS7brdunVTaGiooqKitGHDhmpri4qKVFBQYLcAAADz8li4OX78uMrKyhQSEmI3HhISory8vErXCQ0N1euvv66UlBStWrVKERERioqK0saNG6t8ncTERFmtVtsSFhbm1v0AAAB1i5enG7BYLHaPDcOoMFYuIiJCERERtseRkZE6fPiw5s2bp379+lW6TkJCguLj422PCwoKCDgAAJiYxz65ad68uRo2bFjhU5qjR49W+DSnOr169dL+/furfN7X11dBQUF2CwAAMC+PhRsfHx91795dqampduOpqanq3bu3w9vJzMxUaGiou9sDAAD1lEcvS8XHx+u+++5Tjx49FBkZqddff12HDh1SXFycpPOXlI4cOaKlS5dKkubPn6/w8HB16dJFxcXFWrZsmVJSUpSSkuLJ3QAAAHWIR8PNqFGjdOLECT3xxBPKzc1V165dtW7dOrVr106SlJubazfnTXFxsaZPn64jR47I399fXbp00dq1axUTE+OpXQAAAHWMxTAMw9NN1KaCggJZrVbl5+dz/w0AAPWEM+/fHv/5BQAAAHci3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFPx8nQDZvH1oXzdvmCzDEkWSWseuklXt7V6ui3Uc1Pe/kxr9hTaHt/exU8v3RflwY5gFlk/Fui2VzbpnM7/X+6Hk/qqc5sgT7eFeu6e+Wv1Rd5/H/dpKf3v1FtrvQ+Pf3KzYMECtW/fXn5+furevbs2bdpUbX16erq6d+8uPz8/XX755Vq0aFEtdVq18FlrNez/go0kGZKGLdis8FlrPdkW6rnwWWvtgo0krdlTyHmFixY+a61i/i/YSNI5STGvbOLcwkUJn2UfbCTpizx55LzyaLhZuXKlpk6dqjlz5igzM1N9+/bV0KFDdejQoUrrDxw4oJiYGPXt21eZmZmaPXu2pkyZopSUlFru/L8u9EfjPxZwBecVagrnFmpCXTuvPBpunn/+eY0bN07jx49Xp06dNH/+fIWFhWnhwoWV1i9atEht27bV/Pnz1alTJ40fP15jx47VvHnzarnz874+lO/WOkA6fynKnXVAuawfC9xaB0jnL0W5s84dPBZuiouLtWPHDkVHR9uNR0dHa8uWLZWus3Xr1gr1Q4YMUUZGhkpKSipdp6ioSAUFBXaLu9y+YLNb6wBJFS5FXWwdUO62V6q/7O9sHSCpwqWoi61zB4+Fm+PHj6usrEwhISF24yEhIcrLq/wI5OXlVVpfWlqq48ePV7pOYmKirFarbQkLC3PPDki2e2zcVQcANenchUucqgPqKo/fUGyxWOweG4ZRYexC9ZWNl0tISFB+fr5tOXz48EV2/Lte3FwHADXJ0f/ge/yNAbhIHjuHmzdvroYNG1b4lObo0aMVPp0p17Jly0rrvby81KxZs0rX8fX1VVBQkN3iLmseusmtdYB0/uve7qwDyn04qa9b6wDp/Ne93VnnDh4LNz4+PurevbtSU1PtxlNTU9W7d+9K14mMjKxQv379evXo0UPe3t411mtVHJ3Hhvlu4AxH57Fhvhs4y9F5bJjvBs5wdB6b2pzvxqOfPsbHx+vNN9/UkiVLlJ2drWnTpunQoUOKi4uTdP6S0ujRo231cXFxOnjwoOLj45Wdna0lS5Zo8eLFmj59uqd2QTnPVP/HutDzQGU4r1BTOLdQE+raeWUxym9a8ZAFCxboH//4h3Jzc9W1a1e98MIL6tevnyRpzJgxysnJUVpamq0+PT1d06ZN0549e9SqVSvNnDnTFoYcUVBQIKvVqvz8fLdeomKGYtQEZihGTWGGYtSEmpyh2Jn3b4+Hm9pWU+EGAADUHGfev7kpHgAAmArhBgAAmArhBgAAmArhBgAAmArhBgAAmArhBgAAmArhBgAAmArhBgAAmArhBgAAmIqXpxuobeUTMhcUFHi4EwAA4Kjy921Hfljhkgs3v/76qyQpLCzMw50AAABn/frrr7Jaq//txkvut6XOnTunn376SYGBgbJYLG7ddkFBgcLCwnT48GF+t+oCOFaO41g5jmPlHI6X4zhWjqupY2UYhn799Ve1atVKDRpUf1fNJffJTYMGDdSmTZsafY2goCBOfgdxrBzHsXIcx8o5HC/HcawcVxPH6kKf2JTjhmIAAGAqhBsAAGAqhBs38vX11WOPPSZfX19Pt1Lncawcx7FyHMfKORwvx3GsHFcXjtUld0MxAAAwNz65AQAApkK4AQAApkK4AQAApkK4AQAApkK4cdKCBQvUvn17+fn5qXv37tq0aVO19enp6erevbv8/Px0+eWXa9GiRbXUqec5c6zS0tJksVgqLHv37q3Fjj1j48aNGjZsmFq1aiWLxaL333//gutcqueVs8fqUj2vEhMTdcMNNygwMFDBwcEaMWKE9u3bd8H1LtXzypXjdameWwsXLtQ111xjm6AvMjJSH330UbXreOK8Itw4YeXKlZo6darmzJmjzMxM9e3bV0OHDtWhQ4cqrT9w4IBiYmLUt29fZWZmavbs2ZoyZYpSUlJqufPa5+yxKrdv3z7l5ubalg4dOtRSx55z+vRpXXvttXrllVccqr+Uzytnj1W5S+28Sk9P18MPP6xt27YpNTVVpaWlio6O1unTp6tc51I+r1w5XuUutXOrTZs2euaZZ5SRkaGMjAwNGjRIw4cP1549eyqt99h5ZcBhPXv2NOLi4uzGOnbsaMyaNavS+hkzZhgdO3a0G5swYYLRq1evGuuxrnD2WG3YsMGQZPzyyy+10F3dJclYvXp1tTWX8nn1e44cK86r844ePWpIMtLT06us4bz6L0eOF+fWfzVp0sR48803K33OU+cVn9w4qLi4WDt27FB0dLTdeHR0tLZs2VLpOlu3bq1QP2TIEGVkZKikpKTGevU0V45VuW7duik0NFRRUVHasGFDTbZZb12q59XFuNTPq/z8fElS06ZNq6zhvPovR45XuUv53CorK9OKFSt0+vRpRUZGVlrjqfOKcOOg48ePq6ysTCEhIXbjISEhysvLq3SdvLy8SutLS0t1/PjxGuvV01w5VqGhoXr99deVkpKiVatWKSIiQlFRUdq4cWNttFyvXKrnlSs4r87/knJ8fLxuuukmde3atco6zqvzHD1el/K59fXXX6tRo0by9fVVXFycVq9erc6dO1da66nz6pL7VfCLZbFY7B4bhlFh7EL1lY2bkTPHKiIiQhEREbbHkZGROnz4sObNm6d+/frVaJ/10aV8XjmD80qaNGmSdu/erc2bN1+wlvPK8eN1KZ9bERER2rVrl06dOqWUlBTFxsYqPT29yoDjifOKT24c1Lx5czVs2LDCJw9Hjx6tkErLtWzZstJ6Ly8vNWvWrMZ69TRXjlVlevXqpf3797u7vXrvUj2v3OVSOq8mT56sNWvWaMOGDWrTpk21tZxXzh2vylwq55aPj4+uvPJK9ejRQ4mJibr22mv14osvVlrrqfOKcOMgHx8fde/eXampqXbjqamp6t27d6XrREZGVqhfv369evToIW9v7xrr1dNcOVaVyczMVGhoqLvbq/cu1fPKXS6F88owDE2aNEmrVq3S559/rvbt219wnUv5vHLleFXmUji3KmMYhoqKiip9zmPnVY3ermwyK1asMLy9vY3FixcbWVlZxtSpU43LLrvMyMnJMQzDMGbNmmXcd999tvoffvjBCAgIMKZNm2ZkZWUZixcvNry9vY1//etfntqFWuPssXrhhReM1atXG99++63xzTffGLNmzTIkGSkpKZ7ahVrz66+/GpmZmUZmZqYhyXj++eeNzMxM4+DBg4ZhcF79nrPH6lI9ryZOnGhYrVYjLS3NyM3NtS1nzpyx1XBe/Zcrx+tSPbcSEhKMjRs3GgcOHDB2795tzJ4922jQoIGxfv16wzDqznlFuHHSq6++arRr187w8fExrr/+eruvCsbGxhr9+/e3q09LSzO6detm+Pj4GOHh4cbChQtruWPPceZYPfvss8YVV1xh+Pn5GU2aNDFuuukmY+3atR7ouvaVf6X0j0tsbKxhGJxXv+fssbpUz6vKjpEkIykpyVbDefVfrhyvS/XcGjt2rO2/6y1atDCioqJswcYw6s55ZTGM/7uzBwAAwAS45wYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QbAJW/AgAGaOnWqp9sA4CaEGwBOs1gs1S5jxoyplT6GDRumwYMHV/rc1q1bZbFYtHPnzlrpBUDd4eXpBgDUP7m5ubZ/XrlypebOnat9+/bZxvz9/e3qS0pKauRH8saNG6eRI0fq4MGDateund1zS5Ys0XXXXafrr7/e7a8LoG7jkxsATmvZsqVtsVqtslgstseFhYVq3Lix3n33XQ0YMEB+fn5atmyZHn/8cV133XV225k/f77Cw8PtxpKSktSpUyf5+fmpY8eOWrBgQZV93HbbbQoODlZycrLd+JkzZ7Ry5UqNGzdOJ06c0F133aU2bdooICBAV199td55551q989isej999+3G2vcuLHd6xw5ckSjRo1SkyZN1KxZMw0fPlw5OTm259PS0tSzZ09ddtllaty4sfr06aODBw9W+7oA3INwA6BGzJw5U1OmTFF2draGDBni0DpvvPGG5syZo6eeekrZ2dl6+umn9eijj+qtt96qtN7Ly0ujR49WcnKyfv8zee+9956Ki4t1zz33qLCwUN27d9eHH36ob775Rg8++KDuu+8+/ec//3F5386cOaOBAweqUaNG2rhxozZv3qxGjRrplltuUXFxsUpLSzVixAj1799fu3fv1tatW/Xggw/KYrG4/JoAHMdlKQA1YurUqRo5cqRT6/z973/XP//5T9t67du3V1ZWll577TXFxsZWus7YsWP13HPPKS0tTQMHDpR0/pLUyJEj1aRJEzVp0kTTp0+31U+ePFkff/yx3nvvPd14440u7duKFSvUoEEDvfnmm7bAkpSUpMaNGystLU09evRQfn6+brvtNl1xxRWSpE6dOrn0WgCcR7gBUCN69OjhVP2xY8d0+PBhjRs3Tg888IBtvLS0VFartcr1OnbsqN69e2vJkiUaOHCgvv/+e23atEnr16+XJJWVlemZZ57RypUrdeTIERUVFamoqEiXXXaZazsmaceOHfruu+8UGBhoN15YWKjvv/9e0dHRGjNmjIYMGaKbb75ZgwcP1p133qnQ0FCXXxOA4wg3AGrEH8NDgwYN7C4dSedvNC537tw5SecvTf3xE5WGDRtW+1rjxo3TpEmT9OqrryopKUnt2rVTVFSUJOmf//ynXnjhBc2fP19XX321LrvsMk2dOlXFxcVVbs9isVyw1+7du+t///d/K6zbokULSec/yZkyZYo+/vhjrVy5Uo888ohSU1PVq1evavcFwMUj3ACoFS1atFBeXp4Mw7Bdytm1a5ft+ZCQELVu3Vo//PCD7rnnHqe2feedd+ovf/mLli9frrfeeksPPPCA7TU2bdqk4cOH695775V0Ppjs37+/2stELVq0sPtG2P79+3XmzBnb4+uvv14rV65UcHCwgoKCqtxOt27d1K1bNyUkJCgyMlLLly8n3AC1gBuKAdSKAQMG6NixY/rHP/6h77//Xq+++qo++ugju5rHH39ciYmJevHFF/Xtt9/q66+/VlJSkp5//vlqt92oUSONGjVKs2fP1k8//WQ3z86VV16p1NRUbdmyRdnZ2ZowYYLy8vKq3d6gQYP0yiuvaOfOncrIyFBcXJzdV9nvueceNW/eXMOHD9emTZt04MABpaen6y9/+Yt+/PFHHThwQAkJCdq6dasOHjyo9evX69tvv+W+G6CWEG4A1IpOnTppwYIFevXVV3Xttdfqyy+/tLvRV5LGjx+vN998U8nJybr66qvVv39/JScnq3379hfc/rhx4/TLL79o8ODBatu2rW380Ucf1fXXX68hQ4ZowIABatmypUaMGFHttv75z38qLCxM/fr10913363p06crICDA9nxAQIA2btyotm3bauTIkerUqZPGjh2rs2fPKigoSAEBAdq7d6/uuOMOXXXVVXrwwQc1adIkTZgwwbmDBsAlFuOPF5YBAADqMT65AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApvL/AbUFso3LPf3QAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Label Mapping: {'Crop_class_Garlic': 0, 'Crop_class_Onion1': 1, 'Crop_class_Wheat1': 2, 'Crop_class_Wheat2': 3}\n",
      "Model saved to C:\\Users\\rishi\\OneDrive\\Desktop\\best_model.joblib\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import geopandas as gpd\n",
    "import rasterio\n",
    "import rasterio.mask\n",
    "from shapely.geometry import box\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import train_test_split, GridSearchCV\n",
    "from sklearn.metrics import mean_squared_error, accuracy_score\n",
    "from sklearn.preprocessing import StandardScaler, LabelEncoder\n",
    "from sklearn.feature_selection import SelectKBest, f_classif\n",
    "from skimage.transform import resize\n",
    "from joblib import Parallel, delayed\n",
    "import matplotlib.pyplot as plt\n",
    "import joblib\n",
    "\n",
    "# Define the paths\n",
    "raster_directory = r\"C:\\Users\\rishi\\OneDrive\\Desktop\\sentinel_1\"\n",
    "shapefile_path = r\"C:\\Users\\rishi\\OneDrive\\Desktop\\merged for training\\crop_data_merged.shp\"\n",
    "output_directory = r\"C:\\Users\\rishi\\OneDrive\\Desktop\\sentinel_1\"\n",
    "\n",
    "# Create the output directory if it does not exist\n",
    "os.makedirs(output_directory, exist_ok=True)\n",
    "\n",
    "# Load the shapefile\n",
    "gdf = gpd.read_file(shapefile_path)\n",
    "\n",
    "# Print the CRS of the shapefile\n",
    "print(f\"Shapefile CRS: {gdf.crs}\")\n",
    "\n",
    "# Check if shapefile has the expected 'layer' column\n",
    "assert 'layer' in gdf.columns, \"Shapefile does not contain 'layer' column\"\n",
    "\n",
    "# Define the fixed shape for the features (e.g., 128x128)\n",
    "fixed_shape = (128, 128)\n",
    "\n",
    "def process_raster(filename):\n",
    "    raster_path = os.path.join(raster_directory, filename)\n",
    "    features = []\n",
    "    targets = []\n",
    "\n",
    "    with rasterio.open(raster_path) as src:\n",
    "        print(f\"Raster CRS: {src.crs}\")\n",
    "        if src.crs != gdf.crs:\n",
    "            print(f\"Reprojecting shapefile to match raster CRS\")\n",
    "            gdf_reprojected = gdf.to_crs(src.crs)\n",
    "        else:\n",
    "            gdf_reprojected = gdf\n",
    "\n",
    "        raster_bounds = box(*src.bounds)\n",
    "        print(f\"Raster bounds: {raster_bounds.bounds}\")\n",
    "\n",
    "        intersecting_features = gdf_reprojected[gdf_reprojected.intersects(raster_bounds)]\n",
    "        print(f\"Number of intersecting features: {len(intersecting_features)}\")\n",
    "\n",
    "        if not intersecting_features.empty:\n",
    "            for feature in intersecting_features.itertuples():\n",
    "                geom = feature.geometry\n",
    "                Layer = feature.layer\n",
    "                out_image, out_transform = rasterio.mask.mask(src, [geom], crop=True)\n",
    "                resized_image = resize(out_image[0], fixed_shape, mode='reflect', anti_aliasing=True)\n",
    "                features.append(resized_image.flatten())\n",
    "                targets.append(Layer)\n",
    "    \n",
    "    return features, targets\n",
    "\n",
    "# Use parallel processing to handle multiple files concurrently\n",
    "results = Parallel(n_jobs=-1)(delayed(process_raster)(filename) for filename in os.listdir(raster_directory) if filename.endswith('.tif'))\n",
    "\n",
    "# Combine results from all parallel processes\n",
    "all_features = []\n",
    "all_targets = []\n",
    "for f, t in results:\n",
    "    all_features.extend(f)\n",
    "    all_targets.extend(t)\n",
    "\n",
    "# Debugging: Print the number of features and targets collected\n",
    "print(f\"Number of features collected: {len(all_features)}\")\n",
    "print(f\"Number of targets collected: {len(all_targets)}\")\n",
    "\n",
    "assert len(all_features) == len(all_targets), \"Mismatch between number of images and target values\"\n",
    "\n",
    "X = np.array(all_features)\n",
    "label_encoder = LabelEncoder()\n",
    "y = label_encoder.fit_transform(np.array(all_targets))\n",
    "\n",
    "if len(X) == 0 or len(y) == 0:\n",
    "    raise ValueError(\"No samples collected. Ensure that the raster data and shapefile have corresponding overlaps.\")\n",
    "\n",
    "# Data preprocessing\n",
    "scaler = StandardScaler()\n",
    "X_scaled = scaler.fit_transform(X)\n",
    "\n",
    "# Feature selection\n",
    "selector = SelectKBest(score_func=f_classif, k='all')\n",
    "X_new = selector.fit_transform(X_scaled, y)\n",
    "\n",
    "# Split the data\n",
    "X_train, X_test, y_train, y_test = train_test_split(X_new, y, test_size=0.2, random_state=42)\n",
    "\n",
    "# Model tuning with GridSearchCV\n",
    "param_grid = {\n",
    "    'n_estimators': [100, 200, 300],\n",
    "    'max_depth': [10, 20, 30, None],\n",
    "    'max_features': ['sqrt', 'log2']  # Avoid using 'auto'\n",
    "}\n",
    "\n",
    "model = RandomForestClassifier(random_state=42, class_weight='balanced')\n",
    "grid_search = GridSearchCV(estimator=model, param_grid=param_grid, cv=5, scoring='accuracy')\n",
    "grid_search.fit(X_train, y_train)\n",
    "\n",
    "best_model = grid_search.best_estimator_\n",
    "print(f\"Best parameters found: {grid_search.best_params_}\")\n",
    "\n",
    "# Train the best model\n",
    "best_model.fit(X_train, y_train)\n",
    "y_pred = best_model.predict(X_test)\n",
    "\n",
    "# Calculate Mean Squared Error\n",
    "mse = mean_squared_error(y_test, y_pred)\n",
    "print(f\"Mean Squared Error: {mse}\")\n",
    "\n",
    "# Calculate Accuracy\n",
    "accuracy = accuracy_score(y_test, y_pred)\n",
    "print(f\"Accuracy: {accuracy}\")\n",
    "\n",
    "# Visualize the true vs predicted values\n",
    "plt.scatter(y_test, y_pred)\n",
    "plt.xlabel(\"True Values\")\n",
    "plt.ylabel(\"Predictions\")\n",
    "plt.title(\"True vs Predicted Values\")\n",
    "plt.show()\n",
    "\n",
    "label_mapping = dict(zip(label_encoder.classes_, label_encoder.transform(label_encoder.classes_)))\n",
    "print(\"Label Mapping:\", label_mapping)\n",
    "\n",
    "# Save the model\n",
    "model_path = r\"C:\\Users\\rishi\\OneDrive\\Desktop\\best_model.joblib\"\n",
    "joblib.dump(best_model, model_path)\n",
    "\n",
    "print(f\"Model saved to {model_path}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7b93b08-c956-4510-8fb6-40a9c2306b42",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}