From 0b25407f8cf423aad09f37e2c9cfb77db6e016cc Mon Sep 17 00:00:00 2001
From: Adam Paquette <acp263@nau.edu>
Date: Mon, 7 May 2018 11:53:10 -0700
Subject: [PATCH] Moved all transformations, updated bin script, and notebook.

---
 bin/Socetnet2ISIS.py            |  69 ++++++++
 notebooks/Socet2ISIS.ipynb      |  47 +++--
 plio/spatial/transformations.py | 297 ++++++++++++++++++++++++++++++++
 3 files changed, 387 insertions(+), 26 deletions(-)
 create mode 100644 plio/spatial/transformations.py

diff --git a/bin/Socetnet2ISIS.py b/bin/Socetnet2ISIS.py
index 047e35f..a410606 100644
--- a/bin/Socetnet2ISIS.py
+++ b/bin/Socetnet2ISIS.py
@@ -1,4 +1,73 @@
 import os
+import warnings
 import numpy as np
 
+from plio.examples import get_path
 from plio.io.io_bae import read_atf, read_gpf, read_ipf
+from plio.spatial.transformations import *
+import plio.io.io_controlnetwork as cn
+
+import pandas as pd
+
+# TODO: Change script to potentially handle configuration files
+
+# Setup the at_file and path to cubes
+cub_path = '/Volumes/Blueman/'
+at_file = get_path('CTX_Athabasca_Middle_step0.atf')
+
+# Define ipf mapping to cubs
+image_dict = {'P01_001540_1889_XI_08N204W' : 'P01_001540_1889_XI_08N204W.lev1.cub',
+              'P01_001606_1897_XI_09N203W' : 'P01_001606_1897_XI_09N203W.lev1.cub',
+              'P02_001804_1889_XI_08N204W' : 'P02_001804_1889_XI_08N204W.lev1.cub',
+              'P03_002226_1895_XI_09N203W' : 'P03_002226_1895_XI_09N203W.lev1.cub',
+              'P03_002371_1888_XI_08N204W' : 'P03_002371_1888_XI_08N204W.lev1.cub',
+              'P19_008344_1894_XN_09N203W' : 'P19_008344_1894_XN_09N203W.lev1.cub',
+              'P20_008845_1894_XN_09N203W' : 'P20_008845_1894_XN_09N203W.lev1.cub'}
+
+##
+# End Config
+##
+
+# Read in and setup the atf dict of information
+atf_dict = read_atf(at_file)
+
+# Get the gpf and ipf files using atf dict
+gpf_file = os.path.join(atf_dict['PATH'], atf_dict['GP_FILE']);
+ipf_list = [os.path.join(atf_dict['PATH'], i) for i in atf_dict['IMAGE_IPF']]
+
+# Read in the gpf file and ipf file(s) into seperate dataframes
+gpf_df = read_gpf(gpf_file)
+ipf_df = read_ipf(ipf_list)
+
+# Check for differences between point ids using each dataframes
+# point ids as a reference
+gpf_pt_idx = pd.Index(pd.unique(gpf_df['point_id']))
+ipf_pt_idx = pd.Index(pd.unique(ipf_df['pt_id']))
+
+point_diff = ipf_pt_idx.difference(gpf_pt_idx)
+
+if len(point_diff) != 0:
+    warnings.warn("The following points found in ipf files missing from gpf file: \n\n{}. \
+                  \n\nContinuing, but these points will be missing from the control network".format(list(point_diff)))
+
+# Merge the two dataframes on their point id columns
+socet_df = ipf_df.merge(gpf_df, left_on='pt_id', right_on='point_id')
+
+# Apply the transformations
+apply_transformations(atf_dict, socet_df)
+
+# Define column remap for socet dataframe
+column_remap = {'l.': 'y', 's.': 'x',
+                'res_l': 'LineResidual', 'res_s': 'SampleResidual', 'known': 'Type',
+                'lat_Y_North': 'AprioriY', 'long_X_East': 'AprioriX', 'ht': 'AprioriZ',
+                'sig0': 'AprioriLatitudeSigma', 'sig1': 'AprioriLongitudeSigma', 'sig2': 'AprioriRadiusSigma'}
+
+# Rename the columns using the column remap above
+socet_df.rename(columns = column_remap, inplace=True)
+
+images = pd.unique(socet_df['ipf_file'])
+
+serial_dict = serial_numbers(image_dict, cub_path)
+
+# creates the control network
+cn.to_isis('/Volumes/Blueman/test.net', socet_df, serial_dict)
diff --git a/notebooks/Socet2ISIS.ipynb b/notebooks/Socet2ISIS.ipynb
index 1d7681a..f9ff8e2 100644
--- a/notebooks/Socet2ISIS.ipynb
+++ b/notebooks/Socet2ISIS.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -26,7 +26,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -224,49 +224,44 @@
     "    return socet_df\n",
     "\n",
     "# creates a dict of serial numbers with the cub being the key\n",
