script/regressionCompareL7L8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Band to band regression\n",
    "This notebook reads the batch exported geotiff file of paried Landsat 7 and Landsat 8 images and convert the pixel values to vaex dataframe.\n",
    "Both Ordinary Least Square (OLS) Regression Model and Reduced Major Axis Regression (RMA) Model.\n",
    "RMA was calculated using the python package from https://github.com/OceanOptics/pylr2.\n",
    "\n",
    "Users will need to change the input and output folder. "
   ]
  },
  {
   "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.metrics import mean_squared_error #, r2_score\n",
    "from scipy import stats\n",
    "import seaborn as sns\n",
    "import vaex as vx\n",
    "from pylr2.regress2 import regress2\n",
    "from geopyfsn import getBand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_theme(style=\"darkgrid\", font=\"Arial\")\n",
    "SMALL_SIZE = 12\n",
    "MEDIUM_SIZE = 14\n",
    "BIGGER_SIZE = 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "blue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "folderPath = r\"/data/shunan/data/harmonize_data/201305_08landsat/blue\"\n",
    "L8, L7 = getBand(folderPath)\n",
    "df = pd.DataFrame({'L8': L8, 'L7': L7})\n",
    "df = vx.from_pandas(df)\n",
    "del L7,L8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# plot\n",
    "slope, intercept, r_value, p_value, std_err = stats.linregress(df.L7.values, df.L8.values)\n",
    "slope1, intercept1, r_value1, p_value1, std_err1 = stats.linregress(df.L8.values, df.L7.values)\n",
    "rma_results = regress2(df.L7.values, df.L8.values, _method_type_2=\"reduced major axis\")\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,7))\n",
    "plt.rc('font', size=BIGGER_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=BIGGER_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=BIGGER_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=BIGGER_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=BIGGER_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=BIGGER_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "\n",
    "# plot regression model\n",
    "plt.plot(np.array([0,1]), slope * np.array([0,1]) + intercept, color='red') # ols regression etm+ vs oli\n",
    "plt.plot([0,1], [0,1], color = 'white') # reference line\n",
    "# plt.plot(np.array([0,1]), (np.array([0,1]) - intercept1) / slope1 + intercept, color='peru') # ols regression oli vs etm+ \n",
    "plt.xlim(0, 1)\n",
    "plt.ylim(0, 1)\n",
    "# plt.axis('equal')\n",
    "k = rma_results['slope']\n",
    "b = rma_results['intercept'] \n",
    "plt.plot(np.array([0,1]), k * np.array([0,1]) + b, color='black') # rma regression\n",
    "\n",
    "ax.annotate('n:%.0f' % (len(df.L8.values)), xy=(0.7, 0.1),  xycoords='data',\n",
    "            horizontalalignment='left', verticalalignment='top',\n",
    "            )\n",
    "# scatter plot  \n",
    "df.viz.heatmap(df.L7, df.L8, what=np.log(vx.stat.count()), show=True, colormap=\"viridis\", vmin=0.0, vmax=6,\n",
    "    xlabel='Blue L7 Surface Reflectance', ylabel='Blue L8 Surface Reflectance')\n",
    "ax.set_aspect('equal', 'box')\n",
    "fig.savefig('print/L7L8/Blue.jpg', dpi=300)\n",
    "\n",
    "print('OLS: \\ny={0:.4f}x+{1:.4f}\\nOLS_r:{2:.2f}, p:{3:.3f}, n:{4:.0f}'.format(slope,intercept,r_value,p_value, len(df.L8.values)))\n",
    "print('RMA: \\ny={0:.4f}x+{1:.4f}\\nRMA_r:{2:.2f}'.format(k,b, rma_results['r']))\n",
    "print('RMSE is %.4f' % (mean_squared_error(df.L8.values, df.L7.values, squared=False)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sns.residplot(x=df.L7.values, y=df.L8.values)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.viz.histogram(x=df.L7, label='L7')   \n",
    "df.viz.histogram(x=df.L8, label='L8')   \n",
    "# df.viz.histogram(x=df.L7 * 1.1017 - 0.0084, label='L7adjusted')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "green"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "folderPath = r\"/data/shunan/data/harmonize_data/201305_08landsat/green\"\n",
    "# folderPath = r\"/data/shunan/data/harmonize_data/201305_08landsat/green\"\n",
    "L8, L7 = getBand(folderPath)\n",
    "df = pd.DataFrame({'L8': L8, 'L7': L7})\n",
    "df = vx.from_pandas(df)\n",
    "del L7,L8\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot\n",
    "slope, intercept, r_value, p_value, std_err = stats.linregress(df.L7.values, df.L8.values)\n",
    "slope1, intercept1, r_value1, p_value1, std_err1 = stats.linregress(df.L8.values, df.L7.values)\n",
    "# slope, intercept, r_value, p_value, std_err = stats.linregress(np.log(L7),np.log(L8))\n",
    "\n",
    "# Compute regression type 2\n",
    "rma_results = regress2(df.L7.values, df.L8.values, _method_type_2=\"reduced major axis\")\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,7))\n",
    "plt.rc('font', size=BIGGER_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=BIGGER_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=BIGGER_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=BIGGER_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=BIGGER_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=BIGGER_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "\n",
    "# plot regression model\n",
    "plt.plot(np.array([0,1]), slope * np.array([0,1]) + intercept, color='red') # ols regression etm+ vs oli\n",
    "plt.plot([0,1], [0,1], color = 'white') # reference line\n",
    "# plt.plot(np.array([0,1]), (np.array([0,1]) - intercept1) / slope1 + intercept, color='peru') # ols regression oli vs etm+ \n",
    "plt.xlim(0, 1)\n",
    "plt.ylim(0, 1)\n",
    "# plt.axis('equal')\n",
    "k = rma_results['slope']\n",
    "b = rma_results['intercept'] \n",
    "plt.plot(np.array([0,1]), k * np.array([0,1]) + b, color='black') # rma regression\n",
    "\n",
    "ax.annotate('n:%.0f' % (len(df.L8.values)), xy=(0.7, 0.1),  xycoords='data',\n",
    "            horizontalalignment='left', verticalalignment='top',\n",
    "            )\n",
    "\n",
    "# scatter plot\n",
    "df.viz.heatmap(df.L7, df.L8, what=np.log(vx.stat.count()), show=True, colormap=\"viridis\", vmin=0.0, vmax=6,\n",
    "    xlabel='Green L7 Surface Reflectance', ylabel='Green L8 Surface Reflectance')\n",
    "ax.set_aspect('equal', 'box')    \n",
    "fig.savefig('print/L7L8/Green.jpg', dpi=300)   \n",
    "\n",
    "print('OLS: \\ny={0:.4f}x+{1:.4f}\\nOLS_r:{2:.2f}, p:{3:.3f}, n:{4:.0f}'.format(slope,intercept,r_value,p_value, len(df.L8.values)))\n",
    "print('RMA: \\ny={0:.4f}x+{1:.4f}\\nRMA_r:{2:.2f}'.format(k,b, rma_results['r']))\n",
    "print('RMSE is %.4f' % (mean_squared_error(df.L8.values, df.L7.values, squared=False)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sns.residplot(x=df.L7.values, y=df.L8.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.viz.histogram(x=df.L7, label='L7')   \n",
    "df.viz.histogram(x=df.L8, label='L8')   \n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "red"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "folderPath = r\"/data/shunan/data/harmonize_data/201305_08landsat/red\"\n",
    "L8, L7 = getBand(folderPath)\n",
    "df = pd.DataFrame({'L8': L8, 'L7': L7})\n",
    "df = vx.from_pandas(df)\n",
    "del L7,L8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot\n",
    "slope, intercept, r_value, p_value, std_err = stats.linregress(df.L7.values, df.L8.values)\n",
    "slope1, intercept1, r_value1, p_value1, std_err1 = stats.linregress(df.L8.values, df.L7.values)\n",
    "# slope, intercept, r_value, p_value, std_err = stats.linregress(np.log(L7),np.log(L8))\n",
    "\n",
    "# Compute regression type 2\n",
    "rma_results = regress2(df.L7.values, df.L8.values, _method_type_2=\"reduced major axis\")\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,7))\n",
    "plt.rc('font', size=BIGGER_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=BIGGER_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=BIGGER_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=BIGGER_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=BIGGER_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=BIGGER_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "\n",
    "# plot regression model\n",
    "plt.plot(np.array([0,1]), slope * np.array([0,1]) + intercept, color='red') # ols regression etm+ vs oli\n",
    "plt.plot([0,1], [0,1], color = 'white') # reference line\n",
    "# plt.plot(np.array([0,1]), (np.array([0,1]) - intercept1) / slope1 + intercept, color='peru') # ols regression oli vs etm+ \n",
    "plt.xlim(0, 1)\n",
    "plt.ylim(0, 1)\n",
    "# plt.axis('equal')\n",
    "k = rma_results['slope']\n",
    "b = rma_results['intercept'] \n",
    "plt.plot(np.array([0,1]), k * np.array([0,1]) + b, color='black') # rma regression\n",
    "\n",
    "ax.annotate('n:%.0f' % (len(df.L8.values)), xy=(0.7, 0.1),  xycoords='data',\n",
    "            horizontalalignment='left', verticalalignment='top',\n",
    "            )\n",
    "\n",
    "# scatter plot\n",
    "df.viz.heatmap(df.L7, df.L8, what=np.log(vx.stat.count()), show=True, colormap=\"viridis\", vmin=0.0, vmax=6,\n",
    "    xlabel='Red L7 Surface Reflectance', ylabel='Red L8 Surface Reflectance')\n",
    "ax.set_aspect('equal', 'box')    \n",
    "fig.savefig('print/L7L8/Red.jpg', dpi=300) \n",
    "\n",
    "print('OLS: \\ny={0:.4f}x+{1:.4f}\\nOLS_r:{2:.2f}, p:{3:.3f}, n:{4:.0f}'.format(slope,intercept,r_value,p_value, len(df.L8.values)))\n",
    "print('RMA: \\ny={0:.4f}x+{1:.4f}\\nRMA_r:{2:.2f}'.format(k,b, rma_results['r']))\n",
    "print('RMSE is %.4f' % (mean_squared_error(df.L8.values, df.L7.values, squared=False)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sns.residplot(x=df.L7.values, y=df.L8.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.viz.histogram(x=df.L7, label='L7')   \n",
    "df.viz.histogram(x=df.L8, label='L8')   \n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "nir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "folderPath = r\"/data/shunan/data/harmonize_data/201305_08landsat/nir\"\n",
    "L8, L7 = getBand(folderPath)\n",
    "df = pd.DataFrame({'L8': L8, 'L7': L7})\n",
    "df = vx.from_pandas(df)\n",
    "del L7,L8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot\n",
    "slope, intercept, r_value, p_value, std_err = stats.linregress(df.L7.values, df.L8.values)\n",
    "slope1, intercept1, r_value1, p_value1, std_err1 = stats.linregress(df.L8.values, df.L7.values)\n",
    "# slope, intercept, r_value, p_value, std_err = stats.linregress(np.log(L7),np.log(L8))\n",
    "\n",
    "# Compute regression type 2\n",
    "rma_results = regress2(df.L7.values, df.L8.values, _method_type_2=\"reduced major axis\")\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,7))\n",
    "plt.rc('font', size=BIGGER_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=BIGGER_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=BIGGER_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=BIGGER_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=BIGGER_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=BIGGER_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "\n",
    "# plot regression model\n",
    "plt.plot(np.array([0,1]), slope * np.array([0,1]) + intercept, color='red') # ols regression etm+ vs oli\n",
    "plt.plot([0,1], [0,1], color = 'white') # reference line\n",
    "# plt.plot(np.array([0,1]), (np.array([0,1]) - intercept1) / slope1 + intercept, color='peru') # ols regression oli vs etm+ \n",
    "plt.xlim(0, 1)\n",
    "plt.ylim(0, 1)\n",
    "# plt.axis('equal')\n",
    "k = rma_results['slope']\n",
    "b = rma_results['intercept'] \n",
    "plt.plot(np.array([0,1]), k * np.array([0,1]) + b, color='black') # rma regression\n",
    "\n",
    "ax.annotate('n:%.0f' % (len(df.L8.values)), xy=(0.7, 0.1),  xycoords='data',\n",
    "            horizontalalignment='left', verticalalignment='top',\n",
    "            )\n",
    "\n",
    "# scatter plot\n",
    "df.viz.heatmap(df.L7, df.L8, what=np.log(vx.stat.count()), show=True, colormap=\"viridis\", vmin=0.0, vmax=6,\n",
    "    xlabel='NIR L7 Surface Reflectance', ylabel='NIR L8 Surface Reflectance')\n",
    "ax.set_aspect('equal', 'box')\n",
    "fig.savefig(r'print/L7L8/Nir.jpg', dpi=300)   \n",
    "\n",
    "print('OLS: \\ny={0:.4f}x+{1:.4f}\\nOLS_r:{2:.2f}, p:{3:.3f}, n:{4:.0f}'.format(slope,intercept,r_value,p_value, len(df.L8.values)))\n",
    "print('RMA: \\ny={0:.4f}x+{1:.4f}\\nRMA_r:{2:.2f}'.format(k,b, rma_results['r']))\n",
    "print('RMSE is %.4f' % (mean_squared_error(df.L8.values, df.L7.values, squared=False)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sns.residplot(x=df.L7.values, y=df.L8.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.viz.histogram(x=df.L7, label='L7')   \n",
    "df.viz.histogram(x=df.L8, label='L8')   \n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "swir1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "folderPath = r\"/data/shunan/data/harmonize_data/201305_08landsat/swir1\"\n",
    "L8, L7 = getBand(folderPath)\n",
    "df = pd.DataFrame({'L8': L8, 'L7': L7})\n",
    "df = vx.from_pandas(df)\n",
    "del L7,L8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot\n",
    "slope, intercept, r_value, p_value, std_err = stats.linregress(df.L7.values, df.L8.values)\n",
    "slope1, intercept1, r_value1, p_value1, std_err1 = stats.linregress(df.L8.values, df.L7.values)\n",
    "# slope, intercept, r_value, p_value, std_err = stats.linregress(np.log(L7),np.log(L8))\n",
    "\n",
    "# Compute regression type 2\n",
    "rma_results = regress2(df.L7.values, df.L8.values, _method_type_2=\"reduced major axis\")\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,7))\n",
    "plt.rc('font', size=BIGGER_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=BIGGER_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=BIGGER_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=BIGGER_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=BIGGER_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=BIGGER_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "\n",
    "# plot regression model\n",
    "plt.plot(np.array([0,1]), slope * np.array([0,1]) + intercept, color='red') # ols regression etm+ vs oli\n",
    "plt.plot([0,1], [0,1], color = 'white') # reference line\n",
    "# plt.plot(np.array([0,1]), (np.array([0,1]) - intercept1) / slope1 + intercept, color='peru') # ols regression oli vs etm+ \n",
    "plt.xlim(0, 1)\n",
    "plt.ylim(0, 1)\n",
    "# plt.axis('equal')\n",
    "k = rma_results['slope']\n",
    "b = rma_results['intercept'] \n",
    "plt.plot(np.array([0,1]), k * np.array([0,1]) + b, color='black') # rma regression\n",
    "\n",
    "ax.annotate('n:%.0f' % (len(df.L8.values)), xy=(0.7, 0.1),  xycoords='data',\n",
    "            horizontalalignment='left', verticalalignment='top',\n",
    "            )\n",
    "\n",
    "# scatter plot\n",
    "df.viz.heatmap(df.L7, df.L8, what=np.log(vx.stat.count()), show=True, colormap=\"viridis\", vmin=0.0, vmax=6,\n",
    "    xlabel='SWIR1 L7 Surface Reflectance', ylabel='SWIR1 L8 Surface Reflectance')\n",
    "ax.set_aspect('equal', 'box')    \n",
    "fig.savefig('print/L7L8/SWIR1.jpg', dpi=300)  \n",
    "\n",
    "print('OLS: \\ny={0:.4f}x+{1:.4f}\\nOLS_r:{2:.2f}, p:{3:.3f}, n:{4:.0f}'.format(slope,intercept,r_value,p_value, len(df.L8.values)))\n",
    "print('RMA: \\ny={0:.4f}x+{1:.4f}\\nRMA_r:{2:.2f}'.format(k,b, rma_results['r']))\n",
    "print('RMSE is %.4f' % (mean_squared_error(df.L8.values, df.L7.values, squared=False)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sns.residplot(x=df.L7.values, y=df.L8.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.viz.histogram(x=df.L7, label='L7')   \n",
    "df.viz.histogram(x=df.L8, label='L8')   \n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "swir2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "folderPath = r\"/data/shunan/data/harmonize_data/201305_08landsat/swir2\"\n",
    "L8, L7 = getBand(folderPath)\n",
    "df = pd.DataFrame({'L8': L8, 'L7': L7})\n",
    "df = vx.from_pandas(df)\n",
    "del L7,L8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot\n",
    "slope, intercept, r_value, p_value, std_err = stats.linregress(df.L7.values, df.L8.values)\n",
    "slope1, intercept1, r_value1, p_value1, std_err1 = stats.linregress(df.L8.values, df.L7.values)\n",
    "# slope, intercept, r_value, p_value, std_err = stats.linregress(np.log(L7),np.log(L8))\n",
    "\n",
    "# Compute regression type 2\n",
    "rma_results = regress2(df.L7.values, df.L8.values, _method_type_2=\"reduced major axis\")\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,7))\n",
    "plt.rc('font', size=BIGGER_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=BIGGER_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=BIGGER_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=BIGGER_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=BIGGER_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=BIGGER_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "\n",
    "# plot regression model\n",
    "plt.plot(np.array([0,1]), slope * np.array([0,1]) + intercept, color='red') # ols regression etm+ vs oli\n",
    "plt.plot([0,1], [0,1], color = 'white') # reference line\n",
    "# plt.plot(np.array([0,1]), (np.array([0,1]) - intercept1) / slope1 + intercept, color='peru') # ols regression oli vs etm+ \n",
    "plt.xlim(0, 1)\n",
    "plt.ylim(0, 1)\n",
    "# plt.axis('equal')\n",
    "k = rma_results['slope']\n",
    "b = rma_results['intercept'] \n",
    "plt.plot(np.array([0,1]), k * np.array([0,1]) + b, color='black') # rma regression\n",
    "\n",
    "ax.annotate('n:%.0f' % (len(df.L8.values)), xy=(0.7, 0.1),  xycoords='data',\n",
    "            horizontalalignment='left', verticalalignment='top',\n",
    "            )\n",
    "\n",
    "# scatter plot\n",
    "df.viz.heatmap(df.L7, df.L8, what=np.log(vx.stat.count()), show=True, colormap=\"viridis\", vmin=0.0, vmax=6,\n",
    "    xlabel='SWIR2 L7 Surface Reflectance', ylabel='SWIR2 L8 Surface Reflectance')\n",
    "ax.set_aspect('equal', 'box')    \n",
    "fig.savefig('print/L7L8/SWIR2.jpg', dpi=300)      \n",
    "\n",
    "print('OLS: \\ny={0:.4f}x+{1:.4f}\\nOLS_r:{2:.2f}, p:{3:.3f}, n:{4:.0f}'.format(slope,intercept,r_value,p_value, len(df.L8.values)))\n",
    "print('RMA: \\ny={0:.4f}x+{1:.4f}\\nRMA_r:{2:.2f}'.format(k,b, rma_results['r']))\n",
    "print('RMSE is %.4f' % (mean_squared_error(df.L8.values, df.L7.values, squared=False)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sns.residplot(x=df.L7.values, y=df.L8.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.viz.histogram(x=df.L7, label='L7')   \n",
    "df.viz.histogram(x=df.L8, label='L8')   \n",
    "plt.legend()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.12 ('earth')",
   "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.8.12"
  },
  "orig_nbformat": 2,
  "vscode": {
   "interpreter": {
    "hash": "2545903284817d58326c2b6cd0aa63ace3d30c673b9697d6672225518a7499b7"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}