{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using PLIO to analyze control networks\n",
    "PLIO is a general purpose library for reading data from various sources. In this workshop, we will be using PLIO's ability to read ISIS control networks into a Pandas dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# PLIO uses pysis for some other things. We don't technically need this but it avoids a warning.\n",
    "import os\n",
    "os.environ['ISISROOT'] = '/usgs/cpkgs/anaconda3_linux/envs/isis4.3.0'\n",
    "os.environ['ISISDATA'] = '/usgs/cpkgs/isis3/isis_data'\n",
    "\n",
    "# 3D plotting toolkit for matplotlib\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "# Numerical Python library\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Our networks\n",
    "All of this data was generously provided by Lynn Weller and Mike Bland from their Europa control project.\n",
    "\n",
    "The first network is a very rough starting point. The Galileo images of Europa were put through the [findfeatures](https://isis.astrogeology.usgs.gov/Application/presentation/Tabbed/findfeatures/findfeatures.html) application and then all of the resulting networks were merged together. This network has many known issues including islands, massive residuals, and poor coverage.\n",
    "\n",
    "The second network is the final network containing Galileo and Voyager images of Europa. The issues from the initial network have been resolved and the final point cloud covers the majority of the body."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "galileo_net = '/scratch/jmapel/europa/networks/GLL_FFCombined_thin_SubReg2_Del_2.net'\n",
    "final_net = '/scratch/jmapel/europa/networks/GalileoVoyager_Europa_Merged_2020_CilixFree.net'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The control network dataframe\n",
    "\n",
    "PLIO directly ingests the data from the control network file. Each row in the dataframe is a single control measure and each column is a field from the protobuf control network. The data for control points is stored implicitly in its measures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function is what reads a control network file\n",
    "from plio.io.io_controlnetwork import from_isis\n",
    "\n",
    "galileo_df = from_isis(galileo_net)\n",
    "galileo_df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: How many measures are there in the network? How many points are there in the network? How many images are there in the network?\n",
    "\n",
    "tip: use len(dataframe) to find the number of rows in a dataframe\n",
    "\n",
    "tip: use dataframe[\"columnName\"].nunique() to find the number of unique values in a column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data types\n",
    "The different columns of our dataframe store different types of data. The cell below shows all of the the data types in the dataframe. You can see all of the different possible datatypes for a dataframe in the [pandas docs](https://pandas.pydata.org/pandas-docs/stable/user_guide/basics.html#basics-dtypes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "galileo_df.dtypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the data types are straightforward. For example, the line and sample are 64-bit floats. Let's dig into the more unusual types."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**pointType, measureType, aprioriSurfPointSource, and aprioriRadiusSource** are 64 bit integers, but those integers correspond to enumerations. For example, a pointType of 2 means Free. See the tables below for all of the enumerations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "galileo_df[['pointType', 'measureType', 'aprioriSurfPointSource']].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>**pointType**</center>\n",
    "\n",
    "| Value | Name              |\n",
    "| ----: | :---------------- |\n",
    "| 0     | Tie (obsolete)    |\n",
    "| 1     | Ground (obsolete) |\n",
    "| 2     | Free              |\n",
    "| 3     | Constrained       |\n",
    "| 4     | Fixed             |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>**measureType**</center>\n",
    "\n",
    "| Value | Name               |\n",
    "| ----: | :----------------- |\n",
    "| 0     | Candidate          |\n",
    "| 1     | Manual             |\n",
    "| 2     | RegisteredPixel    |\n",
    "| 3     | RegisteredSubPixel |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>**aprioriSurfPointSource & aprioriRadiusSource **</center>\n",
    "\n",
    "| Value | Name              |\n",
    "| ----: | :---------------- |\n",
    "| 0     | None              |\n",
    "| 1     | User              |\n",
    "| 2     | AverageOfMeasures |\n",
    "| 3     | Reference         |\n",
    "| 4     | Ellipsoid         |\n",
    "| 5     | DEM               |\n",
    "| 6     | Basemap           |\n",
    "| 7     | BundleSolution    |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: Have any measure in this network been sub-pixel registered?\n",
    "\n",
    "tip: look at the measure types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**id, pointChoosername, pointDatetime, aprioriSurfPointSourceFile, aprioriRadiusSourceFile, serialnumber, measureChoosername, and measureDatetime** are all listed as objects but are simply strings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "galileo_df[['id', 'serialnumber', 'pointChoosername', 'pointDatetime', 'measureChoosername', 'measureDatetime']].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**adjustedCovar, pointLog, and measureLog** are more complicated. We will go over adjustedCovar later with the final Euroap network. pointLog is leftover from older network formats and can be ignored. measureLog contains information about the registration of the measure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "galileo_df.loc[1,'measureLog']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data availability\n",
    "Depending on how your network was generated and what processing has been done, many fields will not be set. If a numerical field has a value of 0, then it has not been set. For example, our network has not been bundle adjusted, so there are only a priori ground points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "galileo_df[['aprioriX', 'aprioriY', 'aprioriZ', 'adjustedX', 'adjustedY', 'adjustedZ']].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: Can you find all of the fields that are completely unset in our control network?\n",
    "\n",
    "tip: numerical fields default to 0, strings default to an empty string \"\", and boolean values default to False."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also check which columns are default programmaticaly. The following cell checks if all of the values in a column are a default value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(galileo_df==0).all() | (galileo_df==\"\").all() | (galileo_df==False).all()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Looking at a bundle adjusted control network\n",
    "\n",
    "Our Galileo network is interesting but networks have significantly more useful information in them after bundle adjustment. So, let's take a look at the final Europa network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "final_net_df = from_isis(final_net)\n",
    "final_net_df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: What fields are set in the bundle adjusted network that weren't previously?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyzing the measures\n",
    "The data in a control network dataframe is not always in the format we want to work with. The measure residuals are broken down into the line and sample residuals. The following cell computes the full magnitude of the residuals and adds it to the dataframe under the \"residualMag\" column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "final_net_df['residualMag'] = np.sqrt(final_net_df['sampleResidual']**2 + final_net_df['lineResidual']**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's plot the residuals and see if we can form any theories. The next cell imports matplotlib for plotting tools and then plots the residuals in terms of sample and line residual. Note that the color of points is based on the residual magnitude, whcih should give a nice bullseye effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# This allows us to interact with our plots. This must be set before importing pyplot\n",
    "%matplotlib notebook\n",
    "\n",
    "# General plotting library\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "resid_fig = plt.figure(figsize=(6, 6))\n",
    "resid_ax = resid_fig.add_subplot(111)\n",
    "resid_scatter = resid_ax.scatter(final_net_df['sampleResidual'], final_net_df['lineResidual'], c=final_net_df['residualMag'], marker='+')\n",
    "resid_ax.set_aspect('equal')\n",
    "plt.axhline(0, color='black')\n",
    "plt.axvline(0, color='black')\n",
    "resid_cbar = plt.colorbar(resid_scatter)\n",
    "resid_fig.suptitle('Bundle Adjusted Measure Residuals')\n",
    "resid_ax.set_xlabel('Sample Residual')\n",
    "resid_ax.set_ylabel('Line Residual')\n",
    "resid_cbar.set_label('Residual Magnitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also color our points based on other properties. Let's try and separate the measures out by mission. The serial numbers should help us so let's look at the serial numbers for all of our images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_net_df['serialnumber'].unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each serial number starts with the mission name, which makes separating them out easy. All we need to do is check if the beginning of the serial number matches our mission.\n",
    "\n",
    "The pd.DataFrame.str package allows us to do this type of string comparisons quickly and easily. Here we will use the DataFrame.str.startswith method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "final_galileo_df = final_net_df[final_net_df['serialnumber'].str.startswith('Galileo')]\n",
    "final_voyager1_df = final_net_df[final_net_df['serialnumber'].str.startswith('Voyager1')]\n",
    "final_voyager2_df = final_net_df[final_net_df['serialnumber'].str.startswith('Voyager2')]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's plot the measures and color them based on their mission."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "inst_resid_fig = plt.figure(figsize=(6, 6))\n",
    "inst_resid_ax = inst_resid_fig.add_subplot(111)\n",
    "inst_resid_ax.scatter(final_galileo_df['sampleResidual'], final_galileo_df['lineResidual'], color='Green', marker='+', alpha=0.25, label='Galileo')\n",
    "inst_resid_ax.scatter(final_voyager1_df['sampleResidual'], final_voyager1_df['lineResidual'], color='Red', marker='+', alpha=0.25, label='Voyager1')\n",
    "inst_resid_ax.scatter(final_voyager2_df['sampleResidual'], final_voyager2_df['lineResidual'], color='Blue', marker='+', alpha=0.25, label='Voyager2')\n",
    "inst_resid_ax.set_aspect('equal')\n",
    "plt.axhline(0, color='black')\n",
    "plt.axvline(0, color='black')\n",
    "plt.legend()\n",
    "inst_resid_fig.suptitle('Bundle Adjusted Measure Residuals by Mission')\n",
    "inst_resid_ax.set_xlabel('Sample Residual')\n",
    "inst_resid_ax.set_ylabel('Line Residual')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What can you say about the residuals for the different missions based on our plot?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: What the descriptive statistics for the residual magnitude of the Galileo measures? What about for Voyager 1 and Voyager 2?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_galileo_df['residualMag'].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_voyager1_df['residualMag'].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_voyager2_df['residualMag'].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Do you notice anything interesting about the residual magnitudes for the different instruments? How does this compare to what you noticed with the scatter plot?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can even test if the measure residuals are normally distributed. The following cell performs a chi-squared test to see if the residual magnitudes could reasonably come from a normal distribution. This is important because it will tell us if we have large blunders in our network or systematic error from something like a bad sensor model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Statistics library\n",
    "from scipy import stats\n",
    "\n",
    "alpha = 1e-3 # 99.999% confidence\n",
    "_, normal_test_result = stats.normaltest(final_voyager1_df['residualMag'])\n",
    "print(f'Chi-squared test statistic: {normal_test_result}')\n",
    "if (normal_test_result < alpha):\n",
    "    print(\"The residuals are normally distributed\")\n",
    "else:\n",
    "    print(\"The residuals may not be normally distributed\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyzing the points\n",
    "The information for control points is duplicated for each measure they have. So, the first step in looking at control point data is to extract only the data we want from the dataframe. This will make the dataframe easier to read and it will make things run quicker.\n",
    "\n",
    "To do this, we're going to first extract all of the columns with point data. Then, we're going extract the first measure from each point. After all is said and done, we will have a dataframe with columns related to the point info and only one row for each point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "point_columns = ['id',\n",
    "                 'pointType',\n",
    "                 'pointChoosername',\n",
    "                 'pointDatetime',\n",
    "                 'pointEditLock',\n",
    "                 'pointIgnore',\n",
    "                 'pointJigsawRejected',\n",
    "                 'aprioriSurfPointSource',\n",
    "                 'aprioriSurfPointSourceFile',\n",
    "                 'aprioriRadiusSource',\n",
    "                 'aprioriRadiusSourceFile',\n",
    "                 'latitudeConstrained',\n",
    "                 'longitudeConstrained',\n",
    "                 'radiusConstrained',\n",
    "                 'aprioriX',\n",
    "                 'aprioriY',\n",
    "                 'aprioriZ',\n",
    "                 'aprioriCovar',\n",
    "                 'adjustedX',\n",
    "                 'adjustedY',\n",
    "                 'adjustedZ',\n",
    "                 'adjustedCovar',\n",
    "                 'pointLog']\n",
    "final_points_df = final_net_df[point_columns].drop_duplicates('id')\n",
    "final_points_df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we're going to transform the point data so that it's more useful to us. This cell will take the (X, Y, Z) adjusted ground points and convert them to (lat, lon, radius) using a library called pyproj. pyproj is a very powerful projections library and can do many cartofraphic transformations and projections.\n",
    "\n",
    "**Note: This cell will generate a warning because we are using old pyproj.Proj calls which will eventually need to change. For now we can ignore the warning.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Projection library for switching between rectangular and latitudinal\n",
    "os.environ['PROJ_LIB'] = '/usgs/cpkgs/anaconda3_linux/envs/autocnet/share/proj'\n",
    "import pyproj\n",
    "\n",
    "# Compute the lat/lon/alt\n",
    "europa_radii = [1562600, 1560300, 1559500]\n",
    "ecef = pyproj.Proj(proj='geocent', a=europa_radii[0], b=europa_radii[1], c=europa_radii[2])\n",
    "lla = pyproj.Proj(proj='latlong', a=europa_radii[0], b=europa_radii[1], c=europa_radii[2])\n",
    "lon, lat, alt = pyproj.transform(ecef, lla, final_points_df['adjustedX'].values, final_points_df['adjustedY'].values, final_points_df['adjustedZ'].values, radians=True)\n",
    "\n",
    "# Store the data in the dataframe\n",
    "final_points_df['latitude'] = lat\n",
    "final_points_df['longitude'] = lon\n",
    "final_points_df['altitude'] = alt\n",
    "\n",
    "# We will also want the point radii\n",
    "final_points_df['radius'] = np.sqrt(final_points_df['adjustedX']**2 + final_points_df['adjustedY']**2 + final_points_df['adjustedZ']**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because of how we defined our projection, the latitude and longitude values will be in radians. Also, the longitude will be in 180 postiive East. You can change this by modifying how you use pyproj but that is outside of this workshop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "final_points_df[[\"latitude\", \"longitude\", \"altitude\", \"radius\"]].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: Convert the latitude and longitude from radians to degrees:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to how we computed the residual magnitude, we want to compute the average residual magnitude for each point. The following cell goes back to our original dataframe, computes the mean point by point, and then saves all of the results in our new dataframe.\n",
    "\n",
    "**Note: This cell can take a while to run because it has to re-access the dataframe for every point**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "final_points_df[\"averageResidual\"] = 0\n",
    "for point_id, group in final_net_df.groupby('id'):\n",
    "    final_points_df.loc[final_points_df.id == point_id, \"averageResidual\"] = group['residualMag'].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: What is the 95th percentile for the average residuals?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting the points\n",
    "Now that we have latitudes and longitudes for each point, we can generate some simple plots to look at them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "point_map = plt.figure(figsize=(10, 10))\n",
    "point_ax = point_map.add_subplot(111)\n",
    "point_ax.scatter(final_points_df[\"longitude\"], final_points_df[\"latitude\"], marker='+')\n",
    "point_map.suptitle('Control Points')\n",
    "point_ax.set_xlabel('Longitude')\n",
    "point_ax.set_ylabel('Latitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can also be helpful to color the points based on different values. The following cell draws the same plot but colors each point based on its average residual. Because the residuals are not uniformly distributed we also apply a lograithmic scale to the colors that you can see in the colorbar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "point_resid_map = plt.figure(figsize=(10, 10))\n",
    "point_resid_ax = point_resid_map.add_subplot(111)\n",
    "point_resid_norm = matplotlib.colors.LogNorm(vmax=final_points_df[\"averageResidual\"].max())\n",
    "point_resid_scatter = point_resid_ax.scatter(final_points_df[\"longitude\"], final_points_df[\"latitude\"], c=final_points_df[\"averageResidual\"], alpha=0.5, norm=point_resid_norm, marker='+', cmap=plt.get_cmap('plasma'))\n",
    "point_resid_cbar = plt.colorbar(point_resid_scatter)\n",
    "point_resid_map.suptitle('Control Points')\n",
    "point_resid_ax.set_xlabel('Longitude')\n",
    "point_resid_ax.set_ylabel('Latitude')\n",
    "point_resid_cbar.set_label('Average Residual Magnitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting individual points can be helpful getting a general idea for the distribution of the points, but it can be hard to interpret the data in area where there are many points all ontop of each other. So, let's combine near by points and determine the residual based on the region.\n",
    "\n",
    "To do this, we're going to bin the points into a regular grid across the latitude and longitude and then compute the mean within each bin.\n",
    "\n",
    "**Try changing the grid_step value and re-running the two cells**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "grid_step = 10\n",
    "\n",
    "final_points_df['lonBin'] = final_points_df['longitude'].apply(lambda x: [e for e in range(-180, 180, grid_step) if e <= x][-1])\n",
    "final_points_df['latBin'] = final_points_df['latitude'].apply(lambda x: [e for e in range(-90, 90, grid_step) if e <= x][-1])\n",
    "\n",
    "avg_resid_binned = final_points_df.groupby(['lonBin', 'latBin'])['averageResidual'].mean()\n",
    "\n",
    "filled_data = []\n",
    "for lon_bin in range(-180, 180, grid_step):\n",
    "    for lat_bin in range(-90, 90, grid_step):\n",
    "        try:\n",
    "            filled_data.append(avg_resid_binned.loc[lon_bin, lat_bin])\n",
    "        except:\n",
    "            filled_data.append(0)\n",
    "filled_data = np.array(filled_data).reshape((int(360/grid_step), int(180/grid_step))).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "avg_gridded = plt.figure(figsize=(10, 5))\n",
    "avg_gridded_ax = avg_gridded.add_subplot(111)\n",
    "avg_gridded_plot = avg_gridded_ax.imshow(filled_data, origin='lower', extent= [-180, 180, -90, 90], cmap=plt.get_cmap('plasma'))\n",
    "avg_gridded_ax.scatter(final_points_df[\"longitude\"], final_points_df[\"latitude\"], color='black', marker='+', alpha=0.1)\n",
    "avg_gridded_cbar = plt.colorbar(avg_gridded_plot)\n",
    "avg_gridded.suptitle('Average Residual by lat/lon grid')\n",
    "avg_gridded_ax.set_xlabel('Longitude')\n",
    "avg_gridded_ax.set_ylabel('Latitude')\n",
    "avg_gridded_cbar.set_label('Average Residual Magnitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3D Plotting\n",
    "2D plotting either requires these simple equal area projections or converting to another projection via pyproj. Instead, let's look at our data in true 3D.\n",
    "\n",
    "The following cell plots the same data as before but plots it in 3d instead of just a 2d projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resid_fig_3d = plt.figure(figsize=(10, 10))\n",
    "resid_ax_3d = resid_fig_3d.add_subplot(111, projection='3d')\n",
    "resid_plot_3d = resid_ax_3d.scatter(final_points_df['adjustedX'], final_points_df['adjustedY'], final_points_df['adjustedZ'], c=final_points_df[\"averageResidual\"], alpha=0.5, norm=point_resid_norm, marker='+', cmap=plt.get_cmap('plasma'))\n",
    "resid_cbar_3d = plt.colorbar(resid_plot_3d)\n",
    "resid_fig_3d.suptitle('3D Control Points')\n",
    "resid_cbar_3d.set_label('Average Residual Magnitude (pix)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "AutoCNet PLIO (workshop)",
   "language": "python",
   "name": "autocnet_workshop"
  },
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