-    "def serial_numbers(images, path, extension):\n",
+    "def serial_numbers(image_dict, path):\n",
     "    serial_dict = dict()\n",
-    "    \n",
-    "    for image in images:\n",
-    "        snum = sn.generate_serial_number(os.path.join(path, image + extension))\n",
+    "\n",
+    "    for key in image_dict:\n",
+    "        snum = sn.generate_serial_number(os.path.join(path, image_dict[key]))\n",
     "        snum = snum.replace('Mars_Reconnaissance_Orbiter', 'MRO')\n",
-    "        serial_dict[image] = snum\n",
+    "        serial_dict[key] = snum\n",
     "    return serial_dict"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": null,
    "metadata": {
     "scrolled": false
    },
-   "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "/Users/adampaquette/anaconda/envs/pysat/lib/python3.6/site-packages/ipykernel_launcher.py:173: UserWarning: The following points found in ipf files missing from gpf file: \n",
-      "\n",
-      "['P03_002226_1895_XI_09N203W_15', 'P03_002226_1895_XI_09N203W_16', 'P03_002226_1895_XI_09N203W_17', 'P03_002226_1895_XI_09N203W_18', 'P03_002226_1895_XI_09N203W_19', 'P03_002226_1895_XI_09N203W_20', 'P03_002226_1895_XI_09N203W_21', 'P03_002226_1895_XI_09N203W_22', 'P03_002226_1895_XI_09N203W_24', 'P03_002226_1895_XI_09N203W_26', 'P03_002226_1895_XI_09N203W_30', 'P03_002226_1895_XI_09N203W_31', 'P03_002226_1895_XI_09N203W_32', 'P03_002226_1895_XI_09N203W_34', 'P03_002226_1895_XI_09N203W_36', 'P03_002226_1895_XI_09N203W_37', 'P03_002226_1895_XI_09N203W_44', 'P03_002226_1895_XI_09N203W_48', 'P03_002226_1895_XI_09N203W_49', 'P03_002226_1895_XI_09N203W_56', 'P03_002226_1895_XI_09N203W_57', 'P03_002226_1895_XI_09N203W_61', 'P03_002226_1895_XI_09N203W_62', 'P03_002226_1895_XI_09N203W_63', 'P03_002226_1895_XI_09N203W_65', 'P19_008344_1894_XN_09N203W_4', 'P20_008845_1894_XN_09N203W_15'].                       \n",
-      "\n",
-      "Continuing, but these points will be missing from the control network\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "# Setup stuffs for the cub information namely the path and extension\n",
     "path = '/Volumes/Blueman/'\n",
-    "extension = '.lev1.cub'\n",
-    "prj_file = get_path('CTX_Athabasca_Middle_step0.atf')\n",
+    "atf_file = get_path('CTX_Athabasca_Middle_step0.atf')\n",
+    "\n",
+    "image_dict = {'P01_001540_1889_XI_08N204W' : 'P01_001540_1889_XI_08N204W.lev1.cub',\n",
+    "              'P01_001606_1897_XI_09N203W' : 'P01_001606_1897_XI_09N203W.lev1.cub',\n",
+    "              'P02_001804_1889_XI_08N204W' : 'P02_001804_1889_XI_08N204W.lev1.cub',\n",
+    "              'P03_002226_1895_XI_09N203W' : 'P03_002226_1895_XI_09N203W.lev1.cub',\n",
+    "              'P03_002371_1888_XI_08N204W' : 'P03_002371_1888_XI_08N204W.lev1.cub',\n",
+    "              'P19_008344_1894_XN_09N203W' : 'P19_008344_1894_XN_09N203W.lev1.cub',\n",
+    "              'P20_008845_1894_XN_09N203W' : 'P20_008845_1894_XN_09N203W.lev1.cub'}\n",
     "\n",
-    "socet_df = socet2isis(prj_file)\n",
+    "socet_df = socet2isis(atf_file)\n",
     "\n",
     "images = pd.unique(socet_df['ipf_file'])\n",
     "\n",
