From 65cd7563a17c43249ead250b7d516ff8ccff0d24 Mon Sep 17 00:00:00 2001 From: acpaquette Date: Mon, 25 Nov 2024 17:18:43 -0700 Subject: [PATCH] Added notebook to generate EIS segments Scrub paths --- archive/EIS_spice.ipynb | 547 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 547 insertions(+) create mode 100755 archive/EIS_spice.ipynb diff --git a/archive/EIS_spice.ipynb b/archive/EIS_spice.ipynb new file mode 100755 index 0000000..484992f --- /dev/null +++ b/archive/EIS_spice.ipynb @@ -0,0 +1,547 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "dc00d6d1-7b58-4b36-a52c-e0dfaa57ad01", + "metadata": {}, + "outputs": [], + "source": [ + "# %matplotlib inline\n", + "import os\n", + "isisroot_dir = \"/Path/to/isis_data\"\n", + "os.environ[\"ISISDATA\"] = isisroot_dir\n", + "os.environ[\"ISISROOT\"] = \"/Path/to/isis_root\"\n", + "import csv\n", + "from glob import glob\n", + "import json\n", + "import logging\n", + "import math\n", + "from pathlib import Path\n", + "import subprocess\n", + "\n", + "import numpy as np\n", + "\n", + "import plotly\n", + "import plotly.graph_objs as go\n", + "import matplotlib\n", + "matplotlib.use('qtagg')\n", + "import matplotlib.pyplot as plt\n", + "plotly.offline.init_notebook_mode(connected=True)\n", + "\n", + "import spiceypy as spice\n", + "\n", + "import ale\n", + "from ale.drivers.clipper_drivers import ClipperEISWACPBIsisLabelNaifSpiceDriver\n", + "from ale import isd_generate\n", + "\n", + "from kalasiris import isis\n", + "\n", + "from knoten.utils import reproject" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ab9988b4-9a66-4a16-89be-03f99a7154f1", + "metadata": {}, + "outputs": [], + "source": [ + "# Define various planet radii as spheres\n", + "mars_radii = np.array([3396.2, 3396.2, 3396.2])\n", + "europa_radii = np.array([1560.8, 1560.8, 1560.8])\n", + "\n", + "# Get the EIS cube\n", + "label_path = \"/Path/to/eis/workarea\"\n", + "clipper_file = os.path.join(label_path, \"EIS000XXX_2032116T234928_0000C35F-WAC-PUSHB-IMG_RAW_V1.cub\")\n", + "\n", + "\n", + "# Build kernel set\n", + "kernels = ['base/kernels/pck/pck00009.tpc',\n", + " 'base/kernels/spk/de430.bsp',\n", + " 'base/kernels/spk/mar097.bsp',\n", + " 'base/kernels/lsk/naif0012.tls']\n", + "kernels = [os.path.join(isisdata_root, kernel) for kernel in kernels]\n", + "\n", + "eis_path = \"/Path/to/eis_data/kernels\"\n", + "eis_kernels = [\"sclk/europaclipper_00000.tsc\",\n", + " \"EUROPA_SCLKSCET.00001.tsc\",\n", + " \"ik/clipper_eis_v06.ti\",\n", + " \"fk/clipper_v09.tf\",\n", + " \"iak/eis_data.iak\"]\n", + "eis_kernels = [os.path.join(eis_path, kernel) for kernel in eis_kernels]\n", + "eis_kernels = eis_kernels + kernels\n", + "eis_kernels" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f7adad53-2e5a-40c2-8932-ae3fd3d19de6", + "metadata": {}, + "outputs": [], + "source": [ + "def normalize(a):\n", + " a = np.array(a)\n", + " return a/np.linalg.norm(a)\n", + "\n", + "# Degree coords to center the approach on\n", + "latitude = -5\n", + "longitude = 75\n", + "\n", + "scale = mars_radii[0]/europa_radii[0]\n", + "\n", + "# Height in km(?)\n", + "height = 50 * scale\n", + "\n", + "# Direction of flight use \"east\" and \"west\"\n", + "# \"east\" is flying west toward east\n", + "# \"west\" is flygin east toward west\n", + "direction = \"east\"\n", + "\n", + "# spacecraft speed in km/s\n", + "speed = 4.5 * scale\n", + "\n", + "# +/- time from closest approach.\n", + "# Total time = start + (time_offset * 2)\n", + "time_offset = 500\n", + "# time_offset = time_offset * 2\n", + "\n", + "# Radius definition to use\n", + "radii = mars_radii\n", + "\n", + "# Compute and normalize nadir vector at closest approach\n", + "y, x, z = reproject([longitude, latitude, height], radii[0], radii[2], \"latlong\", \"geocent\")\n", + "spacecraft_point = np.array([x, y, z])\n", + "nadir_vector = normalize(spacecraft_point)\n", + "\n", + "# Create a parallel trajectory that is perpendicular with nadir and parallel\n", + "# to the x, y axis\n", + "normal = np.array([0, 0, 1])\n", + "perpendicular_nadir = np.cross(nadir_vector, normal)\n", + "\n", + "# compute fly_distance in km given speed and expected time from\n", + "# closest approach\n", + "fly_distance = speed * time_offset\n", + "\n", + "# Compute the ending X,Y,Z given by extending the trajectory vector by the fly_distance\n", + "if direction == \"west\":\n", + " start_point = spacecraft_point + (perpendicular_nadir * fly_distance)\n", + " end_point = spacecraft_point - (perpendicular_nadir * fly_distance)\n", + "elif direction == \"east\":\n", + " start_point = spacecraft_point - (perpendicular_nadir * fly_distance)\n", + " end_point = spacecraft_point + (perpendicular_nadir * fly_distance)\n", + "else:\n", + " print(\"Direction \" + direction + \" not defined\")\n", + "\n", + "# Generate X, Y, Z some number of points from start and end\n", + "num_points = 10000\n", + "x_line_space = np.linspace(start_point[0], end_point[0], num_points)\n", + "y_line_space = np.linspace(start_point[1], end_point[1], num_points)\n", + "z_line_space = np.linspace(start_point[2], end_point[2], num_points)\n", + "\n", + "\n", + "# Zip up X, Y, Z line spaces and compute velocities\n", + "positions = np.array(list(zip(x_line_space, y_line_space, z_line_space)))\n", + "velocities = np.array([positions[1] - positions[0]] * num_points)\n", + "velocities = np.array([normalize(velocity) * speed for velocity in velocities])\n", + "times = np.linspace(-time_offset, time_offset, num_points)\n", + "times[int((num_points/2) - 1)], positions[int((num_points/2) - 1)], velocities[int((num_points/2) - 1)]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "79b66a6c-f6f2-4a9c-9e7f-5b69485fcacc", + "metadata": {}, + "outputs": [], + "source": [ + "# Cell examines computed trajectory and look vector against the radii/wireframe\n", + "%matplotlib qt\n", + "\n", + "def WireframeSphere(centre=[0.,0.,0.], radius=[1., 1., 1.],\n", + " n_meridians=20, n_circles_latitude=None):\n", + " \"\"\"\n", + " Create the arrays of values to plot the wireframe of a sphere.\n", + "\n", + " Parameters\n", + " ----------\n", + " centre: array like\n", + " A point, defined as an iterable of three numerical values.\n", + " radius: number\n", + " The radius of the sphere.\n", + " n_meridians: int\n", + " The number of meridians to display (circles that pass on both poles).\n", + " n_circles_latitude: int\n", + " The number of horizontal circles (akin to the Equator) to display.\n", + " Notice this includes one for each pole, and defaults to 4 or half\n", + " of the *n_meridians* if the latter is larger.\n", + "\n", + " Returns\n", + " -------\n", + " sphere_x, sphere_y, sphere_z: arrays\n", + " The arrays with the coordinates of the points to make the wireframe.\n", + " Their shape is (n_meridians, n_circles_latitude).\n", + "\n", + " Examples\n", + " --------\n", + " >>> fig = plt.figure()\n", + " >>> ax = fig.gca(projection='3d')\n", + " >>> ax.set_aspect(\"equal\")\n", + " >>> sphere = ax.plot_wireframe(*WireframeSphere(), color=\"r\", alpha=0.5)\n", + " >>> fig.show()\n", + "\n", + " >>> fig = plt.figure()\n", + " >>> ax = fig.gca(projection='3d')\n", + " >>> ax.set_aspect(\"equal\")\n", + " >>> frame_xs, frame_ys, frame_zs = WireframeSphere()\n", + " >>> sphere = ax.plot_wireframe(frame_xs, frame_ys, frame_zs, color=\"r\", alpha=0.5)\n", + " >>> fig.show()\n", + " \"\"\"\n", + " if n_circles_latitude is None:\n", + " n_circles_latitude = max(n_meridians/2, 4)\n", + " u, v = np.mgrid[0:2*np.pi:n_meridians*1j, 0:np.pi:n_circles_latitude*1j]\n", + " sphere_x = centre[0] + radius[0] * np.cos(u) * np.sin(v)\n", + " sphere_y = centre[1] + radius[1] * np.sin(u) * np.sin(v)\n", + " sphere_z = centre[2] + radius[2] * np.cos(v)\n", + " return sphere_x, sphere_y, sphere_z\n", + "\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot(projection = '3d')\n", + "ax.set_aspect(\"auto\")\n", + "\n", + "x_line_space = positions[:, 0]\n", + "y_line_space = positions[:, 1]\n", + "z_line_space = positions[:, 2]\n", + "\n", + "look_vectors = np.array([normalize(position) * 1000 for position in positions])\n", + "\n", + "# Plot sphere as x, y, z\n", + "frame_xs, frame_ys, frame_zs = WireframeSphere(radius=mars_radii)\n", + "sphere = ax.plot_wireframe(frame_xs, frame_ys, frame_zs, color=\"r\", alpha=0.5)\n", + "ax.plot3D(x_line_space, y_line_space, z_line_space, 'gray', marker='o')\n", + "for i, lv in enumerate(look_vectors):\n", + " ax.plot3D([lv[0], positions[i][0]], [lv[1], positions[i][1]], [lv[2], positions[i][2]], 'orange', marker='o')\n", + "\n", + "# Plot radii and center point in ECEF\n", + "ax.plot3D(mars_radii[0], 0, 0, 'blue', marker='o')\n", + "ax.plot3D(0, mars_radii[1], 0, 'orange', marker='o')\n", + "ax.plot3D(0, 0, mars_radii[2], 'green', marker='o')\n", + "ax.plot3D(0, 0, 0, 'black', marker='o')\n", + "\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e7dedf0d-46da-4590-aedd-9192bc52e8dc", + "metadata": {}, + "outputs": [], + "source": [ + "spk_file = \"/Path/to/custom_clipper_path.spk\"\n", + "if Path(spk_file).is_file():\n", + " os.remove(spk_file)\n", + "\n", + "driver = ClipperEISWACPBIsisLabelNaifSpiceDriver(clipper_file, props={\"nadir\": True, \"kernels\": eis_kernels})\n", + "with driver as active_driver:\n", + " ephem_start = active_driver.ephemeris_start_time\n", + "ephemeris_time = [time + ephem_start for time in times]\n", + "\n", + "handle = spice.spkopn(spk_file, \"SPK\", 512)\n", + "states = pyspiceql.concatStates(positions, velocities)\n", + "\n", + "spice.spkw13(handle,\n", + " -159,\n", + " 499,\n", + " \"IAU_MARS\",\n", + " ephemeris_time[0],\n", + " ephemeris_time[-1],\n", + " \"Custom Trajectory\",\n", + " 3,\n", + " len(ephemeris_time),\n", + " states,\n", + " ephemeris_time)\n", + "\n", + "spice.spkcls(handle)\n", + "eis_kernels.append(spk_file)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d1609e8-db68-4212-921b-62b200272afa", + "metadata": {}, + "outputs": [], + "source": [ + "def interp_3d(time, xyz_array, times):\n", + " f1 = np.interp(time, times, xyz_array[:, 0])\n", + " f2 = np.interp(time, times, xyz_array[:, 1])\n", + " f3 = np.interp(time, times, xyz_array[:, 2])\n", + "\n", + " return np.array([f1, f2, f3])\n", + "\n", + "def compute_segments_with_trajectory(ts,\n", + " radii, \n", + " IFOV, \n", + " block_size, \n", + " smear_fraction, \n", + " tset, \n", + " hmax, \n", + " positions, \n", + " velocities, \n", + " times):\n", + " segment_data = []\n", + " \n", + " # Setup initial conditions\n", + " position = interp_3d(ts, positions, times)\n", + " velocity = interp_3d(ts, velocities, times)\n", + " r = np.linalg.norm(position)\n", + " h = r - radii\n", + " \n", + " while (h > hmax):\n", + " ts = ts + 0.5\n", + " position = interp_3d(ts, positions, times)\n", + " velocity = interp_3d(ts, velocities, times)\n", + " r = np.linalg.norm(position)\n", + " h = r - radii\n", + " \n", + " while h < hmax:\n", + " direction = math.copysign(1, ts)\n", + " position = interp_3d(ts, positions, times)\n", + " velocity = interp_3d(ts, velocities, times)\n", + " r = np.linalg.norm(position)\n", + " h = r - radii\n", + " pixwid = IFOV * h\n", + " v = np.linalg.norm(velocity)\n", + " unitR = normalize(position)\n", + " unitC = normalize(np.cross(position, velocity))\n", + " unitI = normalize(np.cross(unitC, unitR))\n", + " vg = np.dot(velocity, unitI) * (radii/r)\n", + " lts = pixwid / vg\n", + " ltm = lts * (1 + direction * smear_fraction)\n", + " Nb = 0\n", + " smear = smear_fraction\n", + " \n", + " while smear <= smear_fraction:\n", + " Nb = Nb + 1\n", + " te = ts + ((Nb * block_size) - 1) * ltm\n", + " position = interp_3d(te, positions, times)\n", + " velocity = interp_3d(te, velocities, times)\n", + " r = np.linalg.norm(position)\n", + " v = np.linalg.norm(velocity)\n", + " h = r - radii\n", + " pixwid = IFOV * h\n", + " unitR = normalize(position)\n", + " unitC = normalize(np.cross(position, velocity))\n", + " unitI = normalize(np.cross(unitC, unitR))\n", + " vg = np.dot(velocity, unitI) * (radii/r)\n", + " \n", + " lte = pixwid / vg\n", + " smear = abs(lte - ltm) / ltm\n", + " \n", + " Nb = max(1, Nb - 1) \n", + " te = ts + ((Nb * block_size) - 1) * ltm\n", + " data_set = [ts, te, Nb]\n", + " segment_data.append(data_set)\n", + " ts = te + tset\n", + "\n", + " return segment_data\n", + "\n", + "def compute_segments_const_velocity(ts, radii, speed, IFOV, block_size, smear_fraction, tset, height, hmax):\n", + " segment_data = []\n", + " \n", + " v = speed\n", + " s = v * ts\n", + " r0 = radii + height\n", + " r = math.sqrt(r0**2 + s**2)\n", + " h = r - radii\n", + " \n", + " while (h > hmax):\n", + " ts = ts + 0.5\n", + " s = v * ts\n", + " r = math.sqrt(r0**2 + s**2)\n", + " h = r - radii\n", + " \n", + " while h < hmax:\n", + " direction = math.copysign(1, ts)\n", + " s = v * ts\n", + " r = math.sqrt(r0**2 + s**2)\n", + " h = r - radii\n", + " pixwid = IFOV * h\n", + " vg = v * (radii/r0) / (1 + math.pow((s/r0), 2))\n", + " lts = pixwid / vg\n", + " ltm = lts * (1 + direction * smear_fraction)\n", + " Nb = 0\n", + " smear = smear_fraction\n", + " while smear <= smear_fraction:\n", + " Nb = Nb + 1\n", + " te = ts + ((Nb * block_size) - 1) * ltm\n", + " s = v * te\n", + " r = math.sqrt(r0**2 + s**2)\n", + " h = r - radii\n", + " pixwid = IFOV * h\n", + " vg = v * (radii/r0) / (1 + math.pow((s/r0), 2))\n", + " lte = pixwid / vg\n", + " smear = abs(lte - ltm) / ltm\n", + "\n", + " Nb = max(1, Nb - 1)\n", + " te = ts + ((Nb * block_size) - 1) * ltm\n", + " data_set = [ts, te, Nb]\n", + " segment_data.append(data_set)\n", + " ts = te + tset\n", + " return segment_data\n", + "\n", + "R = mars_radii[0]\n", + "IFOV = 2.18e-4\n", + "block_size = 256\n", + "f = 0.05\n", + "tset = 0.3\n", + "hmax = 1000 * scale\n", + "\n", + "segments_from_trajectory = compute_segments_with_trajectory(-time_offset, R, IFOV, block_size, f, tset, hmax, positions, velocities, times)\n", + "segments_from_const_velocity = compute_segments_const_velocity(-time_offset, R, speed, IFOV, block_size, f, tset, height, hmax)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30aa43ca-490d-4689-bd9f-d659ad636c65", + "metadata": {}, + "outputs": [], + "source": [ + "# Define the fore, nadir, and aft sensor lines\n", + "detector_offsets = [0, 1023, 2047]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0a492226-c1d6-47c1-88bd-e8f20d9b836e", + "metadata": {}, + "outputs": [], + "source": [ + "path = \"/Path/to/eis_segments/\"\n", + "for detector_offset in detector_offsets:\n", + " i = 0\n", + " driver = ClipperEISWACPBIsisLabelNaifSpiceDriver(clipper_file, props={\"nadir\": True, \"kernels\": eis_kernels})\n", + " with driver as active_driver:\n", + " ephem_start = active_driver.ephemeris_start_time\n", + " \n", + " for segment_data in segments_from_trajectory:\n", + " image_path = os.path.join(path, f\"eis_segment_{detector_offset}_{i:02}.cub\")\n", + " isis.makecube(to=image_path, pixels=\"VALUE\", value=0.0, samples=4096, lines=int(segment_data[-1] * block_size), bands=1)\n", + " isis.copylabel(from_=image_path, \n", + " source=clipper_file, \n", + " instrument=False,\n", + " bandbin=True,\n", + " kernels=True,\n", + " mapping=False,\n", + " radiometry=False,\n", + " polygon=False,\n", + " camstats=False)\n", + " try:\n", + " isis.editlab(from_=image_path, options=\"DELG\", grpname=\"AlphaCube\")\n", + " except:\n", + " pass\n", + " isis.editlab(from_=image_path, options=\"ADDG\", grpname=\"Instrument\")\n", + " \n", + " # Most of the following PVL keyword edits are hard coded and may require changes \n", + " # to the \"value\" argument to produce the expected/correct label\n", + " isis.editlab(from_=image_path, options=\"SETKEY\", grpname=\"Instrument\", keyword=\"SpacecraftName\", value=\"Clipper\")\n", + " isis.editlab(from_=image_path, options=\"SETKEY\", grpname=\"Instrument\", keyword=\"InstrumentId\", value=\"WAC-PUSHBROOM\")\n", + " isis.editlab(from_=image_path, options=\"SETKEY\", grpname=\"Instrument\", keyword=\"TargetName\", value=\"Mars\")\n", + " start_time = spice.et2utc(ephem_start + segment_data[0], 'C', 14)\n", + " isis.editlab(from_=image_path, options=\"SETKEY\", grpname=\"Instrument\", keyword=\"StartTime\", value=start_time)\n", + " exposure_duration = segment_data[1] - segment_data[0]\n", + " isis.editlab(from_=image_path, options=\"SETKEY\", grpname=\"Instrument\", keyword=\"ExposureDuration\", value=exposure_duration, units=\"s\")\n", + " isis.editlab(from_=image_path, options=\"SETKEY\", grpname=\"Instrument\", keyword=\"DetectorOffset\", value=detector_offset)\n", + " \n", + " # Create CSV for time table and write using csv2table\n", + " csv_file_name = \"temp.csv\"\n", + " with open(csv_file_name, 'w', newline='\\n') as csvfile:\n", + " fieldnames = [\"EphemerisTime\", \"ExposureTime\", \"LineStart\"]\n", + " writer = csv.DictWriter(csvfile, fieldnames=fieldnames)\n", + " \n", + " writer.writeheader()\n", + " writer.writerow({\"EphemerisTime\": segment_data[0] + ephem_start, \n", + " \"ExposureTime\": exposure_duration/(segment_data[-1] * block_size), \n", + " \"LineStart\": 1})\n", + " isis.csv2table(csv=csv_file_name, tableName=\"LineScanTimes\", to=image_path, coltypes=\"(Double, Double, Integer)\")\n", + " try:\n", + " isis.spiceinit(from_=image_path,\n", + " spk=\"/Users/acpaquette/custom_clipper_path.spk\",\n", + " CKNADIR=True,\n", + " CKRECON=False,\n", + " ik=\"/Users/acpaquette/eis_data/data/kernels/ik/clipper_eis_v05.ti\",\n", + " sclk=\"/Users/acpaquette/eis_data/data/kernels/sclk/europaclipper_00000.tsc\",\n", + " iak=\"/Users/acpaquette/eis_data/data/kernels/iak/eis_data.iak\",\n", + " pck=\"/Users/acpaquette/eis_data/data/kernels/pck/pck00009.tpc\",\n", + " shape=\"USER\",\n", + " model=\"/Volumes/workDrive/isis_data/base/dems/molaMarsPlanetaryRadius0005.cub\")\n", + " except subprocess.CalledProcessError as err:\n", + " print('Had an ISIS error:')\n", + " print(' '.join(err.cmd))\n", + " print(err.stdout)\n", + " print(err.stderr)\n", + " raise err\n", + " \n", + " data_image_path = os.path.join(path, f\"eis_segment_{detector_offset}_{i:02}.data.cub\")\n", + " pixel_data = \"/Volumes/workDrive/molaCenterShade.resamp.cub\"\n", + " if i <= 19:\n", + " pixel_data = \"/Volumes/workDrive/molaRightShade.resamp.cub\"\n", + " elif i >= 60:\n", + " pixel_data = \"/Volumes/workDrive/molaLeftShade.resamp.cub\"\n", + " isis.map2cam(from_=pixel_data,\n", + " match=image_path,\n", + " to=data_image_path)\n", + " i += 1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30dd955e-5c28-4bee-b5df-98c2f062a676", + "metadata": {}, + "outputs": [], + "source": [ + "for detector_offset in detector_offsets:\n", + " files = glob(path, os.path.append(f\"eis_segment_{detector_offset}_*.data.cub\"))\n", + " files.sort()\n", + " \n", + " for file in files:\n", + " out_json = os.path.splitext(file)[0] + \".json\"\n", + " isd_generate.file_to_isd(file, out_json, log_level=logging.INFO, kernels = eis_kernels, only_naif_spice=True, nadir=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "08086add-a50f-41ff-aa2c-ccb918bb97ea", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.20" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}