diff --git a/notebooks/Socet2ISIS.ipynb b/notebooks/Socet2ISIS.ipynb index f9ff8e2dd784c0c36eca6a96d8133d847e531657..a0e578889f8f89a93fed0d8d1f25978b9cc08d4c 100644 --- a/notebooks/Socet2ISIS.ipynb +++ b/notebooks/Socet2ISIS.ipynb @@ -16,10 +16,9 @@ "import math\n", "import pyproj\n", "\n", - "# sys.path.insert(0, \"/home/tthatcher/Desktop/Projects/Plio/plio\")\n", - "\n", "from plio.examples import get_path\n", "from plio.io.io_bae import read_gpf, read_ipf\n", + "from collections import defaultdict\n", "import plio.io.io_controlnetwork as cn\n", "import plio.io.isis_serial_number as sn" ] @@ -35,45 +34,50 @@ "# .apf file (should be the same as all others) \n", "def read_atf(atf_file):\n", " with open(atf_file) as f:\n", - "\n", - " files = []\n", - " ipf = []\n", - " sup = []\n", - " files_dict = []\n", " \n", - " # Grabs every PRJ, GPF, SUP, and IPF image from the ATF file\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", - " 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':\n", - " files.append(line)\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", - " files = np.array(files)\n", + " # Gets the base filepath\n", + " files['basepath'] = os.path.dirname(os.path.abspath(atf_file))\n", " \n", - " # Creates appropriate arrays for certain files in the right format\n", - " for file in files:\n", - " file = file.strip()\n", - " file = file.split(' ')\n", - "\n", - " # Grabs all the IPF files\n", - " if file[1].endswith('.ipf'):\n", - " ipf.append(file[1])\n", - "\n", - " # Grabs all the SUP files\n", - " if file[1].endswith('.sup'):\n", - " sup.append(file[1])\n", - "\n", - " files_dict.append(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'] = ipf\n", + " files_dict['IMAGE_IPF'] = files['.ipf']\n", " \n", " # Sets the value of IMAGE_SUP to all SUP images\n", - " files_dict['IMAGE_SUP'] = sup\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'] = os.path.dirname(os.path.abspath(atf_file))\n", + " files_dict['PATH'] = files['basepath']\n", " \n", " return files_dict\n", "\n", @@ -163,25 +167,126 @@ " coord_360 = to_360(ocentric_coord)\n", " return coord_360\n", "\n", - "def body_fix(record, semi_major, semi_minor):\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", - " lon, lat, height = pyproj.transform(lla, ecef, record['long_X_East'], record['lat_Y_North'], record['ht'])\n", - " return lon, lat, height\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", + "\n", + "def ignore_toggle(record):\n", + " if record['stat'] == 0:\n", + " return True\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", + " Parameters\n", + " ----------\n", + " lat : float\n", + " A point's latitude in degrees\n", + " lon : float\n", + " A point's longitude in degrees\n", + " rad : float\n", + " The radius (z-value) of the point in meters\n", + " latsigma : float\n", + " The desired latitude accuracy in meters (Default 10.0)\n", + " lonsigma : float\n", + " The desired longitude accuracy in meters (Default 10.0)\n", + " radsigma : float\n", + " The desired radius accuracy in meters (Defualt: 15.0)\n", + " semimajor_axis : float\n", + " The semi-major or equitorial radius in meters (Default: 1737400.0 - Moon)\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 rectcov\n", + "# return np.array([[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]])\n", + "\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_transformations(atf_dict, df):\n", - " prj_file = os.path.join(atf_dict['PATH'], atf_dict['PROJECT'].split('\\\\')[-1])\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['lat_Y_North'] = df.apply(lat_ISIS_coord, semi_major = eRadius, semi_minor = pRadius, axis=1)\n", - " df['long_X_East'] = df.apply(lon_ISIS_coord, semi_major = eRadius, semi_minor = pRadius, axis=1)\n", - " df['long_X_East'], df['lat_Y_North'], df['ht'] = zip(*df.apply(body_fix, semi_major = eRadius, semi_minor = pRadius, axis = 1))\n", - " \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", + " \n", "def socet2isis(prj_file):\n", " # Read in and setup the atf dict of information\n", " atf_dict = read_atf(prj_file)\n", @@ -213,9 +318,10 @@ " \n", " # Define column remap for socet dataframe\n", " column_remap = {'l.': 'y', 's.': 'x',\n", - " 'res_l': 'LineResidual', 'res_s': 'SampleResidual', 'known': 'Type',\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", + " '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_remap, inplace=True)\n", @@ -224,44 +330,42 @@ " return socet_df\n", "\n", "# creates a dict of serial numbers with the cub being the key\n", - "def serial_numbers(image_dict, path):\n", + "def serial_numbers(images, path, extension):\n", " serial_dict = dict()\n", - "\n", - " for key in image_dict:\n", - " snum = sn.generate_serial_number(os.path.join(path, image_dict[key]))\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[key] = snum\n", - " return serial_dict" + " serial_dict[image] = snum\n", + " return serial_dict\n", + "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { - "scrolled": false + "scrolled": true }, "outputs": [], "source": [ "# Setup stuffs for the cub information namely the path and extension\n", - "path = '/Volumes/Blueman/'\n", - "atf_file = get_path('CTX_Athabasca_Middle_step0.atf')\n", + "path = '/path/where/cub/files/are/'\n", + "\n", + "# Extension of your cub files\n", + "extension = '.something.cub'\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", + "# Path to atf file\n", + "atf_file = ('/path/to/atf/file')\n", "\n", "socet_df = socet2isis(atf_file)\n", "\n", "images = pd.unique(socet_df['ipf_file'])\n", "\n", - "serial_dict = serial_numbers(image_dict, path)\n", + "serial_dict = serial_numbers(images, path, extension)\n", "\n", "# creates the control network\n", - "cn.to_isis('/Volumes/Blueman/banana.net', socet_df, serial_dict)" + "cn.to_isis('/path/you/want/the/cnet/to/be/in/cn.net', socet_df, serial_dict)" ] }, { @@ -288,7 +392,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.3" + "version": "3.6.4" } }, "nbformat": 4, diff --git a/plio/io/io_controlnetwork.py b/plio/io/io_controlnetwork.py index 87b877a604648018d03daedd4f8695890c28ea40..a52d58097723627958bf7bbc54c847e57c8f60e0 100644 --- a/plio/io/io_controlnetwork.py +++ b/plio/io/io_controlnetwork.py @@ -1,6 +1,7 @@ from time import gmtime, strftime import pandas as pd +import numpy as np import pvl from plio.io import ControlNetFileV0002_pb2 as cnf @@ -189,19 +190,19 @@ class IsisStore(object): if version == 2: point_attrs = [i for i in cnf._CONTROLPOINTFILEENTRYV0002.fields_by_name if i != 'measures'] measure_attrs = [i for i in cnf._CONTROLPOINTFILEENTRYV0002_MEASURE.fields_by_name] - + cols = point_attrs + measure_attrs cp = cnf.ControlPointFileEntryV0002() self._handle.seek(header_start_byte) pbuf_header = cnf.ControlNetFileHeaderV0002() pbuf_header.ParseFromString(self._handle.read(header_bytes)) - + self._handle.seek(point_start_byte) cp = cnf.ControlPointFileEntryV0002() pts = [] for s in pbuf_header.pointMessageSizes: - cp.ParseFromString(self._handle.read(s)) + cp.ParseFromString(self._handle.read(s)) pt = [getattr(cp, i) for i in point_attrs if i != 'measures'] for measure in cp.measures: @@ -267,24 +268,24 @@ class IsisStore(object): # As per protobuf docs for assigning to a repeated field. if attr == 'aprioriCovar': arr = g.iloc[0]['aprioriCovar'] - point_spec.aprioriCovar.extend(arr.ravel().tolist()) + if isinstance(arr, np.ndarray): + arr = arr.ravel().tolist() + + point_spec.aprioriCovar.extend(arr) else: setattr(point_spec, attr, attrtype(g.iloc[0][attr])) - point_spec.type = 2 # Hardcoded to free + point_spec.type = 2 # Hardcoded to free this is bad # The reference index should always be the image with the lowest index point_spec.referenceIndex = 0 - # A single extend call is cheaper than many add calls to pack points measure_iterable = [] - for node_id, m in g.iterrows(): measure_spec = point_spec.Measure() # For all of the attributes, set if they are an dict accessible attr of the obj. for attr, attrtype in self.measure_attrs: if attr in g.columns: setattr(measure_spec, attr, attrtype(m[attr])) - measure_spec.serialnumber = serials[m.image_index] measure_spec.sample = m.x measure_spec.line = m.y @@ -298,7 +299,6 @@ class IsisStore(object): point_message = point_spec.SerializeToString() point_sizes.append(point_spec.ByteSize()) point_messages.append(point_message) - return point_messages, point_sizes def create_buffer_header(self, networkid, targetname,