-    "serial_dict = serial_numbers(images, path, extension)\n",
+    "serial_dict = serial_numbers(image_dict, path)\n",
     "\n",
     "# creates the control network\n",
-    "cn.to_isis('/Volumes/Blueman/cn.net', socet_df, serial_dict)"
+    "cn.to_isis('/Volumes/Blueman/banana.net', socet_df, serial_dict)"
    ]
   },
   {
diff --git a/plio/spatial/transformations.py b/plio/spatial/transformations.py
new file mode 100644
index 0000000..18160ee
--- /dev/null
+++ b/plio/spatial/transformations.py
@@ -0,0 +1,297 @@
+import os
+import math
+import pyproj
+
+import plio.io.isis_serial_number as sn
+
+def line_sample_size(record, path):
+    """
+    Converts columns l. and s. to sample size, line size, and generates an
+    image index
+
+    Parameters
+    ----------
+    record : object
+             Pandas series object
+
+    path : str
+           Path to the associated sup files for a socet project
+
+    Returns
+    -------
+    : list
+      A list of sample_size, line_size, and img_index
+    """
+    with open(os.path.join(path, record['ipf_file'] + '.sup')) as f:
+        for i, line in enumerate(f):
+            if i == 2:
+                img_index = line.split('\\')
+                img_index = img_index[-1].strip()
+                img_index = img_index.split('.')[0]
+
+            if i == 3:
+                line_size = line.split(' ')
+                line_size = line_size[-1].strip()
+                assert int(line_size) > 0, "Line number {} from {} is a negative number: Invalid Data".format(line_size, record['ipf_file'])
+
+            if i == 4:
+                sample_size = line.split(' ')
+                sample_size = sample_size[-1].strip()
+                assert int(sample_size) > 0, "Sample number {} from {} is a negative number: Invalid Data".format(sample_size, record['ipf_file'])
+                break
+
+
+        line_size = int(line_size)/2.0 + record['l.'] + 1
+        sample_size = int(sample_size)/2.0 + record['s.'] + 1
+        return sample_size, line_size, img_index
+
+# converts known to ISIS keywords
+def known(record):
+    """
+    Converts the known field from a socet dataframe into the
+    isis point_type column
+
+    Parameters
+    ----------
+    record : object
+             Pandas series object
+
+    Returns
+    -------
+    : str
+      String representation of a known field
+    """
+    if record['known'] == 0:
+        return 'Free'
+
+    elif record['known'] == 1 or record['known'] == 2 or record['known'] == 3:
+        return 'Constrained'
+
+# converts +/- 180 system to 0 - 360 system
+def to_360(num):
+    """
+    Transforms a given number into 0 - 360 space
+
+    Parameters
+    ----------
+    num : int
+          A given integer
+
+    Returns
+    -------
+    : int
+      num moduloed by 360
+    """
+    return num % 360
+
+def oc2og(dlat, dMajorRadius, dMinorRadius):
+    """
+    Ocentric to ographic latitudes
+
+    Parameters
+    ----------
+    dlat : float
+           Latitude to convert
+
+    dMajorRadius : float
+                   Radius from the center of the body to the equater
+
+    dMinorRadius : float
+                   Radius from the pole to the center of mass
+
+    Returns
+    -------
+    dlat : float
+           Converted latitude into ographic space
+    """
+    try:
+        dlat = math.radians(dlat)
+        dlat = math.atan(((dMajorRadius / dMinorRadius)**2) * (math.tan(dlat)))
+        dlat = math.degrees(dlat)
+    except:
+        print ("Error in oc2og conversion")
+    return dlat
+
+def og2oc(dlat, dMajorRadius, dMinorRadius):
+    """
+    Ographic to ocentric latitudes
+
+    Parameters
+    ----------
+    dlat : float
+           Latitude to convert
+
+    dMajorRadius : float
+                   Radius from the center of the body to the equater
+
+    dMinorRadius : float
+                   Radius from the pole to the center of mass
+
+    Returns
+    -------
+    dlat : float
+           Converted latitude into ocentric space
+    """
+    try:
+        dlat = math.radians(dlat)
+        dlat = math.atan((math.tan(dlat) / ((dMajorRadius / dMinorRadius)**2)))
+        dlat = math.degrees(dlat)
+    except:
+        print ("Error in og2oc conversion")
+    return dlat
+
+# gets eRadius and pRadius from a .prj file
+def get_axis(file):
+    """
+    Gets eRadius and pRadius from a .prj file
+
+    Parameters
+    ----------
+    file : str
+           file with path to a given socet project file
+
+    Returns
+    -------
+    : list
+      A list of the eRadius and pRadius of the project file
+    """
+    with open(file) as f:
+        from collections import defaultdict
+
+        files = defaultdict(list)
+
+        for line in f:
+
+            ext = line.strip().split(' ')
+            files[ext[0]].append(ext[-1])
+
+        eRadius = float(files['A_EARTH'][0])
+        pRadius = eRadius * (1 - float(files['E_EARTH'][0]))
+
+        return eRadius, pRadius
+
+def lat_ISIS_coord(record, semi_major, semi_minor):
+    """
+    Function to convert lat_Y_North to ISIS_lat
+
+    Parameters
+    ----------
+    record : object
+             Pandas series object
+
+    semi_major : float
+                 Radius from the center of the body to the equater
+
+    semi_minor : float
+                 Radius from the pole to the center of mass
+
+    Returns
+    -------
+    coord_360 : float
+                Converted latitude into ocentric space, and mapped
+                into 0 to 360
+    """
+    ocentric_coord = og2oc(record['lat_Y_North'], semi_major, semi_minor)
+    coord_360 = to_360(ocentric_coord)
+    return coord_360
+
+def lon_ISIS_coord(record, semi_major, semi_minor):
+    """
+    Function to convert long_X_East to ISIS_lon
+
+    Parameters
+    ----------
+    record : object
+             Pandas series object
+
+    semi_major : float
+                 Radius from the center of the body to the equater
+
+    semi_minor : float
+                 Radius from the pole to the center of mass
+
+    Returns
+    -------
+    coord_360 : float
+                Converted longitude into ocentric space, and mapped
+                into 0 to 360
+    """
+    ocentric_coord = og2oc(record['long_X_East'], semi_major, semi_minor)
+    coord_360 = to_360(ocentric_coord)
+    return coord_360
+
+def body_fix(record, semi_major, semi_minor):
+    """
+    Transforms latitude, longitude, and height of a socet point into
+    a body fixed point
+
+    Parameters
+    ----------
+    record : object
+             Pandas series object
+
+    semi_major : float
+                 Radius from the center of the body to the equater
+
+    semi_minor : float
+                 Radius from the pole to the center of mass
+
+    Returns
+    -------
+    : list
+      Body fixed lat, lon and height coordinates as lat, lon, ht
+
+    """
+    ecef = pyproj.Proj(proj='geocent', a=semi_major, b=semi_minor)
+    lla = pyproj.Proj(proj='latlon', a=semi_major, b=semi_minor)
+    lon, lat, height = pyproj.transform(lla, ecef, record['long_X_East'], record['lat_Y_North'], record['ht'])
+    return lon, lat, height
+
+def apply_transformations(atf_dict, df):
+    """
+    Takes a atf dictionary and a socet dataframe and applies the necessary
+    transformations to convert that dataframe into a isis compatible
+    dataframe
+
+    Parameters
+    ----------
+    atf_dict : dict
+               Dictionary containing information from an atf file
+
+    df : object
+         Pandas dataframe object
+
+    """
+    prj_file = os.path.join(atf_dict['PATH'], atf_dict['PROJECT'].split('\\')[-1])
+
+    eRadius, pRadius = get_axis(prj_file)
+
+    df['s.'], df['l.'], df['image_index'] = (zip(*df.apply(line_sample_size, path = atf_dict['PATH'], axis=1)))
+    df['known'] = df.apply(known, axis=1)
+    df['lat_Y_North'] = df.apply(lat_ISIS_coord, semi_major = eRadius, semi_minor = pRadius, axis=1)
+    df['long_X_East'] = df.apply(lon_ISIS_coord, semi_major = eRadius, semi_minor = pRadius, axis=1)
+    df['long_X_East'], df['lat_Y_North'], df['ht'] = zip(*df.apply(body_fix, semi_major = eRadius, semi_minor = pRadius, axis = 1))
+
+def serial_numbers(image_dict, path):
+    """
+    Creates a dict of serial numbers with the cub being the key
+
+    Parameters
+    ----------
+    images : list
+
+    path : str
+
+    extension : str
+
+    Returns
+    -------
+    serial_dict : dict
+    """
+    serial_dict = dict()
+
+    for key in image_dict:
+        snum = sn.generate_serial_number(os.path.join(path, image_dict[key]))
+        snum = snum.replace('Mars_Reconnaissance_Orbiter', 'MRO')
+        serial_dict[key] = snum
+    return serial_dict
-- 
GitLab