RADAR preparation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import statements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime\n",
    "from scipy import stats\n",
    "import time\n",
    "import pickle\n",
    "import pandas as pd\n",
    "import gpxpy\n",
    "from matplotlib import pyplot as plt\n",
    "import pyproj\n",
    "import georaster as gr\n",
    "import numpy as np\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## File structure\n",
    "The relations of each file to the specific flight are noted here.\n",
    "\n",
    "**radar_data** specifies the interpreted shot count - return time data. **tracks** specifies the track gpx (or custom pickled format) name and **wpt_sync** specifies which method to use for synchronizing the RADAR and GNSS time: using waypoint/shot count or GNSS time / shot count. **waypoints** specifies the gpx waypoint file to use, if that type of synchronization is chosen. **sync_file** specifies the relation between waypoint number, or GNSS time, and shot count. **doubled** is a boolean value for if the x-values in **radar_data** are stretched to twice their length (to help the previous interpretation part).\n",
    "\n",
    "A folder named **input** should be present in the same folder as the script. In this folder, five subfolders: /a_radar_data, /b_shot_vs_wpt, /b_waypoints, /c_tracks should be present, as well as **data_structure.csv**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>radar_data</th>\n",
       "      <th>tracks</th>\n",
       "      <th>wpt_sync</th>\n",
       "      <th>waypoints</th>\n",
       "      <th>sync_file</th>\n",
       "      <th>doubled</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>name</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>erik1</th>\n",
       "      <td>erik1_out_171212.txt</td>\n",
       "      <td>MARMA17C.gpx</td>\n",
       "      <td>True</td>\n",
       "      <td>Marma17C Waypoints.gpx</td>\n",
       "      <td>erik1_waypoints.csv</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>erik2</th>\n",
       "      <td>erik2_out_171212.txt</td>\n",
       "      <td>MARMA17C.gpx</td>\n",
       "      <td>True</td>\n",
       "      <td>Marma17C Waypoints.gpx</td>\n",
       "      <td>erik2_waypoints.csv</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>erik3</th>\n",
       "      <td>erik3_out_171212.txt</td>\n",
       "      <td>MARMA17C.gpx</td>\n",
       "      <td>True</td>\n",
       "      <td>Marma17C Waypoints.gpx</td>\n",
       "      <td>erik3_waypoints.csv</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>klara2</th>\n",
       "      <td>klara2_out_171212.txt</td>\n",
       "      <td>2014_radar.gpx</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "      <td>klara2_times.csv</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>anders1</th>\n",
       "      <td>anders1_out_171220.txt</td>\n",
       "      <td>anders1_track.p</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "      <td>anders1_times.csv</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     radar_data           tracks  wpt_sync  \\\n",
       "name                                                         \n",
       "erik1      erik1_out_171212.txt     MARMA17C.gpx      True   \n",
       "erik2      erik2_out_171212.txt     MARMA17C.gpx      True   \n",
       "erik3      erik3_out_171212.txt     MARMA17C.gpx      True   \n",
       "klara2    klara2_out_171212.txt   2014_radar.gpx     False   \n",
       "anders1  anders1_out_171220.txt  anders1_track.p     False   \n",
       "\n",
       "                      waypoints            sync_file  doubled  \n",
       "name                                                           \n",
       "erik1    Marma17C Waypoints.gpx  erik1_waypoints.csv    False  \n",
       "erik2    Marma17C Waypoints.gpx  erik2_waypoints.csv    False  \n",
       "erik3    Marma17C Waypoints.gpx  erik3_waypoints.csv    False  \n",
       "klara2                      NaN     klara2_times.csv     True  \n",
       "anders1                     NaN    anders1_times.csv     True  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "structure = pd.read_csv(\"input/data_structure.csv\", index_col=0)\n",
    "\n",
    "structure  # Preview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step A: Read RADAR data\n",
    "The RADAR data are structured as shot number (relative time), return time and return type. This step reads the radar data csv table and formats it to a Pandas DataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def a_read_radar(in_data, flight):\n",
    "    \"\"\"Reads radar csv tables and formats them to a DataFrame\n",
    "    :param flights: List of flights to process.\n",
    "    \"\"\"\n",
    "    doubled = structure.loc[flight, \"doubled\"]\n",
    "    \n",
    "    values = pd.read_csv(\"input/a_radar_data/\" + structure.loc[flight, \"radar_data\"])\n",
    "    \n",
    "    values.sort_values(\"x\", inplace=True)\n",
    "    values.y = values.y * -1  # The y-values should be positive\n",
    "\n",
    "    if doubled:  # If x-axis was doubled during interpretation\n",
    "        values.x = values.x * 0.5\n",
    "\n",
    "    return values[[\"x\", \"y\", \"typ\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Data preview**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "      <th>typ</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3.785150</td>\n",
       "      <td>136.966568</td>\n",
       "      <td>mark</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3.995588</td>\n",
       "      <td>135.359590</td>\n",
       "      <td>mark</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.206026</td>\n",
       "      <td>133.752611</td>\n",
       "      <td>mark</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.416463</td>\n",
       "      <td>132.145632</td>\n",
       "      <td>mark</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4.626901</td>\n",
       "      <td>130.538653</td>\n",
       "      <td>mark</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          x           y   typ\n",
       "0  3.785150  136.966568  mark\n",
       "1  3.995588  135.359590  mark\n",
       "2  4.206026  133.752611  mark\n",
       "3  4.416463  132.145632  mark\n",
       "4  4.626901  130.538653  mark"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "preview_flight = \"anders1\"\n",
    "a_preview_data = a_read_radar(in_data=None, flight=preview_flight)\n",
    "a_preview_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step B: Synchronize shot count with GNSS time\n",
    "The step uses the RADAR 'x-values' (shot count) and a waypoint file's respective timestamps, or hard coded GNSS times, to synchronize the RADAR data with the track file. The function uses linear regression to correlate these two parameters, which therefore accounts for possible inaccuracies in the **sync_file** due to human error. If this relation gives a Pearson correlation (if e.g. a faulty point exists in the sync file) lower than 0.9, a warning will be given. Having a good correlation is integral to georeferencing the data.\n",
    "\n",
    "The step also changes the 'typ' column to three separate columns representing return time (or lack thereof) of each respective type. The data are then resampled to 1s intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def b_shot_sync(in_data, flight, evaluation=False):\n",
    "    \"\"\"Uses a 'X-value to waypoint/time stamp' synchronisation file,\n",
    "    alternatively with a waypoint GPX-file, to synchronise shot count with the GNSS time.\"\"\"\n",
    "    \n",
    "    # Read structure table\n",
    "    waypoint_gpx = structure.loc[flight, \"waypoints\"]  # Name of waypoint gpx file, if needed\n",
    "    sync_file = structure.loc[flight, \"sync_file\"]  # Synchronisation file\n",
    "    wpt_sync = structure.loc[flight, \"wpt_sync\"]  # Boolean: If sync should be done with waypoints\n",
    "\n",
    "    # If synchronisation should be done with waypoints\n",
    "    if wpt_sync:\n",
    "        gpx = gpxpy.parse(open(f\"input/b_waypoints/{waypoint_gpx}\"))  # Load waypoints gpx\n",
    "        \n",
    "        # Get waypoints (name as int) and respective time (in Unix time)\n",
    "        wpts = pd.Series(\n",
    "            [time.mktime(wpt.time.timetuple()) \n",
    "             for wpt in gpx.waypoints], [int(wpt.name) for wpt in gpx.waypoints])\n",
    "        \n",
    "        # Load csv with x values for respective waypoint\n",
    "        wpt_x = pd.read_csv(f\"input/b_shot_vs_wpt/{sync_file}\", squeeze=True, index_col=0)  \n",
    "        x_times = wpt_x.apply(lambda x: wpts[x])  # Use waypoint names (indices) to get respective time\n",
    "\n",
    "    # If synchronisation should be done using time stamps\n",
    "    else:\n",
    "        x_times = pd.read_csv(\"input/b_shot_vs_wpt/\" + sync_file, squeeze=True, index_col=0)\n",
    "        \n",
    "        # Seconds from DateTime string\n",
    "        x_times = x_times.apply(\n",
    "            lambda x: time.mktime(datetime.datetime.strptime(x, \"%d/%m/%Y %H:%M:%S\").timetuple()))  \n",
    "\n",
    "    # If evaluation is True, return a graph of x-values in relation to time (Unix-time)\n",
    "    if evaluation:\n",
    "        x_times.plot(marker=\"+\")\n",
    "        plt.title(flight)\n",
    "        plt.show()\n",
    "    \n",
    "    # Check correlation between x-values and time, and warn if poor\n",
    "    corr = np.corrcoef(x_times, x_times.index)[1, 0]\n",
    "    if corr < 0.9:\n",
    "        print(f\"WARNING: {flight} sync file shows poor correlation (r={round(corr, 4)})\")\n",
    "\n",
    "    # Time where x=0; Regression between X-values and Unix-times.\n",
    "    # Extract m-value (where x=0) and convert back to DateTime.\n",
    "    x_0 = datetime.datetime.fromtimestamp(stats.linregress(x_times.index, list(x_times))[1])\n",
    "\n",
    "    # Make new column \"time\" using the x_0 time added with shot number / 2 (shoots twice a second)\n",
    "    in_data[\"time\"] = in_data[\"x\"].apply(lambda x: x_0 + datetime.timedelta(seconds=x / 2))\n",
    "    in_data = in_data.set_index(\"time\")[[\"y\", \"typ\"]]  # Set time as index, and only keep y and typ\n",
    "    \n",
    "\n",
    "    # Resample the data in 1s intervals, and make columns with depth for each signal.\n",
    "    types = pd.DataFrame(\n",
    "        {typ: typ_df.resample(\"1S\").mean().iloc[:, 0] for typ, typ_df in in_data.groupby(\"typ\")})\n",
    "    \n",
    "    return types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Data preview**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>fuffens</th>\n",
       "      <th>glacier</th>\n",
       "      <th>mark</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>time</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:18:57</th>\n",
       "      <td>146.004631</td>\n",
       "      <td>54.606126</td>\n",
       "      <td>200.920350</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:18:58</th>\n",
       "      <td>144.795451</td>\n",
       "      <td>53.805664</td>\n",
       "      <td>209.217710</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:18:59</th>\n",
       "      <td>143.586272</td>\n",
       "      <td>53.005203</td>\n",
       "      <td>217.389633</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:19:00</th>\n",
       "      <td>142.377093</td>\n",
       "      <td>52.204742</td>\n",
       "      <td>225.561556</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:19:01</th>\n",
       "      <td>141.908656</td>\n",
       "      <td>51.404280</td>\n",
       "      <td>232.990576</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        fuffens    glacier        mark\n",
       "time                                                  \n",
       "2012-04-16 09:18:57  146.004631  54.606126  200.920350\n",
       "2012-04-16 09:18:58  144.795451  53.805664  209.217710\n",
       "2012-04-16 09:18:59  143.586272  53.005203  217.389633\n",
       "2012-04-16 09:19:00  142.377093  52.204742  225.561556\n",
       "2012-04-16 09:19:01  141.908656  51.404280  232.990576"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b_preview_data = b_shot_sync(a_preview_data, flight=preview_flight)\n",
    "b_preview_data.dropna().head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step C: Georeference the RADAR data using a GNSS track\n",
    "With the correct time of the shots being known in the previous step, the GNSS and RADAR data can be joined using the time index. The GNSS track has to have the same frequency however, so it is consequently resampled to 1s intervals, using mean upsampling or linear interpolation for downsampling, or if gaps exist in the data. The GPX track input WGS 1984 coordinates are also projected into the specified projected system (epsg:3006 = SWEREF 99TM)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def c_georeference(in_data, flight):\n",
    "    \"\"\"Uses a 'track GPX' file to georeference the synchronised input data\"\"\"\n",
    "    \n",
    "    gpx_name = structure.loc[flight, \"tracks\"]\n",
    "    \n",
    "    if gpx_name.endswith(\".gpx\"):  # Load gpx file as DataFrame\n",
    "        gpx = gpxpy.parse(open(\"input/c_tracks/\" + gpx_name))\n",
    "\n",
    "        # Extract data for DataFrame\n",
    "        columns = ['Longitude', 'Latitude', 'Altitude']\n",
    "        gpx_df = pd.DataFrame(columns=columns)\n",
    "        for track in gpx.tracks:\n",
    "            for segment in track.segments:\n",
    "                for i, point in enumerate(segment.points):\n",
    "                    # Append to gpx_df\n",
    "                    gpx_df.loc[point.time, columns] = point.longitude, point.latitude, point.elevation  \n",
    "    \n",
    "    elif gpx_name.endswith(\".p\"):  # Exception made for 'anders1' flight. Loads prepared pickled dataframe\n",
    "        gpx_df = pickle.load(open(\"input/c_tracks/\" + gpx_name, \"rb\"))\n",
    "\n",
    "\n",
    "    # Project geographic coordinates to SWEREF99TM\n",
    "    in_proj = pyproj.Proj(init=\"epsg:4326\")  # WGS1984\n",
    "    out_proj = pyproj.Proj(init=\"epsg:3006\") # SWEREF99TM\n",
    "    \n",
    "    for i, row in gpx_df.iterrows():\n",
    "        x, y = pyproj.transform(in_proj, out_proj, row.Longitude, row.Latitude)\n",
    "        gpx_df.loc[i, \"Easting\"] = x\n",
    "        gpx_df.loc[i, \"Northing\"] = y\n",
    "    \n",
    "    # Mean downsampling (if freq < 1s), and interpolated gpx (if gpx is sampled over >1s)\n",
    "    gpx_i = gpx_df.resample(\"1S\").agg(lambda x: np.mean(x, axis=0)).interpolate()  \n",
    "   \n",
    "    joined = in_data.join(gpx_i)  # Join gpx data with radar data \n",
    "    \n",
    "    return joined[[\"Easting\", \"Northing\", \"Altitude\", \"mark\", \"glacier\", \"fuffens\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Data preview**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Easting</th>\n",
       "      <th>Northing</th>\n",
       "      <th>Altitude</th>\n",
       "      <th>mark</th>\n",
       "      <th>glacier</th>\n",
       "      <th>fuffens</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>time</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:18:57</th>\n",
       "      <td>654603.710526</td>\n",
       "      <td>7.556587e+06</td>\n",
       "      <td>1315.20</td>\n",
       "      <td>200.920350</td>\n",
       "      <td>54.606126</td>\n",
       "      <td>146.004631</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:18:58</th>\n",
       "      <td>654579.098396</td>\n",
       "      <td>7.556585e+06</td>\n",
       "      <td>1319.55</td>\n",
       "      <td>209.217710</td>\n",
       "      <td>53.805664</td>\n",
       "      <td>144.795451</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:18:59</th>\n",
       "      <td>654554.486266</td>\n",
       "      <td>7.556582e+06</td>\n",
       "      <td>1323.90</td>\n",
       "      <td>217.389633</td>\n",
       "      <td>53.005203</td>\n",
       "      <td>143.586272</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:19:00</th>\n",
       "      <td>654530.237682</td>\n",
       "      <td>7.556580e+06</td>\n",
       "      <td>1327.50</td>\n",
       "      <td>225.561556</td>\n",
       "      <td>52.204742</td>\n",
       "      <td>142.377093</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:19:01</th>\n",
       "      <td>654505.989098</td>\n",
       "      <td>7.556577e+06</td>\n",
       "      <td>1331.10</td>\n",
       "      <td>232.990576</td>\n",
       "      <td>51.404280</td>\n",
       "      <td>141.908656</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           Easting      Northing  Altitude        mark  \\\n",
       "time                                                                     \n",
       "2012-04-16 09:18:57  654603.710526  7.556587e+06   1315.20  200.920350   \n",
       "2012-04-16 09:18:58  654579.098396  7.556585e+06   1319.55  209.217710   \n",
       "2012-04-16 09:18:59  654554.486266  7.556582e+06   1323.90  217.389633   \n",
       "2012-04-16 09:19:00  654530.237682  7.556580e+06   1327.50  225.561556   \n",
       "2012-04-16 09:19:01  654505.989098  7.556577e+06   1331.10  232.990576   \n",
       "\n",
       "                       glacier     fuffens  \n",
       "time                                        \n",
       "2012-04-16 09:18:57  54.606126  146.004631  \n",
       "2012-04-16 09:18:58  53.805664  144.795451  \n",
       "2012-04-16 09:18:59  53.005203  143.586272  \n",
       "2012-04-16 09:19:00  52.204742  142.377093  \n",
       "2012-04-16 09:19:01  51.404280  141.908656  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c_preview_data = c_georeference(b_preview_data, flight=preview_flight)\n",
    "c_preview_data.dropna().head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step D: Shot return time to altitude conversion\n",
    "The step takes the return time values and converts them into reflection distance, and subsequently to elevation using the GNSS altitudes. Since the signal travels in a vastly different speed through air vs. in ice, this is accounted for using two different speeds for each respective setting.\n",
    "\n",
    "The distance to the glacier surface $d_g$ (if there is any at that specific shot) is first measured using air speed $v_a$:\n",
    "\n",
    "$$ d_g = t_g \\times v_a \\times \\frac{1}{512}$$\n",
    "\n",
    "The englacial signal distance $d_e$, if relevant, is then defined using the following formula:\n",
    "\n",
    "$$ d_e = d_g + v_e (t_e - t_g) \\times \\frac{1}{512}$$\n",
    "\n",
    "Where $v_e$ is the englacial signal velocity and $t_e$ is the englacial signal time. All distances are divided by 512 to get the distance in meters.\n",
    "Glacier travel time can also be expressed through glacier distance and air velocity, which is easier in the algorithm:\n",
    "\n",
    "$$t_g = \\frac{d_g \\times 512}{v_a} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def d_depth_from_y(in_data, flight):\n",
    "    \"\"\"Takes return time (y) values and converts them to GNSS derived altitude\"\"\"\n",
    "\n",
    "    # RADAR signal velocities\n",
    "    v0 = 300  # In air\n",
    "    v1 = 170  # In ice\n",
    "\n",
    "    def depth_ice(time, glacier_dist):  # In cases with glacier penetration (including air travel time)\n",
    "        t_g = (glacier_dist * 512) / v0\n",
    "        return glacier_dist + (v1 * (time - t_g)) / 512\n",
    "\n",
    "    def depth_air(val):  # In cases without glacier penetration.\n",
    "        return (val / 512) * v0\n",
    "\n",
    "    # Split dataset to data with or without glacier penetration\n",
    "    glacier_data = in_data[in_data[\"glacier\"].notnull()].copy()\n",
    "    mark_data =    in_data[in_data[\"glacier\"].isnull()].copy()\n",
    "    \n",
    "    # Remove any fuffens if it doesn't have a glacier value (needed for depth calculation)\n",
    "    mark_data.loc[:, \"fuffens\"] = np.NaN\n",
    "    \n",
    "    # Calculate depth for values featuring no glacier penetration (only through air)\n",
    "    glacier_data[\"glacier\"] = depth_air(glacier_data[\"glacier\"])\n",
    "    mark_data[\"mark\"]    =    depth_air(mark_data[\"mark\"])\n",
    "    \n",
    "    # Calculate depth for values featuring glacier penetration    \n",
    "    for col in [\"mark\", \"fuffens\"]:\n",
    "        glacier_data[col] = depth_ice(glacier_data[col], glacier_data[\"glacier\"])\n",
    "\n",
    "    # Concatenate the two split DataFrames\n",
    "    data = pd.concat([glacier_data, mark_data]).sort_index()\n",
    "\n",
    "    # Convert depths to altitudes by suptracting GNSS altitude with depth\n",
    "    data[[\"mark\", \"fuffens\", \"glacier\"]] = data[[\"mark\", \"fuffens\", \"glacier\"]].apply(\n",
    "        lambda x: data[\"Altitude\"] - x)\n",
    "\n",
    "    return data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Data preview**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Easting</th>\n",
       "      <th>Northing</th>\n",
       "      <th>Altitude</th>\n",
       "      <th>mark</th>\n",
       "      <th>glacier</th>\n",
       "      <th>fuffens</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>time</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:18:57</th>\n",
       "      <td>654603.710526</td>\n",
       "      <td>7.556587e+06</td>\n",
       "      <td>1315.20</td>\n",
       "      <td>1234.623328</td>\n",
       "      <td>1283.204223</td>\n",
       "      <td>1252.857063</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:18:58</th>\n",
       "      <td>654579.098396</td>\n",
       "      <td>7.556585e+06</td>\n",
       "      <td>1319.55</td>\n",
       "      <td>1236.421588</td>\n",
       "      <td>1288.023243</td>\n",
       "      <td>1257.811791</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:18:59</th>\n",
       "      <td>654554.486266</td>\n",
       "      <td>7.556582e+06</td>\n",
       "      <td>1323.90</td>\n",
       "      <td>1238.261496</td>\n",
       "      <td>1292.842264</td>\n",
       "      <td>1262.766518</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:19:00</th>\n",
       "      <td>654530.237682</td>\n",
       "      <td>7.556580e+06</td>\n",
       "      <td>1327.50</td>\n",
       "      <td>1239.351404</td>\n",
       "      <td>1296.911284</td>\n",
       "      <td>1266.971246</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:19:01</th>\n",
       "      <td>654505.989098</td>\n",
       "      <td>7.556577e+06</td>\n",
       "      <td>1331.10</td>\n",
       "      <td>1240.687980</td>\n",
       "      <td>1300.980304</td>\n",
       "      <td>1270.930024</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           Easting      Northing  Altitude         mark  \\\n",
       "time                                                                      \n",
       "2012-04-16 09:18:57  654603.710526  7.556587e+06   1315.20  1234.623328   \n",
       "2012-04-16 09:18:58  654579.098396  7.556585e+06   1319.55  1236.421588   \n",
       "2012-04-16 09:18:59  654554.486266  7.556582e+06   1323.90  1238.261496   \n",
       "2012-04-16 09:19:00  654530.237682  7.556580e+06   1327.50  1239.351404   \n",
       "2012-04-16 09:19:01  654505.989098  7.556577e+06   1331.10  1240.687980   \n",
       "\n",
       "                         glacier      fuffens  \n",
       "time                                           \n",
       "2012-04-16 09:18:57  1283.204223  1252.857063  \n",
       "2012-04-16 09:18:58  1288.023243  1257.811791  \n",
       "2012-04-16 09:18:59  1292.842264  1262.766518  \n",
       "2012-04-16 09:19:00  1296.911284  1266.971246  \n",
       "2012-04-16 09:19:01  1300.980304  1270.930024  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d_preview_data = d_depth_from_y(c_preview_data, flight=preview_flight)\n",
    "d_preview_data.dropna().head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step E: Altitude correction using DEM\n",
    "The GNSS altitude measurements considered to be unreliable, and the altitudes of every point featuring a glacier surface is corrected using a DEM. The average offset is also applied onto data points with no glacier surface. The ground is not used, due to many erroneous ground points seemingly existing, maybe due to the helicopter flying close to mountain sides.\n",
    "\n",
    "A possible error in this might be that every point is corrected using a single DEM, and if large melt has occurred within the study period, e.g. 2012-2017, there will be an equally large error if the englacial properties have not changed accordingly. A solution would be to use multiple DEM's, possibly with linear interpolation between the two. The errors might however not be large enough to make a difference, but they need to be considered."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def e_correct_alt(in_data, flight):\n",
    "\n",
    "    glacier_data = in_data[in_data[\"glacier\"].notnull()].copy()\n",
    "    mark_data =    in_data[in_data[\"glacier\"].isnull()].copy()\n",
    "\n",
    "    dem = gr.SingleBandRaster(\"input/Marma17_DEM_031017.tif\")\n",
    "\n",
    "    differences = []\n",
    "    for i, row in glacier_data.iterrows():\n",
    "        # Difference between GNSS Altitude + Depth and DEM Elevation\n",
    "        diff = dem.value_at_coords(row[\"Easting\"], row[\"Northing\"]) - row[\"glacier\"]  \n",
    "        differences.append(diff)\n",
    "\n",
    "        glacier_data.loc[i, [\"Altitude\", \"glacier\", \"mark\", \"fuffens\"]] += diff\n",
    "\n",
    "    # Apply mean offset to ground with no glacier present.\n",
    "    mark_data[[\"Altitude\", \"mark\"]] += np.mean(differences)\n",
    "\n",
    "    data = pd.concat([glacier_data, mark_data]).sort_index()\n",
    "\n",
    "    return data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Data preview**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Easting</th>\n",
       "      <th>Northing</th>\n",
       "      <th>Altitude</th>\n",
       "      <th>mark</th>\n",
       "      <th>glacier</th>\n",
       "      <th>fuffens</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>time</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:18:57</th>\n",
       "      <td>654603.710526</td>\n",
       "      <td>7.556587e+06</td>\n",
       "      <td>1394.685596</td>\n",
       "      <td>1314.108925</td>\n",
       "      <td>1362.689819</td>\n",
       "      <td>1332.342660</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:18:58</th>\n",
       "      <td>654579.098396</td>\n",
       "      <td>7.556585e+06</td>\n",
       "      <td>1397.876000</td>\n",
       "      <td>1314.747587</td>\n",
       "      <td>1366.349243</td>\n",
       "      <td>1336.137791</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:18:59</th>\n",
       "      <td>654554.486266</td>\n",
       "      <td>7.556582e+06</td>\n",
       "      <td>1400.907956</td>\n",
       "      <td>1315.269452</td>\n",
       "      <td>1369.850220</td>\n",
       "      <td>1339.774474</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:19:00</th>\n",
       "      <td>654530.237682</td>\n",
       "      <td>7.556580e+06</td>\n",
       "      <td>1404.080171</td>\n",
       "      <td>1315.931575</td>\n",
       "      <td>1373.491455</td>\n",
       "      <td>1343.551417</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2012-04-16 09:19:01</th>\n",
       "      <td>654505.989098</td>\n",
       "      <td>7.556577e+06</td>\n",
       "      <td>1407.221136</td>\n",
       "      <td>1316.809116</td>\n",
       "      <td>1377.101440</td>\n",
       "      <td>1347.051160</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           Easting      Northing     Altitude         mark  \\\n",
       "time                                                                         \n",
       "2012-04-16 09:18:57  654603.710526  7.556587e+06  1394.685596  1314.108925   \n",
       "2012-04-16 09:18:58  654579.098396  7.556585e+06  1397.876000  1314.747587   \n",
       "2012-04-16 09:18:59  654554.486266  7.556582e+06  1400.907956  1315.269452   \n",
       "2012-04-16 09:19:00  654530.237682  7.556580e+06  1404.080171  1315.931575   \n",
       "2012-04-16 09:19:01  654505.989098  7.556577e+06  1407.221136  1316.809116   \n",
       "\n",
       "                         glacier      fuffens  \n",
       "time                                           \n",
       "2012-04-16 09:18:57  1362.689819  1332.342660  \n",
       "2012-04-16 09:18:58  1366.349243  1336.137791  \n",
       "2012-04-16 09:18:59  1369.850220  1339.774474  \n",
       "2012-04-16 09:19:00  1373.491455  1343.551417  \n",
       "2012-04-16 09:19:01  1377.101440  1347.051160  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e_preview_data = e_correct_alt(d_preview_data, flight=preview_flight)\n",
    "e_preview_data.dropna().head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Main application loop\n",
    "Here all of the above steps are executed. It takes a list of flight names (RADAR series names), alternatively a single name or \"all\" (every name in the structure file). By default it writes a csv for each flight, and returns a dictionary {name: DataFrame}. If the **unified** parameter is True, it writes a single csv with the flight names as a separate column, and returns a respective DataFrame.\n",
    "\n",
    "The **tag** parameter sets a specific name for the output file (e.g. tag=\"newest\" -> unified_newest.csv). The default value is the date as YYMMDD of execution.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "auto_tag = datetime.datetime.today().strftime(\"%Y%m%d\")[2:]  # YYMMDD\n",
    "\n",
    "def prepare_radar(flights, save=False, tag=auto_tag, unified=False):\n",
    "\n",
    "    # If 'all' is entered, use every flight in the structure table indices\n",
    "    if flights == \"all\":\n",
    "        flights = list(structure.index.values)\n",
    "    # If input is a single flight as a string, convert it to a list for compatibility\n",
    "    elif isinstance(flights, str):\n",
    "        flights = [flights]\n",
    "\n",
    "    # Properties of the different steps\n",
    "    steps = [\"a\", \"b\", \"c\", \"d\", \"e\"]\n",
    "    functions = [a_read_radar, b_shot_sync, c_georeference, d_depth_from_y, e_correct_alt]\n",
    "    step_function = dict(zip(steps, functions))\n",
    "\n",
    "    ###################\n",
    "    # Perform steps\n",
    "    ###################\n",
    "    data = {}\n",
    "    flight_data = pd.DataFrame()  # Needed to not raise an error at A\n",
    "    for flight in flights:\n",
    "        for step in steps:\n",
    "            flight_data = step_function[step](flight_data, flight)\n",
    "        \n",
    "        # Add flight data to data dictionary\n",
    "        data[flight] = flight_data\n",
    "        \n",
    "        # Save data\n",
    "        if save and not unified:\n",
    "            flight_data.to_csv(f\"export/{flight}_{tag}.csv\")\n",
    "        \n",
    "        print(f\"{flight} finished.\")\n",
    "\n",
    "    # Option to return a DataFrame with a flight column, instead of dictionary with flight as a key\n",
    "    if unified:\n",
    "        all_data = pd.DataFrame()\n",
    "        for name, df in data.items():\n",
    "            df[\"Name\"] = name\n",
    "            all_data = pd.concat([all_data, df])\n",
    "        all_data.to_csv(f\"export/unified_{tag}.csv\")\n",
    "        \n",
    "        return all_data\n",
    "    \n",
    "    if len(flights) == 1:\n",
    "        return flight_data\n",
    "    else:\n",
    "        return data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Data preview**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "erik1 finished.\n",
      "erik2 finished.\n",
      "erik3 finished.\n",
      "klara2 finished.\n",
      "anders1 finished.\n",
      "Files in 'export/' folder: ['clouds', 'unified_180301.csv']\n"
     ]
    }
   ],
   "source": [
    "data = prepare_radar(\"all\", save=True, unified=True)\n",
    "print(\"Files in 'export/' folder:\", os.listdir(\"export/\")[1:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Point cloud preparation\n",
    "This function prepares three csv's (one for each return type) and prepares them to be easily read by point cloud softwares, such as CloudCompare. Color is also given appropriately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def export_points(tag=auto_tag):\n",
    "    data = pd.read_csv(f\"export/unified_{tag}.csv\", index_col=0)  # Reads unified csv with specified tag\n",
    "    \n",
    "    # Convert data names to integers (sorted alphabetically)\n",
    "    names = data[\"Name\"].unique().tolist()\n",
    "    names.sort()\n",
    "    data[\"Name\"] = data[\"Name\"].apply(lambda x: names.index(x))\n",
    "    \n",
    "    \n",
    "    colors = [(190, 121, 74), (56, 57, 255), (248, 57, 4)]  # Brown, blue and red for respective data type\n",
    "    \n",
    "    # Create one file for each data type\n",
    "    data.drop(\"Altitude\", axis=1, inplace=True)\n",
    "    for i, col in enumerate([\"mark\", \"glacier\", \"fuffens\"]):\n",
    "        csv = data.reset_index(drop=True)\n",
    "        \n",
    "        # Add color columns\n",
    "        for c_i, c in enumerate([\"r\", \"g\", \"b\"]):\n",
    "            csv[c] = colors[i][c_i]\n",
    "        \n",
    "        # Write csv with correctly ordered columns\n",
    "        csv = csv[[\"Easting\", \"Northing\", col, \"r\", \"g\", \"b\", \"Name\"]].dropna()\n",
    "        csv.to_csv(f\"export/clouds/radar_{col}_{tag}.csv\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Data Preview**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Files in 'export/clouds' folder: ['radar_fuffens_180301.csv', 'radar_glacier_180301.csv', 'radar_mark_180301.csv']\n"
     ]
    }
   ],
   "source": [
    "export_points()\n",
    "print(\"Files in 'export/clouds' folder:\", os.listdir(\"export/clouds\")[1:])"
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}