inerpolating_points_1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import netCDF4\n",
    "#import thread"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#this is the training data for our models.\n",
    "filepath='/home/vedanshi/Documents/DSML/IMD work/HARv2_d10km_d_2d_prcp_nc_2023.nc'\n",
    "data=xr.open_dataset(filepath,engine='netcdf4')\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.lon.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.prcp_nc.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.prcp_nc[102].T.plot.surface()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.prcp_nc[0, :, :].T.plot.surface()\n",
    "data.prcp_nc[0].max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "visualizing latitude and longitude as south_north and west_east"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(data.lon.values))\n",
    "print(len(data.lat.values))\n",
    "df1=pd.DataFrame(data.lon)\n",
    "df1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#checking for duplicates \n",
    "print(sum(df1.duplicated()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df2=pd.DataFrame(data.lat)\n",
    "df2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#checking for duplicated\n",
    "print(sum(df2.duplicated()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df2.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "extracting latitude and longitude values for day one, i.e these values represent the data recorded on respective lon and lat on 1st day, \n",
    "#### Plan :\n",
    "- make the method just for 1 day and then itrate over for the next 365 days.\n",
    "- then interpolate for that day."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_lat = data.lat[:,0].values\n",
    "latitude = np.array(new_lat)\n",
    "# Flatten the 2D array to a 1D array\n",
    "#flattened_lat = new_lat.flatten()\n",
    "# Iterate over the flattened array\n",
    "#for lat in flattened_lat:\n",
    " #   latitude.append(lat)\n",
    "df3=pd.DataFrame(latitude)\n",
    "df3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_lon = data.lon[0,:].values\n",
    "longitude = np.array(new_lon)\n",
    "# Flatten the 2D array to a 1D array\n",
    "#flattened_lon = new_lon.flatten()\n",
    "# Iterate over the flattened array\n",
    "#for lon in flattened_lon:\n",
    " #   longitude.append(lon)\n",
    "df4=pd.DataFrame(longitude)\n",
    "df4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "now i have to extract the precipitaion values in mm somehow form the prcp_nc variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "# Extracting precipitation data directly as a NumPy array\n",
    "prcp_values = data['prcp_nc'][0, :, :].values\n",
    "# Flatten the 2D array to a 1D array\n",
    "prcp_values_flat = prcp_values.flatten()\n",
    "# Now, prcp_values_flat contains all precipitation values as a 1D array\n",
    "print(prcp_values_flat.max()==data.prcp_nc[0].max())\n",
    "print((len(latitude)*len(longitude)==len(prcp_values_flat)))\n",
    "print(len(prcp_values_flat))\n",
    "print(len(latitude)*len(longitude))\n",
    "prcp_values_flat.max()\n",
    "# We have now gotten all the prcp values, now we have to fill them in a dataframe, but the dimenssions dont match "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = prcp_values = data['prcp_nc'].shape\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Making a function of the above process so as it can be used when we need to extract the infromation of the remaining days."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# i is the index at which all the data is to be exctracted, we will the above data only.\n",
    "def extraction(i):\n",
    "    lati=data.lat[:,i]\n",
    "    lat=np.array(lati)\n",
    "    longi=data.lon[i,:]\n",
    "    lon=np.array(longi)\n",
    "    prcp=data['prcp_nc'][i,:,:].values\n",
    "    prcp_falt=prcp.flatten()\n",
    "    return prcp_falt,lon,lat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Maping of precipitaion vlaues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "indexr=[]\n",
    "for i in range(252):\n",
    "    indexr.append(i)\n",
    "indexr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a list of column labels based on the length of latitude\n",
    "columns_labels = list(range(len(longitude)))\n",
    "# Create an empty DataFrame with the correct number of columns\n",
    "temp_df = pd.DataFrame(index=indexr, columns=columns_labels)\n",
    "print(len(temp_df))\n",
    "temp_df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#populating the data frame.\n",
    "for i in range(381):\n",
    "    for j in range(252):\n",
    "        temp_df[i][j]=prcp_values_flat[i*j]\n",
    "temp_df.index=latitude\n",
    "temp_df.columns=longitude\n",
    "temp_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Define latitude and longitude ranges\n",
    "lat_range = np.array(range(0, 91, 10))  # Positive latitude range (0 to 90)\n",
    "lon_range = np.array(range(0, 181, 10))  # Positive longitude range (0 to 180)\n",
    "\n",
    "# Plot the grid without any data\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.grid(True)  # Plot grid lines\n",
    "plt.xticks(lon_range)  # Set x-axis ticks to longitude values\n",
    "plt.yticks(lat_range)  # Set y-axis ticks to latitude values\n",
    "plt.xlabel('Longitude')\n",
    "plt.ylabel('Latitude')\n",
    "plt.title('Grid without Data (Positive Latitude and Longitude)')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a 2D grid of points for the HAR data\n",
    "har_lat, har_lon = np.meshgrid(latitude, longitude)\n",
    "len(har_lon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inter=pd.read_excel(r'/home/vedanshi/Documents/DSML/IMD work/avalanche_glacier sites.xlsx')\n",
    "inter.columns=[0,1,2,3,4,5,6,7,8,9,10]\n",
    "inter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "to_be_found_lat=inter[1:][0]\n",
    "to_be_found_lon=inter[1:][1]\n",
    "to_be_found_lat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "to_be_found_lon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# making grid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "x_values and y_values are the lat_renge and lon_range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_grid=np.empty((252,381),dtype=object)\n",
    "# 252 - lat values.\n",
    "# 381 - lon values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to calculate grid parameters\n",
    "def calculate_grid_parameters(coord):\n",
    "    d = np.diff(coord)[0] / 2\n",
    "    print(d)\n",
    "    cell_boundaries = np.concatenate([[coord[0] - d], coord + d])\n",
    "    return cell_boundaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to create grid cell taking arguments from x_values and y_values\n",
    "def initialize_and_populate_grid(x_values, y_values): #x_values and y_value = latitude, longitude.\n",
    "    ngrid_cell = len(x_values) * len(y_values)\n",
    "    grid_cell = [{} for i in range(ngrid_cell)]\n",
    "    dy = y_values[1] - y_values[0]\n",
    "    k = 0\n",
    "    \n",
    "    latcell = calculate_grid_parameters(x_values)\n",
    "    loncell = calculate_grid_parameters(y_values)\n",
    "    \n",
    "    for i, lat1  in enumerate(x_values):\n",
    "        for j, lon1 in enumerate(y_values):\n",
    "            x_coords = [round(latcell[i], 2), round(latcell[i], 2), round(latcell[i+1], 2), round(latcell[i+1], 2), round(latcell[i], 2)]\n",
    "            y_coords = [round(loncell[j], 2), round(loncell[j+1], 2), round(loncell[j+1], 2), round(loncell[j], 2), round(loncell[j], 2)] if dy > 0 else \\\n",
    "                       [round(loncell[j+1], 2), round(loncell[j], 2), round(loncell[j], 2), round(loncell[j+1], 2), round(loncell[j+1], 2)]\n",
    "            \n",
    "            grid_cell[k] = {'X': x_coords, 'Y': y_coords, 'Xc': round(lat1, 2), 'Yc': round(lon1, 2)}\n",
    "            k += 1\n",
    "    return grid_cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if __name__ == \"__main__\":\n",
    "    grid_cell = initialize_and_populate_grid(latitude, longitude)\n",
    "grid_cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_inter=[]\n",
    "for i in range(1,21):\n",
    "    interp=[to_be_found_lat[i],to_be_found_lon[i]]\n",
    "    temp_inter.append(interp)\n",
    "to_be_inter=pd.DataFrame(temp_inter)\n",
    "to_be_inter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Flatten the arrays to use in griddata\n",
    "prcp_points = np.column_stack((har_lat.ravel(), har_lon.ravel()))\n",
    "prcp=prcp_values_flat.ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Shape of prcp_points:\", prcp_points.shape)\n",
    "print(\"Shape of prcp_values_flat:\", prcp_values_flat.shape)\n",
    "print(\"Shape of to_be_inter:\", to_be_inter.shape)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ineterpolating values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roundoff_lat=latitude.round(decimals=1)\n",
    "roundoff_lon=longitude.round(decimals=1)\n",
    "roundoff_lat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roundoff_lon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "the values rounded off above could be used to extrac the vlaues of y_true."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def accuracy(y_true,y_interpolated):\n",
    "    # Calculate MAE\n",
    "    mae = np.mean(np.abs(y_interpolated - y_true))\n",
    "    # Calculate MSE\n",
    "    mse = np.mean((y_interpolated - y_true)**2)\n",
    "    # Calculate RMSE\n",
    "    rmse = np.sqrt(mse)\n",
    "    # Calculate R^2 score\n",
    "    mean_true = np.mean(y_true)\n",
    "    r2 = 1 - np.sum((y_true - y_interpolated)**2) / np.sum((y_true - mean_true)**2)\n",
    "    print(\"MAE:\", mae)\n",
    "    print(\"MSE:\", mse)\n",
    "    print(\"RMSE:\", rmse)\n",
    "    print(\"R^2 Score:\", r2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def interpolate(coords,method_used):\n",
    "    #entered coords must be upto accurcay of only 1 digit after the decimal.\n",
    "    import scipy.interpolate as sci\n",
    "    if method_used=='linear' or method_used=='bilinear':\n",
    "        if method_used=='bilinear':\n",
    "            return sci.griddata(prcp_points,prcp_values_flat,coords,method='nearest') #bilinear\n",
    "        else:\n",
    "            return sci.griddata(prcp_points,prcp_values_flat,coords,method='linear') #linear\n",
    "    else :\n",
    "            inter=sci.NearestNDInterpolator(prcp_points,prcp_values_flat,rescale=False,tree_options=None) #nearest neighbour\n",
    "            return inter(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interpolate(to_be_inter,'bilinear')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interpolate(to_be_inter,'linear')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interpolate(to_be_inter,'knn')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}