Add essential functions from old transformations.py into socetnet2isi… (#90)

* Add essential functions from old transformations.py into socetnet2isis and isisnet2socet bin scripts

* Remove plio.spatial import statement from isisnet2socet

bin/isisnet2socet

@@ -2,14 +2,18 @@
 import os
 import math
 import argparse
+import warnings
+import pvl
+import math
+import pyproj
 
+import numpy as np
 import pandas as pd
 
 from plio.io.io_bae import save_gpf, save_ipf
 import plio.io.io_controlnetwork as cn
 import plio.io.isis_serial_number as sn
-from plio.spatial.transformations import apply_isis_transformations
-from plio.utils.utils import split_all_ext
+from plio.utils.utils import find_in_dict, split_all_ext
 
 def parse_args():
     parser = argparse.ArgumentParser()
@@ -25,6 +29,181 @@ def parse_args():
 
     return parser.parse_args()
 
+def reverse_known(record):
+    """
+    Converts the known field from an isis dataframe into the
+    socet known column
+
+    Parameters
+    ----------
+    record : object
+             Pandas series object
+
+    Returns
+    -------
+    : str
+      String representation of a known field
+    """
+    lookup = {0:0,
+              2:0,
+              1:3,
+              3:3,
+              4:3}
+    record_type = record['known']
+    return lookup[record_type]
+
+def reproject(record, semi_major, semi_minor, source_proj, dest_proj, **kwargs):
+    """
+    Thin wrapper around PyProj's Transform() function to transform 1 or more three-dimensional
+    point from one coordinate system to another. If converting between Cartesian
+    body-centered body-fixed (BCBF) coordinates and Longitude/Latitude/Altitude coordinates,
+    the values input for semi-major and semi-minor axes determine whether latitudes are
+    planetographic or planetocentric and determine the shape of the datum for altitudes.
+    If semi_major == semi_minor, then latitudes are interpreted/created as planetocentric
+    and altitudes are interpreted/created as referenced to a spherical datum.
+    If semi_major != semi_minor, then latitudes are interpreted/created as planetographic
+    and altitudes are interpreted/created as referenced to an ellipsoidal datum.
+
+    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
+
+    source_proj : str
+                         Pyproj string that defines a projection space ie. 'geocent'
+
+    dest_proj : str
+                      Pyproj string that defines a project space ie. 'latlon'
+
+    Returns
+    -------
+    : list
+      Transformed coordinates as y, x, z
+
+    """
+    source_pyproj = pyproj.Proj(proj = source_proj, a = semi_major, b = semi_minor)
+    dest_pyproj = pyproj.Proj(proj = dest_proj, a = semi_major, b = semi_minor)
+
+    y, x, z = pyproj.transform(source_pyproj, dest_pyproj, record[0], record[1], record[2], **kwargs)
+    return y, x, z
+
+def fix_sample_line(record, serial_dict, cub_dict):
+    """
+    Extracts the sample, line data from a cube and computes deviation from the
+    center of the image
+
+    Parameters
+    ----------
+    record : dict
+             Dict containing the key serialnumber, l., and s.
+
+    serial_dict : dict
+                  Maps serial numbers to images
+
+    cub_dict : dict
+                     Maps basic cub names to their assocated absoluate path cubs
+
+    Returns
+    -------
+    new_line : int
+               new line deviation from the center
+
+    new_sample : int
+                 new sample deviation from the center
+
+    """
+    # Cube location to load
+    cube = pvl.load(cub_dict[serial_dict[record['serialnumber']]])
+    line_size = find_in_dict(cube, 'Lines')
+    sample_size = find_in_dict(cube, 'Samples')
+
+    new_line = record['l.'] - (int(line_size / 2.0)) - 1
+    new_sample = record['s.'] - (int(sample_size / 2.0)) - 1
+
+    return new_line, new_sample
+
+def ignore_toggle(record):
+    """
+    Maps the stat column in a record to 0 or 1 based on True or False
+
+    Parameters
+    ----------
+    record : dict
+             Dict containing the key stat
+    """
+    if record['stat'] == True:
+        return 0
+    else:
+        return 1
+
+def apply_isis_transformations(df, eRadius, pRadius, serial_dict, cub_dict):
+    """
+    Takes an ISIS3 control network dataframe and applies the necessary
+    transformations to convert that dataframe into a Socet Set-compatible
+    dataframe
+
+    Parameters
+    ----------
+    df : object
+         Pandas dataframe object
+
+    eRadius : float
+              Equitorial radius of the target body
+
+    pRadius : float
+              Polar radius of the target body
+
+    serial_dict : dict
+                  Dictionary mapping serials as keys to images as the values
+
+    cub_dict : str
+               Dictionary mapping the basename of IPF files as keys to image cube names as values
+
+    """
+    # Convert from geocentered coords (x, y, z), to lat lon coords (latitude, longitude, alltitude)
+    ecef = np.array([[df['long_X_East']], [df['lat_Y_North']], [df['ht']]])
+    lla = reproject(ecef, semi_major = eRadius, semi_minor = pRadius,
+                           source_proj = 'geocent', dest_proj = 'latlong')
+
+    df['long_X_East'], df['lat_Y_North'], df['ht'] = lla[0][0], lla[1][0], lla[2][0]
+
+    # Convert longitude and latitude from degrees to radians
+    df['long_X_East'] = df['long_X_East'].apply(np.radians)
+    df['lat_Y_North'] = df['lat_Y_North'].apply(np.radians)
+
+    # Update the stat fields and add the val field as it is just a clone of stat
+    df['stat'] = df.apply(ignore_toggle, axis = 1)
+    df['val'] = df['stat']
+
+    # Update the known field, add the ipf_file field for saving, and
+    # update the line, sample using data from the cubes
+    df['known'] = df.apply(reverse_known, axis = 1)
+    df['ipf_file'] = df['serialnumber'].apply(lambda serial_number: serial_dict[serial_number])
+    df['l.'], df['s.'] = zip(*df.apply(fix_sample_line, serial_dict = serial_dict,
+                                       cub_dict = cub_dict, axis = 1))
+
+    # Add dummy for generic value setting
+    x_dummy = lambda x: np.full(len(df), x)
+
+    df['sig0'] = x_dummy(1)
+    df['sig1'] = x_dummy(1)
+    df['sig2'] = x_dummy(1)
+
+    df['res0'] = x_dummy(0)
+    df['res1'] = x_dummy(0)
+    df['res2'] = x_dummy(0)
+
+    df['fid_x'] = x_dummy(0)
+    df['fid_y'] = x_dummy(0)
+
+    df['no_obs'] = x_dummy(1)
+    df['fid_val'] = x_dummy(0)
 
 def main(args):
     # Create cub dict to map ipf to cub

bin/socetnet2isis

@@ -3,14 +3,17 @@ import os
 import sys
 import argparse
 import warnings
+import pvl
+import math
+import pyproj
 
+import numpy as np
 import pandas as pd
 
 from plio.io.io_bae import read_atf, read_gpf, read_ipf
 import plio.io.io_controlnetwork as cn
 import plio.io.isis_serial_number as sn
-from plio.spatial.transformations import apply_socet_transformations
-from plio.utils.utils import split_all_ext
+from plio.utils.utils import find_in_dict, split_all_ext
 
 def parse_args():
     parser = argparse.ArgumentParser()
@@ -23,7 +26,262 @@ def parse_args():
 
     return parser.parse_args()
 
+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
+
+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 * math.sqrt(1 - (float(files['E_EARTH'][0]) ** 2))
+
+        return eRadius, pRadius
+
+def reproject(record, semi_major, semi_minor, source_proj, dest_proj, **kwargs):
+    """
+    Thin wrapper around PyProj's Transform() function to transform 1 or more three-dimensional
+    point from one coordinate system to another. If converting between Cartesian
+    body-centered body-fixed (BCBF) coordinates and Longitude/Latitude/Altitude coordinates,
+    the values input for semi-major and semi-minor axes determine whether latitudes are
+    planetographic or planetocentric and determine the shape of the datum for altitudes.
+    If semi_major == semi_minor, then latitudes are interpreted/created as planetocentric
+    and altitudes are interpreted/created as referenced to a spherical datum.
+    If semi_major != semi_minor, then latitudes are interpreted/created as planetographic
+    and altitudes are interpreted/created as referenced to an ellipsoidal datum.
+
+    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
+
+    source_proj : str
+                         Pyproj string that defines a projection space ie. 'geocent'
+
+    dest_proj : str
+                      Pyproj string that defines a project space ie. 'latlon'
+
+    Returns
+    -------
+    : list
+      Transformed coordinates as y, x, z
+
+    """
+    source_pyproj = pyproj.Proj(proj = source_proj, a = semi_major, b = semi_minor)
+    dest_pyproj = pyproj.Proj(proj = dest_proj, a = semi_major, b = semi_minor)
+
+    y, x, z = pyproj.transform(source_pyproj, dest_pyproj, record[0], record[1], record[2], **kwargs)
+    return y, x, z
+
+# TODO: Does isis cnet need a convariance matrix for sigmas? Even with a static matrix of 1,1,1,1
+def compute_sigma_covariance_matrix(lat, lon, rad, latsigma, lonsigma, radsigma, semimajor_axis):
+
+    """
+    Given geospatial coordinates, desired accuracy sigmas, and an equitorial radius, compute a 2x3
+    sigma covariange matrix.
+    Parameters
+    ----------
+    lat : float
+          A point's latitude in degrees
+
+    lon : float
+          A point's longitude in degrees
+
+    rad : float
+          The radius (z-value) of the point in meters
+
+    latsigma : float
+               The desired latitude accuracy in meters (Default 10.0)
+
+    lonsigma : float
+               The desired longitude accuracy in meters (Default 10.0)
+
+    radsigma : float
+               The desired radius accuracy in meters (Defualt: 15.0)
+
+    semimajor_axis : float
+                     The semi-major or equitorial radius in meters (Default: 1737400.0 - Moon)
+    Returns
+    -------
+    rectcov : ndarray
+              (2,3) covariance matrix
+    """
+    lat = math.radians(lat)
+    lon = math.radians(lon)
+
+    # SetSphericalSigmasDistance
+    scaled_lat_sigma = latsigma / semimajor_axis
+
+    # This is specific to each lon.
+    scaled_lon_sigma = lonsigma * math.cos(lat) / semimajor_axis
+
+    # SetSphericalSigmas
+    cov = np.eye(3,3)
+    cov[0,0] = math.radians(scaled_lat_sigma) ** 2
+    cov[1,1] = math.radians(scaled_lon_sigma) ** 2
+    cov[2,2] = radsigma ** 2
+
+    # Approximate the Jacobian
+    j = np.zeros((3,3))
+    cosphi = math.cos(lat)
+    sinphi = math.sin(lat)
+    cos_lmbda = math.cos(lon)
+    sin_lmbda = math.sin(lon)
+    rcosphi = rad * cosphi
+    rsinphi = rad * sinphi
+    j[0,0] = -rsinphi * cos_lmbda
+    j[0,1] = -rcosphi * sin_lmbda
+    j[0,2] = cosphi * cos_lmbda
+    j[1,0] = -rsinphi * sin_lmbda
+    j[1,1] = rcosphi * cos_lmbda
+    j[1,2] = cosphi * sin_lmbda
+    j[2,0] = rcosphi
+    j[2,1] = 0.
+    j[2,2] = sinphi
+    mat = j.dot(cov)
+    mat = mat.dot(j.T)
+    rectcov = np.zeros((2,3))
+    rectcov[0,0] = mat[0,0]
+    rectcov[0,1] = mat[0,1]
+    rectcov[0,2] = mat[0,2]
+    rectcov[1,0] = mat[1,1]
+    rectcov[1,1] = mat[1,2]
+    rectcov[1,2] = mat[2,2]
+
+    return rectcov
+
+def compute_cov_matrix(record, semimajor_axis):
+    cov_matrix = compute_sigma_covariance_matrix(record['lat_Y_North'], record['long_X_East'], record['ht'], record['sig0'], record['sig1'], record['sig2'], semimajor_axis)
+    return cov_matrix.ravel().tolist()
+
+def stat_toggle(record):
+    if record['stat'] == 0:
+        return True
+    else:
+        return False
+
+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
+    """
+
+    lookup = {0: 'Free',
+              1: 'Constrained',
+              2: 'Constrained',
+              3: 'Constrained'}
+    return lookup[record['known']]
+
+def apply_socet_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'])
+
+    eRadius, pRadius = get_axis(prj_file)
+
+    # Convert longitude and latitude from radians to degrees
+    df['long_X_East'] = df['long_X_East'].apply(np.degrees)
+    df['lat_Y_North'] = df['lat_Y_North'].apply(np.degrees)
+
+    lla = np.array([[df['long_X_East']], [df['lat_Y_North']], [df['ht']]])
+
+    ecef = reproject(lla, semi_major = eRadius, semi_minor = pRadius,
+                              source_proj = 'latlon', dest_proj = 'geocent')
+
+    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['long_X_East'] = ecef[0][0]
+    df['lat_Y_North'] = ecef[1][0]
+    df['ht'] = ecef[2][0]
+    df['aprioriCovar'] = df.apply(compute_cov_matrix, semimajor_axis = eRadius, axis=1)
+    df['stat'] = df.apply(stat_toggle, axis=1)
+    
 def main(args):
     # Setup the at_file, path to cubes, and control network out path
     at_file = args.at_file