Skip to content
Snippets Groups Projects
Commit c2587e91 authored by Tyler Thatcher's avatar Tyler Thatcher
Browse files

Added covariance matrix

parent 88b4acee
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
notebooks/Socet2ISIS.ipynb
+ 10
− 4
View file @ c2587e91
%% Cell type:code id: tags:
``` python
import os
import sys
from functools import singledispatch
import warnings
import pandas as pd
import numpy as np
import math
import pyproj
# Path to local plio if wanted
sys.path.insert(0, "/home/tthatcher/Desktop/Projects/Plio/plio")
sys.path.insert(0, "/path/to/plio")
from plio.examples import get_path
from collections import defaultdict
from plio.io.io_bae import read_gpf, read_ipf
import plio.io.io_controlnetwork as cn
import plio.io.isis_serial_number as sn
```
%% Cell type:code id: tags:
``` python
# Reads a .atf file and outputs all of the
# .ipf, .gpf, .sup, .prj, and path to locate the
# .apf file (should be the same as all others)
def read_atf(atf_file):
with open(atf_file) as f:
files = []
ipf = []
sup = []
files_dict = []
# 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)
files = np.array(files)
# Creates appropriate arrays for certain files in the right format
for file in files:
file = file.strip()
file = file.split(' ')
# Grabs all the IPF files
if file[1].endswith('.ipf'):
ipf.append(file[1])
# Grabs all the SUP files
if file[1].endswith('.sup'):
sup.append(file[1])
files_dict.append(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
# Sets the value of IMAGE_SUP to all SUP images
files_dict['IMAGE_SUP'] = sup
# Sets the value of PATH to the path of the ATF file
files_dict['PATH'] = os.path.dirname(os.path.abspath(atf_file))
return files_dict
# converts columns l. and s. to isis
def line_sample_size(record, path):
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):
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):
return num % 360
# ocentric to ographic latitudes
def oc2og(dlat, dMajorRadius, dMinorRadius):
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
# ographic to ocentric latitudes
def og2oc(dlat, dMajorRadius, dMinorRadius):
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):
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
# function to convert lat_Y_North to ISIS_lat
def lat_ISIS_coord(record, semi_major, semi_minor):
ocentric_coord = og2oc(record['lat_Y_North'], semi_major, semi_minor)
coord_360 = to_360(ocentric_coord)
return coord_360
# function to convert long_X_East to ISIS_lon
def lon_ISIS_coord(record, semi_major, semi_minor):
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, inverse=False):
"""
Parameters
----------
record : ndarray
(n,3) where columns are x, y, height or lon, lat, alt
"""
ecef = pyproj.Proj(proj='geocent', a=semi_major, b=semi_minor)
lla = pyproj.Proj(proj='latlon', a=semi_major, b=semi_minor)
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 ignore_toggle(record):
if record['stat'] == 0:
return True
else:
return False
# TODO: Does isis cnet need a convariance matrix for sigmas? Even with a static matrix of 1,1,1,1 no output
# 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)
coslambda = math.cos(lon)
sinlambda = math.sin(lon)
rcosphi = rad * cosphi
rsinphi = rad * sinphi
j[0,0] = -rsinphi * coslambda
j[0,1] = -rcosphi * sinlambda
j[0,2] = cosphi * coslambda
j[1,0] = -rsinphi * sinlambda
j[1,1] = rcosphi * coslambda
j[1,2] = cosphi * sinlambda
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
# return np.array([[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]])
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()
# applys transformations to columns
def apply_transformations(atf_dict, df):
prj_file = os.path.join(atf_dict['PATH'], atf_dict['PROJECT'].split('\\')[-1])
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['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['ignore'] = df.apply(ignore_toggle, axis=1)
def socet2isis(prj_file):
# Read in and setup the atf dict of information
atf_dict = read_atf(prj_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',
'sig_l': 'linesigma', 'sig_s': 'samplesigma'}
# Rename the columns using the column remap above
socet_df.rename(columns = column_remap, inplace=True)
# Return the socet dataframe to be converted to a control net
return socet_df
# creates a dict of serial numbers with the cub being the key
def serial_numbers(images, path, extension):
serial_dict = dict()
for image in images:
snum = sn.generate_serial_number(os.path.join(path, image + extension))
snum = snum.replace('Mars_Reconnaissance_Orbiter', 'MRO')
serial_dict[image] = snum
return serial_dict
```
%% Cell type:code id: tags:
``` python
# path = '/home/tthatcher/Desktop/Projects/plio_imgs/quest_imgs/'
# extension = '.lev1.cub'
# atf_file = ('/home/tthatcher/Desktop/Projects/plio_imgs/quest_imgs/CTX_Athabasca_Middle_step0.atf')
# Setup stuffs for the cub information namely the path and extension
path = '/home/tthatcher/Desktop/Projects/plio_imgs/correct_imgs/'
path = '/path/where/cub/files/are/'
# Extension of your cub files
extension = '.8bit.cub'
atf_file = ('/where/y')
# Path to atf file
atf_file = ('/path/to/atf/file')
socet_df = socet2isis(atf_file)
images = pd.unique(socet_df['ipf_file'])
serial_dict = serial_numbers(images, path, extension)
# creates the control network
cn.to_isis('/path/you/want/the/cnet/to/be/in/cn.net', socet_df, serial_dict)
# cn.to_isis('/home/tthatcher/Desktop/Projects/plio_imgs/quest_imgs/cn.net', socet_df, serial_dict)
```
%% Cell type:code id: tags:
``` python
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment