diff --git a/notebooks/Socet2ISIS.ipynb b/notebooks/Socet2ISIS.ipynb index 1d7681a240afaf6e83d53f2321d4150482ce5fd5..345a9422612be775c2b2e575ef8021661f7eb666 100644 --- a/notebooks/Socet2ISIS.ipynb +++ b/notebooks/Socet2ISIS.ipynb @@ -16,7 +16,8 @@ "import math\n", "import pyproj\n", "\n", - "# sys.path.insert(0, \"/home/tthatcher/Desktop/Projects/Plio/plio\")\n", + "# Path to local plio if wanted\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", @@ -26,7 +27,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -163,12 +164,107 @@ " 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 no output\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", @@ -176,12 +272,18 @@ " \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 +315,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", @@ -231,42 +334,37 @@ " 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" + " return serial_dict\n", + "\n" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 6, "metadata": { - "scrolled": false + "scrolled": true }, - "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": [ + "# path = '/home/tthatcher/Desktop/Projects/plio_imgs/quest_imgs/'\n", + "# extension = '.lev1.cub'\n", + "# atf_file = ('/home/tthatcher/Desktop/Projects/plio_imgs/quest_imgs/CTX_Athabasca_Middle_step0.atf')\n", + "\n", "# 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", + "path = '/home/tthatcher/Desktop/Projects/plio_imgs/correct_imgs/'\n", + "extension = '.8bit.cub'\n", + "atf_file = ('/where/y')\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", "\n", "# creates the control network\n", - "cn.to_isis('/Volumes/Blueman/cn.net', socet_df, serial_dict)" + "cn.to_isis('/path/you/want/the/cnet/to/be/in/cn.net', socet_df, serial_dict)\n", + "\n", + "# cn.to_isis('/home/tthatcher/Desktop/Projects/plio_imgs/quest_imgs/cn.net', socet_df, serial_dict)" ] }, { @@ -293,7 +391,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.3" + "version": "3.6.4" } }, "nbformat": 4,