{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os \n",
    "os.environ[\"ISISROOT\"] = \"/usgs/cpkgs/anaconda3_linux/envs/isis3.9.0\"\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from plio.io.io_gdal import GeoDataset \n",
    "import numpy as np\n",
    "\n",
    "from autocnet.cg import change_detection as cd\n",
    "from autocnet.utils.utils import bytescale \n",
    "\n",
    "from plio.io.io_gdal import GeoDataset\n",
    "import kalasiris as isis\n",
    "from subprocess import CalledProcessError\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "roi1_proj = \"/scratch/krodriguez/hirise/ESP_031022_1415_RED1.stitched.norm.34517_37165.proj.cub\"\n",
    "roi2_proj = \"/scratch/krodriguez/hirise/ESP_014038_1415_RED0.stitched.norm.30421_33069.proj.cub\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Run to create the subsolar azimuth band.  There is currently no option to output directly as array\n",
    "try:\n",
    "    isis.phocube(from_=roi1_proj, to='out.cub', subsolargroundazimuth=True)\n",
    "except CalledProcessError as e:\n",
    "    print(' '.join(e.cmd))\n",
    "    print(e.stdout\n",
    "    print(e.stderr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ssa = GeoDataset('out.cub')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ssa_array = ssa.read_array()\n",
    "ssa_slice = ssa_array[700:1200,1100:1500]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x2b32db435d50>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 254820 characters]\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "roi1_proj_geo = GeoDataset(roi1_proj)\n",
    "roi2_proj_geo = GeoDataset(roi2_proj)\n",
    "isis_null = isis.specialpixels.Real.Null\n",
    "\n",
    "proj_array1 = roi1_proj_geo.read_array()\n",
    "proj_array2 = roi2_proj_geo.read_array()\n",
    "\n",
    "proj_array1[proj_array1 == isis_null] = np.nan\n",
    "proj_array2[proj_array2 == isis_null] = np.nan\n",
    "\n",
    "plt.figure(0, figsize=(10, 10))\n",
    "plt.imshow(proj_array1, alpha=.4)\n",
    "plt.imshow(proj_array2, alpha=.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "slice1 = proj_array1[700:1200,1100:1500]\n",
    "slice2 = proj_array2[700:1200,1100:1500]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "t1,bdiff = cd.blob_detector(slice1, slice2, np.mean(ssa_slice), cd.image_diff, max_sigma=30, \n",
    "                         num_sigma=5, threshold=.45, angle_tolerance=20, overlap=.5, dist_upper_bound=10, subtractive=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 489424 characters]\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(0, figsize=(10, 10))\n",
    "axis = fig.add_subplot(1,1,1)\n",
    "plt.imshow(bdiff, alpha=1)\n",
    "for pt1 in t1.geometry:\n",
    "    x, y = pt1.xy\n",
    "    c = plt.Circle((x[0], y[0]), 1, color='red', linewidth=2, fill=False)\n",
    "    axis.add_patch(c)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "AutocnetDev",
   "language": "python",
   "name": "autocnetdev"
  },
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
