diff --git a/bin/socet2isis b/bin/socet2isis index a32aa0fab4d7982fd28a3596f33a2dd6fe061279..d89a6f94c16dad6f1a0514bb8d9828fcd8dd1dcb 100644 --- a/bin/socet2isis +++ b/bin/socet2isis @@ -74,10 +74,10 @@ def main(args): # 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'} + '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) diff --git a/notebooks/Socet2ISIS.ipynb b/notebooks/Socet2ISIS.ipynb index a75564baa800ed5f4ac3e93f8935b1fb843635fe..b060d18a23d1e74383bfa298c88bb5eccde2c375 100644 --- a/notebooks/Socet2ISIS.ipynb +++ b/notebooks/Socet2ISIS.ipynb @@ -2,22 +2,26 @@ "cells": [ { "cell_type": "code", - "execution_count": 56, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "import sys\n", - "from functools import singledispatch\n", - "import warnings\n", + "import csv\n", + "import pvl\n", "\n", "import pandas as pd\n", "import numpy as np\n", "import math\n", "import pyproj\n", + "from functools import singledispatch\n", + "import warnings\n", "\n", "from plio.examples import get_path\n", - "from plio.io.io_bae import read_gpf, read_ipf\n", + "from plio.io.io_bae import read_gpf, read_ipf, read_atf\n", + "from plio.utils.utils import find_in_dict\n", + "from plio.spatial.transformations import *\n", "from collections import defaultdict\n", "import plio.io.io_controlnetwork as cn\n", "import plio.io.isis_serial_number as sn" @@ -25,283 +29,27 @@ }, { "cell_type": "code", - "execution_count": 85, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# Reads a .atf file and outputs all of the \n", - "# .ipf, .gpf, .sup, .prj, and path to locate the \n", - "# .apf file (should be the same as all others) \n", - "def read_atf(atf_file):\n", - " with open(atf_file) as f:\n", - " \n", - " # Extensions of files we want\n", - " files_ext = ['.prj', '.sup', '.ipf', '.gpf']\n", - " files_dict = []\n", - " files = defaultdict(list)\n", - "\n", - " for line in f:\n", - " ext = os.path.splitext(line)[-1].strip()\n", - " \n", - " # Check is needed for split as all do not have a space\n", - " if ext in files_ext:\n", - " \n", - " # If it is the .prj file, it strips the directory away and grabs file name\n", - " if ext == '.prj':\n", - " files[ext].append(line.strip().split(' ')[1].split('\\\\')[-1])\n", - " \n", - " # If the ext is in the list of files we care about, it addes to the dict\n", - " files[ext].append(line.strip().split(' ')[-1])\n", - " \n", - " else:\n", - " \n", - " # Adds to the dict even if not in files we care about\n", - " files[ext.strip()].append(line)\n", - " \n", - " # Gets the base filepath\n", - " files['basepath'] = os.path.dirname(os.path.abspath(atf_file))\n", - " \n", - " # Creates a dict out of file lists for GPF, PRJ, IPF, and ATF\n", - " files_dict = (dict(files_dict))\n", - " \n", - " # Sets the value of IMAGE_IPF to all IPF images\n", - " files_dict['IMAGE_IPF'] = files['.ipf']\n", - " \n", - " # Sets the value of IMAGE_SUP to all SUP images\n", - " files_dict['IMAGE_SUP'] = files['.sup']\n", - " \n", - " # Sets value for GPF file\n", - " files_dict['GP_FILE'] = files['.gpf'][0]\n", - " \n", - " # Sets value for PRJ file\n", - " files_dict['PROJECT'] = files['.prj'][0]\n", - " \n", - " # Sets the value of PATH to the path of the ATF file\n", - " files_dict['PATH'] = files['basepath']\n", - " \n", - " return files_dict\n", - "\n", - "# converts columns l. and s. to isis\n", - "# no transform applied\n", - "def line_sample_size(record, path):\n", - " with open(os.path.join(path, record['ipf_file'] + '.sup')) as f:\n", - " for i, line in enumerate(f):\n", - " if i == 2:\n", - " img_index = line.split('\\\\')\n", - " img_index = img_index[-1].strip()\n", - " img_index = img_index.split('.')[0]\n", - " \n", - " if i == 3:\n", - " line_size = line.split(' ')\n", - " line_size = line_size[-1].strip()\n", - " assert int(line_size) > 0, \"Line number {} from {} is a negative number: Invalid Data\".format(line_size, record['ipf_file'])\n", - " \n", - " if i == 4:\n", - " sample_size = line.split(' ')\n", - " sample_size = sample_size[-1].strip()\n", - " assert int(sample_size) > 0, \"Sample number {} from {} is a negative number: Invalid Data\".format(sample_size, record['ipf_file'])\n", - " break\n", - " \n", - " \n", - " line_size = int(line_size)/2.0 + record['l.'] + 1\n", - " sample_size = int(sample_size)/2.0 + record['s.'] + 1\n", - " return sample_size, line_size, img_index\n", - " \n", - "# converts known to ISIS keywords\n", - "# transform\n", - "def known(record):\n", - " if record['known'] == 0:\n", - " return 'Free'\n", - " \n", - " elif record['known'] == 1 or record['known'] == 2 or record['known'] == 3:\n", - " return 'Constrained'\n", - " \n", - "# converts +/- 180 system to 0 - 360 system\n", - "def to_360(num):\n", - " return num % 360\n", - "\n", - "# ocentric to ographic latitudes\n", - "# transform but unsure how to handle\n", - "def oc2og(dlat, dMajorRadius, dMinorRadius):\n", - " try: \n", - " dlat = math.radians(dlat)\n", - " dlat = math.atan(((dMajorRadius / dMinorRadius)**2) * (math.tan(dlat)))\n", - " dlat = math.degrees(dlat)\n", - " except:\n", - " print (\"Error in oc2og conversion\")\n", - " return dlat\n", - "\n", - "# ographic to ocentric latitudes\n", - "# transform but unsure how to handle\n", - "def og2oc(dlat, dMajorRadius, dMinorRadius):\n", - " try:\n", - " dlat = math.radians(dlat)\n", - " dlat = math.atan((math.tan(dlat) / ((dMajorRadius / dMinorRadius)**2)))\n", - " dlat = math.degrees(dlat)\n", - " except:\n", - " print (\"Error in og2oc conversion\")\n", - " return dlat\n", - "\n", - "# gets eRadius and pRadius from a .prj file\n", - "def get_axis(file):\n", - " with open(file) as f:\n", - " from collections import defaultdict\n", - "\n", - " files = defaultdict(list)\n", - " \n", - " for line in f:\n", - " \n", - " ext = line.strip().split(' ')\n", - " files[ext[0]].append(ext[-1])\n", - " \n", - " eRadius = float(files['A_EARTH'][0])\n", - " pRadius = eRadius * (1 - float(files['E_EARTH'][0]))\n", - " \n", - " return eRadius, pRadius\n", - " \n", - "# function to convert lat_Y_North to ISIS_lat\n", - "def lat_ISIS_coord(record, semi_major, semi_minor):\n", - " ocentric_coord = og2oc(record['lat_Y_North'], semi_major, semi_minor)\n", - " coord_360 = to_360(ocentric_coord)\n", - " return coord_360\n", - "\n", - "# function to convert long_X_East to ISIS_lon\n", - "def lon_ISIS_coord(record, semi_major, semi_minor):\n", - " ocentric_coord = og2oc(record['long_X_East'], semi_major, semi_minor)\n", - " coord_360 = to_360(ocentric_coord)\n", - " return coord_360\n", - "\n", - "def body_fix(record, semi_major, semi_minor, inverse=False):\n", - " \"\"\"\n", - " Parameters\n", - " ----------\n", - " record : ndarray\n", - " (n,3) where columns are x, y, height or lon, lat, alt\n", - " \"\"\"\n", - " \n", - " ecef = pyproj.Proj(proj='geocent', a=semi_major, b=semi_minor)\n", - " lla = pyproj.Proj(proj='latlon', a=semi_major, b=semi_minor)\n", - " \n", - " if inverse:\n", - " lon, lat, height = pyproj.transform(ecef, lla, record[0], record[1], record[2])\n", - " return lon, lat, height\n", - " else:\n", - " y, x, z = pyproj.transform(lla, ecef, record[0], record[1], record[2])\n", - " return y, x, z\n", + "def socet2isis(at_file, cub_file_path, cub_ipf_map, target_name, outpath=None):\n", + " # Setup the at_file, path to cubes, and control network out path\n", + " at_file = at_file\n", + " cnet_out = os.path.split(os.path.splitext(at_file)[0])[1]\n", + " cub_path = cub_file_path\n", "\n", - "def ignore_toggle(record):\n", - " if record['stat'] == 0:\n", - " return True\n", + " if(outpath):\n", + " outpath = outpath\n", " else:\n", - " return False\n", - "\n", - "# TODO: Does isis cnet need a convariance matrix for sigmas? Even with a static matrix of 1,1,1,1 \n", - "def compute_sigma_covariance_matrix(lat, lon, rad, latsigma, lonsigma, radsigma, semimajor_axis):\n", - " \n", - " \"\"\"\n", - " Given geospatial coordinates, desired accuracy sigmas, and an equitorial radius, compute a 2x3\n", - " sigma covariange matrix.\n", - " \n", - " Parameters\n", - " ----------\n", - " lat : float\n", - " A point's latitude in degrees\n", - " \n", - " lon : float\n", - " A point's longitude in degrees\n", - " \n", - " rad : float\n", - " The radius (z-value) of the point in meters\n", - " \n", - " latsigma : float\n", - " The desired latitude accuracy in meters (Default 10.0)\n", - " \n", - " lonsigma : float\n", - " The desired longitude accuracy in meters (Default 10.0)\n", - " \n", - " radsigma : float\n", - " The desired radius accuracy in meters (Defualt: 15.0)\n", - " \n", - " semimajor_axis : float\n", - " The semi-major or equitorial radius in meters (Default: 1737400.0 - Moon)\n", - " \n", - " Returns\n", - " -------\n", - " rectcov : ndarray\n", - " (2,3) covariance matrix\n", - " \"\"\"\n", - " \n", - " lat = math.radians(lat)\n", - " lon = math.radians(lon)\n", - " \n", - " # SetSphericalSigmasDistance\n", - " scaled_lat_sigma = latsigma / semimajor_axis\n", - "\n", - " # This is specific to each lon.\n", - " scaled_lon_sigma = lonsigma * math.cos(lat) / semimajor_axis\n", - " \n", - " # SetSphericalSigmas\n", - " cov = np.eye(3,3)\n", - " cov[0,0] = scaled_lat_sigma ** 2\n", - " cov[1,1] = scaled_lon_sigma ** 2\n", - " cov[2,2] = radsigma ** 2\n", - " \n", - " # Approximate the Jacobian\n", - " j = np.zeros((3,3))\n", - " cosphi = math.cos(lat)\n", - " sinphi = math.sin(lat)\n", - " coslambda = math.cos(lon)\n", - " sinlambda = math.sin(lon)\n", - " rcosphi = rad * cosphi\n", - " rsinphi = rad * sinphi\n", - " j[0,0] = -rsinphi * coslambda\n", - " j[0,1] = -rcosphi * sinlambda\n", - " j[0,2] = cosphi * coslambda\n", - " j[1,0] = -rsinphi * sinlambda\n", - " j[1,1] = rcosphi * coslambda\n", - " j[1,2] = cosphi * sinlambda\n", - " j[2,0] = rcosphi\n", - " j[2,1] = 0.\n", - " j[2,2] = sinphi\n", - " mat = j.dot(cov)\n", - " mat = mat.dot(j.T)\n", - " rectcov = np.zeros((2,3))\n", - " rectcov[0,0] = mat[0,0]\n", - " rectcov[0,1] = mat[0,1]\n", - " rectcov[0,2] = mat[0,2]\n", - " rectcov[1,0] = mat[1,1]\n", - " rectcov[1,1] = mat[1,2]\n", - " rectcov[1,2] = mat[2,2]\n", - " \n", - " return np.array(rectcov)\n", - "# return np.array([[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]])\n", - "\n", + " outpath = os.path.split(at_file)[0]\n", "\n", - "def compute_cov_matrix(record, semimajor_axis):\n", - " cov_matrix = compute_sigma_covariance_matrix(record['lat_Y_North'], record['long_X_East'], record['ht'], record['sig0'], record['sig1'], record['sig2'], semimajor_axis)\n", - " return cov_matrix.ravel().tolist()\n", - "\n", - "# applys transformations to columns\n", - "def apply_two_isis_transformations(atf_dict, df):\n", - " prj_file = os.path.join(atf_dict['PATH'], atf_dict['PROJECT'])\n", - " \n", - " eRadius, pRadius = get_axis(prj_file)\n", - " \n", - " lla = np.array([[df['long_X_East']], [df['lat_Y_North']], [df['ht']]])\n", - " \n", - " ecef = body_fix(lla, semi_major = eRadius, semi_minor = pRadius, inverse=False)\n", - " \n", - " df['s.'], df['l.'], df['image_index'] = (zip(*df.apply(line_sample_size, path = atf_dict['PATH'], axis=1)))\n", - " df['known'] = df.apply(known, axis=1)\n", - " df['long_X_East'] = ecef[0][0]\n", - " df['lat_Y_North'] = ecef[1][0]\n", - " df['ht'] = ecef[2][0] \n", - " df['aprioriCovar'] = df.apply(compute_cov_matrix, semimajor_axis = eRadius, axis=1)\n", - "# df['ignore'] = df.apply(ignore_toggle, axis=1)\n", + " with open(cub_ipf_map) as cub_ipf_map:\n", + " reader = csv.reader(cub_ipf_map, delimiter = ',')\n", + " image_dict = dict([(row[0], row[1]) for row in reader])\n", " \n", - "def socet2isis(prj_file):\n", " # Read in and setup the atf dict of information\n", - " atf_dict = read_atf(prj_file)\n", + " atf_dict = read_atf(at_file)\n", " \n", " # Get the gpf and ipf files using atf dict\n", " gpf_file = os.path.join(atf_dict['PATH'], atf_dict['GP_FILE']);\n", @@ -326,66 +74,160 @@ " socet_df = ipf_df.merge(gpf_df, left_on='pt_id', right_on='point_id')\n", " \n", " # Apply the transformations\n", - "# apply_two_isis_transformations(atf_dict, socet_df)\n", + " apply_transformations(atf_dict, socet_df)\n", " \n", " # Define column remap for socet dataframe\n", - "# column_map = {'pt_id': 'id', 'l.': 'y', 's.': 'x',\n", - "# 'res_l': 'lineResidual', 'res_s': 'sampleResidual', 'known': 'Type',\n", - "# 'lat_Y_North': 'aprioriY', 'long_X_East': 'aprioriX', 'ht': 'aprioriZ',\n", - "# 'sig0': 'aprioriLatitudeSigma', 'sig1': 'aprioriLongitudeSigma', 'sig2': 'aprioriRadiusSigma',\n", - "# 'sig_l': 'linesigma', 'sig_s': 'samplesigma'}\n", + " column_map = {'pt_id': 'id', 'l.': 'y', 's.': 'x',\n", + " 'res_l': 'lineResidual', 'res_s': 'sampleResidual', 'known': 'Type',\n", + " 'lat_Y_North': 'aprioriY', 'long_X_East': 'aprioriX', 'ht': 'aprioriZ',\n", + " 'sig0': 'aprioriLatitudeSigma', 'sig1': 'aprioriLongitudeSigma', 'sig2': 'aprioriRadiusSigma',\n", + " 'sig_l': 'linesigma', 'sig_s': 'samplesigma'}\n", " \n", " # Rename the columns using the column remap above\n", - "# socet_df.rename(columns = column_map, inplace=True)\n", + " socet_df.rename(columns = column_map, inplace=True)\n", " \n", - " # Return the socet dataframe to be converted to a control net\n", - " return socet_df\n", + " serial_dict = serial_numbers(image_dict, cub_path)\n", "\n", - "# creates a dict of serial numbers with the cub being the key\n", - "def serial_numbers(images, path, extension):\n", - " serial_dict = dict()\n", - " \n", - " for image in images:\n", - " snum = sn.generate_serial_number(os.path.join(path, image + extension))\n", - " snum = snum.replace('Mars_Reconnaissance_Orbiter', 'MRO')\n", - " serial_dict[image] = snum\n", - " return serial_dict" + " # creates the control network\n", + " cn.to_isis(os.path.join(outpath, cnet_out + '.net'), socet_df, serial_dict, targetname = targetname)\n", + " return socet_df" ] }, { "cell_type": "code", - "execution_count": 86, + "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Setup stuffs for the cub information namely the path and extension\n", - "path = '/home/acpaquette/repos/plio/test_cubes'\n", + "cub_path = '/Volumes/Blueman/'\n", "targetname = 'Mars'\n", - "# Extension of your cub files\n", - "extension = '.8bit.cub'\n", + "cub_map = '/Users/adampaquette/repos/plio/plio/examples/SocetSet/cub_map2.csv'\n", "\n", "# Path to atf file\n", - "atf_file = ('/home/acpaquette/repos/plio/plio/examples/SocetSet/Relative.atf')\n", + "atf_file = ('/Users/adampaquette/repos/plio/plio/examples/SocetSet/Relative.atf')\n", "\n", - "socet_df = socet2isis(atf_file)\n", + "socet_df = socet2isis(atf_file, cub_path, cub_map, targetname)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def reverse_known(record):\n", + " \"\"\"\n", + " Converts the known field from a socet dataframe into the\n", + " isis point_type column\n", "\n", - "# images = pd.unique(socet_df['ipf_file'])\n", + " Parameters\n", + " ----------\n", + " record : object\n", + " Pandas series object\n", + "\n", + " Returns\n", + " -------\n", + " : str\n", + " String representation of a known field\n", + " \"\"\"\n", + " record_type = record['known']\n", + " if record_type == 0 or record_type == 2:\n", + " return 0\n", + "\n", + " elif record_type == 1 or record_type == 3 or record_type == 4:\n", + " return 3\n", + " \n", + "def lat_socet_coord(record, semi_major, semi_minor):\n", + " \"\"\"\n", + " Function to convert lat_Y_North to ISIS_lat\n", + "\n", + " Parameters\n", + " ----------\n", + " record : object\n", + " Pandas series object\n", + "\n", + " semi_major : float\n", + " Radius from the center of the body to the equater\n", + "\n", + " semi_minor : float\n", + " Radius from the pole to the center of mass\n", + "\n", + " Returns\n", + " -------\n", + " coord_360 : float\n", + " Converted latitude into ocentric space, and mapped\n", + " into 0 to 360\n", + " \"\"\"\n", + " ographic_coord = oc2og(record['lat_Y_North'], semi_major, semi_minor)\n", + " return ((ographic_coord + 180) % 360) - 180\n", + "\n", + "def lon_socet_coord(record, semi_major, semi_minor):\n", + " \"\"\"\n", + " Function to convert lat_Y_North to ISIS_lat\n", + "\n", + " Parameters\n", + " ----------\n", + " record : object\n", + " Pandas series object\n", + "\n", + " semi_major : float\n", + " Radius from the center of the body to the equater\n", + "\n", + " semi_minor : float\n", + " Radius from the pole to the center of mass\n", "\n", - "# serial_dict = serial_numbers(images, path, extension)\n", + " Returns\n", + " -------\n", + " coord_360 : float\n", + " Converted latitude into ocentric space, and mapped\n", + " into 0 to 360\n", + " \"\"\"\n", + " ographic_coord = oc2og(record['long_X_East'], semi_major, semi_minor)\n", + " return ((ographic_coord + 180) % 360) - 180\n", + "\n", + "def fix_sample_line(record, serial_dict, extension, cub_path):\n", + " # Cube location to load\n", + " cube = pvl.load(os.path.join(cub_path, serial_dict[record['serialnumber']] + extension))\n", + " line_size = find_in_dict(cube, 'Lines')\n", + " sample_size = find_in_dict(cube, 'Samples')\n", "\n", - "# creates the control network\n", - "# cn.to_isis('/home/acpaquette/repos/plio/plio/examples/SocetSet/cn.net', socet_df, serial_dict, targetname = targetname)" + " new_line = record['l.'] - (int(line_size)/2.0) - 1\n", + " new_sample = record['s.'] - (int(sample_size)/2.0) - 1\n", + " return new_line, new_sample\n", + "\n", + "def ignore_toggle(record):\n", + " if record['stat'] == True:\n", + " return 0\n", + " else:\n", + " return 1" ] }, { "cell_type": "code", - "execution_count": 116, - "metadata": {}, + "execution_count": null, + "metadata": { + "scrolled": false + }, "outputs": [], "source": [ - "return_df = cn.from_isis(\"/home/acpaquette/repos/plio/plio/examples/SocetSet/cn.net\")\n", + "return_df = cn.from_isis(\"/Users/adampaquette/repos/plio/plio/examples/SocetSet/Relative.net\")\n", + "\n", + "eRadius = 3.39619000000000e+006\n", + "pRadius = eRadius * (1 - 1.08339143554195e-001)\n", + "\n", + "adjusted_flag = False\n", + "\n", + "cub_path = '/Volumes/Blueman/'\n", + "extension = '.8bit.cub'\n", + "cub_list = ['D06_029601_1846_XN_04N224W', \n", + " 'F05_037684_1857_XN_05N224W']\n", + "\n", + "# \n", + "cub_dict = {i: i + extension for i in cub_list}\n", + "serial_dict = {sn.generate_serial_number(os.path.join(cub_path, i + extension)): i for i in cub_list}\n", "\n", "columns = []\n", "column_index = []\n", @@ -395,571 +237,82 @@ " column_index.append(i)\n", " columns.append(column)\n", "\n", - "return_df = return_df.iloc[:, column_index]" + "return_df = return_df.iloc[:, column_index]\n", + "\n", + "column_map = {'id': 'pt_id', 'line': 'l.', 'sample': 's.', \n", + " 'lineResidual': 'res_l', 'sampleResidual': 'res_s', 'type': 'known', \n", + " 'aprioriLatitudeSigma': 'sig0', 'aprioriLongitudeSigma': 'sig1', 'aprioriRadiusSigma': 'sig2', \n", + " 'linesigma': 'sig_l', 'samplesigma': 'sig_s', 'ignore': 'stat'}\n", + "\n", + "if adjusted_flag:\n", + " column_map['adjustedY'] = 'lat_Y_North'\n", + " column_map['adjustedX'] = 'long_X_East'\n", + " column_map['adjustedZ'] = 'ht'\n", + "else:\n", + " column_map['aprioriY'] = 'lat_Y_North'\n", + " column_map['aprioriX'] = 'long_X_East'\n", + " column_map['aprioriZ'] = 'ht'\n", + "\n", + "return_df.rename(columns = column_map, inplace=True)\n", + "\n", + "return_df['known'] = return_df.apply(reverse_known, axis = 1)\n", + "return_df['ipf_file'] = return_df['serialnumber'].apply(lambda serial_number: serial_dict[serial_number])\n", + "return_df['l.'], return_df['s.'] = zip(*return_df.apply(fix_sample_line, serial_dict = serial_dict, \n", + " extension = extension, \n", + " cub_path = cub_path, axis = 1))\n", + "\n", + "ecef = np.array([[return_df['long_X_East']], [return_df['lat_Y_North']], [return_df['ht']]])\n", + "lla = body_fix(ecef, semi_major = eRadius, semi_minor = pRadius, inverse=True)\n", + "return_df['long_X_East'], return_df['lat_Y_North'], return_df['ht'] = lla[0][0], lla[1][0], lla[2][0]\n", + "\n", + "return_df['lat_Y_North'] = return_df.apply(lat_socet_coord, semi_major=eRadius, semi_minor=pRadius, axis = 1)\n", + "return_df['long_X_East'] = return_df.apply(lon_socet_coord, semi_major=eRadius, semi_minor=pRadius, axis = 1)\n", + "\n", + "return_df['stat'] = return_df.apply(ignore_toggle, axis = 1)\n", + "return_df['val'] = return_df['stat']\n", + "\n", + "# Add dumby\n", + "x_dummy = lambda x: np.full(len(return_df), x)\n", + "\n", + "return_df['sig0'] = x_dummy(0)\n", + "return_df['sig1'] = x_dummy(0)\n", + "return_df['sig2'] = x_dummy(0)\n", + "\n", + "return_df['res0'] = x_dummy(0)\n", + "return_df['res1'] = x_dummy(0)\n", + "return_df['res2'] = x_dummy(0)\n", + "\n", + "return_df['fid_x'] = x_dummy(0)\n", + "return_df['fid_y'] = x_dummy(0)\n", + "\n", + "return_df['no_obs'] = x_dummy(1)\n", + "return_df['fid_val'] = x_dummy(0)" ] }, { "cell_type": "code", - "execution_count": 117, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "column_map = {'pt_id': 'id', 'l.': 'y', 's.': 'x',\n", - " 'res_l': 'lineResidual', 'res_s': 'sampleResidual', 'known': 'Type',\n", - " 'lat_Y_North': 'aprioriY', 'long_X_East': 'aprioriX', 'ht': 'aprioriZ',\n", - " 'sig0': 'aprioriLatitudeSigma', 'sig1': 'aprioriLongitudeSigma', 'sig2': 'aprioriRadiusSigma',\n", - " 'sig_l': 'linesigma', 'sig_s': 'samplesigma'}\n", - "\n", - "column_map = {k: v for v, k in column_map.items()}\n", - "return_df.rename(columns = column_map, inplace=True)\n", - "return_df.drop(['chooserName', 'datetime', 'referenceIndex', 'jigsawRejected', 'editLock', 'aprioriSurfPointSource', 'aprioriSurfPointSourceFile','aprioriRadiusSource', 'aprioriRadiusSourceFile'] , axis = 1, inplace=True)" + "return_df[['long_X_East', 'lat_Y_North', 'ht']].iloc[2]" ] }, { "cell_type": "code", - "execution_count": 129, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
lat_Y_Northlong_X_Eastht
0139525.2307493.390974e+064506.496945
1139525.2307493.390974e+064506.496945
2139489.2780453.390969e+064516.454802
3139489.2780453.390969e+064516.454802
4139823.4897973.390990e+064536.274914
5139823.4897973.390990e+064536.274914
6139772.7380043.390936e+064518.050219
7139772.7380043.390936e+064518.050219
8139575.9148153.390952e+063816.666542
9139575.9148153.390952e+063816.666542
10139614.7562963.390953e+063791.232717
11139614.7562963.390953e+063791.232717
12139912.0413743.390914e+063875.608660
13139912.0413743.390914e+063875.608660
14139909.4520333.390930e+063845.361327
15139909.4520333.390930e+063845.361327
16139669.8268493.391120e+063270.672620
17139669.8268493.391120e+063270.672620
18139694.5170173.391205e+063289.744506
19139694.5170173.391205e+063289.744506
20139968.7933383.391126e+063274.711397
21139968.7933383.391126e+063274.711397
22139979.2007803.391138e+063298.297228
23139979.2007803.391138e+063298.297228
24139688.0312173.391041e+064253.956077
25139688.0312173.391041e+064253.956077
26139686.9108233.391089e+064216.743792
27139686.9108233.391089e+064216.743792
28139786.2052843.390979e+063579.127600
29139786.2052843.390979e+063579.127600
30139785.0109973.391002e+063546.549796
31139785.0109973.391002e+063546.549796
\n", - "
" - ], - "text/plain": [ - " lat_Y_North long_X_East ht\n", - "0 139525.230749 3.390974e+06 4506.496945\n", - "1 139525.230749 3.390974e+06 4506.496945\n", - "2 139489.278045 3.390969e+06 4516.454802\n", - "3 139489.278045 3.390969e+06 4516.454802\n", - "4 139823.489797 3.390990e+06 4536.274914\n", - "5 139823.489797 3.390990e+06 4536.274914\n", - "6 139772.738004 3.390936e+06 4518.050219\n", - "7 139772.738004 3.390936e+06 4518.050219\n", - "8 139575.914815 3.390952e+06 3816.666542\n", - "9 139575.914815 3.390952e+06 3816.666542\n", - "10 139614.756296 3.390953e+06 3791.232717\n", - "11 139614.756296 3.390953e+06 3791.232717\n", - "12 139912.041374 3.390914e+06 3875.608660\n", - "13 139912.041374 3.390914e+06 3875.608660\n", - "14 139909.452033 3.390930e+06 3845.361327\n", - "15 139909.452033 3.390930e+06 3845.361327\n", - "16 139669.826849 3.391120e+06 3270.672620\n", - "17 139669.826849 3.391120e+06 3270.672620\n", - "18 139694.517017 3.391205e+06 3289.744506\n", - "19 139694.517017 3.391205e+06 3289.744506\n", - "20 139968.793338 3.391126e+06 3274.711397\n", - "21 139968.793338 3.391126e+06 3274.711397\n", - "22 139979.200780 3.391138e+06 3298.297228\n", - "23 139979.200780 3.391138e+06 3298.297228\n", - "24 139688.031217 3.391041e+06 4253.956077\n", - "25 139688.031217 3.391041e+06 4253.956077\n", - "26 139686.910823 3.391089e+06 4216.743792\n", - "27 139686.910823 3.391089e+06 4216.743792\n", - "28 139786.205284 3.390979e+06 3579.127600\n", - "29 139786.205284 3.390979e+06 3579.127600\n", - "30 139785.010997 3.391002e+06 3546.549796\n", - "31 139785.010997 3.391002e+06 3546.549796" - ] - }, - "execution_count": 129, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ - "return_df[['lat_Y_North', 'long_X_East', 'ht']]" + "return_df[['long_X_East', 'lat_Y_North', 'ht']].iloc[2]" ] }, { "cell_type": "code", - "execution_count": 128, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
lat_Y_Northlong_X_Eastht
00.0957082.356167-2342.889214
10.0957082.356167-2342.889214
20.0959202.355564-2349.638414
30.0959202.355564-2349.638414
40.0963392.361186-2314.316425
50.0963392.361186-2314.316425
60.0959542.360368-2370.502882
70.0959542.360368-2370.502882
80.0810582.357037-2363.989968
90.0810582.357037-2363.989968
100.0805182.357691-2360.922571
110.0805182.357691-2360.922571
120.0823112.362733-2388.123298
130.0823112.362733-2388.123298
140.0816682.362678-2371.973499
150.0816682.362678-2371.973499
160.0694582.358505-2193.309629
170.0694582.358505-2193.309629
180.0698612.358862-2106.769773
190.0698612.358862-2106.769773
200.0695432.363543-2174.971745
210.0695432.363543-2174.971745
220.0700442.363710-2162.103231
230.0700442.363710-2162.103231
240.0903422.358866-2269.610862
250.0903422.358866-2269.610862
260.0895502.358814-2222.328983
270.0895502.358814-2222.328983
280.0760122.360565-2328.281125
290.0760122.360565-2328.281125
300.0753202.360529-2305.362047
310.0753202.360529-2305.362047
\n", - "
" - ], - "text/plain": [ - " lat_Y_North long_X_East ht\n", - "0 0.095708 2.356167 -2342.889214\n", - "1 0.095708 2.356167 -2342.889214\n", - "2 0.095920 2.355564 -2349.638414\n", - "3 0.095920 2.355564 -2349.638414\n", - "4 0.096339 2.361186 -2314.316425\n", - "5 0.096339 2.361186 -2314.316425\n", - "6 0.095954 2.360368 -2370.502882\n", - "7 0.095954 2.360368 -2370.502882\n", - "8 0.081058 2.357037 -2363.989968\n", - "9 0.081058 2.357037 -2363.989968\n", - "10 0.080518 2.357691 -2360.922571\n", - "11 0.080518 2.357691 -2360.922571\n", - "12 0.082311 2.362733 -2388.123298\n", - "13 0.082311 2.362733 -2388.123298\n", - "14 0.081668 2.362678 -2371.973499\n", - "15 0.081668 2.362678 -2371.973499\n", - "16 0.069458 2.358505 -2193.309629\n", - "17 0.069458 2.358505 -2193.309629\n", - "18 0.069861 2.358862 -2106.769773\n", - "19 0.069861 2.358862 -2106.769773\n", - "20 0.069543 2.363543 -2174.971745\n", - "21 0.069543 2.363543 -2174.971745\n", - "22 0.070044 2.363710 -2162.103231\n", - "23 0.070044 2.363710 -2162.103231\n", - "24 0.090342 2.358866 -2269.610862\n", - "25 0.090342 2.358866 -2269.610862\n", - "26 0.089550 2.358814 -2222.328983\n", - "27 0.089550 2.358814 -2222.328983\n", - "28 0.076012 2.360565 -2328.281125\n", - "29 0.076012 2.360565 -2328.281125\n", - "30 0.075320 2.360529 -2305.362047\n", - "31 0.075320 2.360529 -2305.362047" - ] - }, - "execution_count": 128, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "socet_df[['lat_Y_North', 'long_X_East', 'ht']]" - ] + "outputs": [], + "source": [] }, { "cell_type": "code", @@ -985,7 +338,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.4" + "version": "3.6.3" } }, "nbformat": 4, diff --git a/plio/examples/SocetSet/cub_map.csv b/plio/examples/SocetSet/cub_map1.csv similarity index 100% rename from plio/examples/SocetSet/cub_map.csv rename to plio/examples/SocetSet/cub_map1.csv diff --git a/plio/examples/SocetSet/cub_map2.csv b/plio/examples/SocetSet/cub_map2.csv new file mode 100644 index 0000000000000000000000000000000000000000..a9401a26740e13f6e657c74be5d60bab898ec5c4 --- /dev/null +++ b/plio/examples/SocetSet/cub_map2.csv @@ -0,0 +1,2 @@ +D06_029601_1846_XN_04N224W,D06_029601_1846_XN_04N224W.8bit.cub +F05_037684_1857_XN_05N224W,F05_037684_1857_XN_05N224W.8bit.cub diff --git a/plio/io/io_bae.py b/plio/io/io_bae.py index 2810a8e19e91a8c9ad2472b2adb3d92783ff1faf..16b2a68186a00ea1b387833f356f7bd546d707db 100644 --- a/plio/io/io_bae.py +++ b/plio/io/io_bae.py @@ -1,6 +1,7 @@ -import json import re import os +import json +from collections import defaultdict from functools import singledispatch import numpy as np @@ -290,44 +291,48 @@ def read_atf(atf_file): project """ with open(atf_file) as f: - - files = [] - ipf = [] - sup = [] + # Extensions of files we want + files_ext = ['.prj', '.sup', '.ipf', '.gpf'] files_dict = [] + files = defaultdict(list) - # Grabs every PRJ, GPF, SUP, and IPF image from the ATF file for line in f: - if line[-4:-1] == 'prj' or line[-4:-1] == 'gpf' or line[-4:-1] == 'sup' or line[-4:-1] == 'ipf' or line[-4:-1] == 'atf': - files.append(line) + ext = os.path.splitext(line)[-1].strip() - files = np.array(files) + # Check is needed for split as all do not have a space + if ext in files_ext: - # Creates appropriate arrays for certain files in the right format - for file in files: - file = file.strip() - file = file.split(' ') + # If it is the .prj file, it strips the directory away and grabs file name + if ext == '.prj': + files[ext].append(line.strip().split(' ')[1].split('\\')[-1]) - # Grabs all the IPF files - if file[1].endswith('.ipf'): - ipf.append(file[1]) + # If the ext is in the list of files we care about, it addes to the dict + files[ext].append(line.strip().split(' ')[-1]) - # Grabs all the SUP files - if file[1].endswith('.sup'): - sup.append(file[1]) + else: - files_dict.append(file) + # Adds to the dict even if not in files we care about + files[ext.strip()].append(line) + + # Gets the base filepath + files['basepath'] = os.path.dirname(os.path.abspath(atf_file)) # Creates a dict out of file lists for GPF, PRJ, IPF, and ATF files_dict = (dict(files_dict)) # Sets the value of IMAGE_IPF to all IPF images - files_dict['IMAGE_IPF'] = ipf + files_dict['IMAGE_IPF'] = files['.ipf'] # Sets the value of IMAGE_SUP to all SUP images - files_dict['IMAGE_SUP'] = sup + files_dict['IMAGE_SUP'] = files['.sup'] + + # Sets value for GPF file + files_dict['GP_FILE'] = files['.gpf'][0] + + # Sets value for PRJ file + files_dict['PROJECT'] = files['.prj'][0] # Sets the value of PATH to the path of the ATF file - files_dict['PATH'] = os.path.dirname(os.path.abspath(atf_file)) + files_dict['PATH'] = files['basepath'] return files_dict diff --git a/plio/spatial/transformations.py b/plio/spatial/transformations.py index efc4ad609424a1c88e4ca778370b647249655e4c..7d4e7903f5352678a21e9646901a99b930e010bc 100644 --- a/plio/spatial/transformations.py +++ b/plio/spatial/transformations.py @@ -2,6 +2,8 @@ import os import math import pyproj +import numpy as np + import plio.io.isis_serial_number as sn def line_sample_size(record, path): @@ -217,7 +219,7 @@ def lon_ISIS_coord(record, semi_major, semi_minor): coord_360 = to_360(ocentric_coord) return coord_360 -def body_fix(record, semi_major, semi_minor): +def body_fix(record, semi_major, semi_minor, inverse = False): """ Transforms latitude, longitude, and height of a socet point into a body fixed point @@ -241,8 +243,19 @@ def body_fix(record, semi_major, semi_minor): """ 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 + + if inverse: + lon, lat, height = pyproj.transform(ecef, lla, record[0], record[1], record[2]) + return lon, lat, height + else: + y, x, z = pyproj.transform(lla, ecef, record[0], record[1], record[2]) + return y, x, z + +def stat_toggle(record): + if record['stat'] == 0: + return True + else: + return False def apply_transformations(atf_dict, df): """ @@ -259,15 +272,21 @@ def apply_transformations(atf_dict, df): Pandas dataframe object """ - prj_file = os.path.join(atf_dict['PATH'], atf_dict['PROJECT'].split('\\')[-1]) + prj_file = os.path.join(atf_dict['PATH'], atf_dict['PROJECT']) eRadius, pRadius = get_axis(prj_file) + lla = np.array([[df['long_X_East']], [df['lat_Y_North']], [df['ht']]]) + + ecef = body_fix(lla, semi_major = eRadius, semi_minor = pRadius, inverse=False) + 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)) + 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 serial_numbers(image_dict, path): """ @@ -290,3 +309,85 @@ def serial_numbers(image_dict, path): for key in image_dict: serial_dict[key] = sn.generate_serial_number(os.path.join(path, image_dict[key])) return serial_dict + +# 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] = scaled_lat_sigma ** 2 + cov[1,1] = 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()