Added notebook to generate EIS segments (91b38b3c) · Commits · aflab / astrogeology / Knoten

archive/EIS_spice.ipynb

0 → 100755

+520 −0

Original line number	Original line	Diff line number	Diff line
			%% Cell type:code id:dc00d6d1-7b58-4b36-a52c-e0dfaa57ad01 tags:

			``` python
			# %matplotlib inline
			import os
			isisroot_dir = "/Path/to/isis_data"
			os.environ["ISISDATA"] = isisroot_dir
			os.environ["ISISROOT"] = "/Path/to/isis_root"
			import csv
			from glob import glob
			import json
			import logging
			import math
			from pathlib import Path
			import subprocess

			import numpy as np

			import plotly
			import plotly.graph_objs as go
			import matplotlib
			matplotlib.use('qtagg')
			import matplotlib.pyplot as plt
			plotly.offline.init_notebook_mode(connected=True)

			import spiceypy as spice

			import ale
			from ale.drivers.clipper_drivers import ClipperEISWACPBIsisLabelNaifSpiceDriver
			from ale import isd_generate

			from kalasiris import isis

			from knoten.utils import reproject
			```

			%% Cell type:code id:ab9988b4-9a66-4a16-89be-03f99a7154f1 tags:

			``` python
			# Define various planet radii as spheres
			mars_radii = np.array([3396.2, 3396.2, 3396.2])
			europa_radii = np.array([1560.8, 1560.8, 1560.8])

			# Get the EIS cube
			label_path = "/Path/to/eis/workarea"
			clipper_file = os.path.join(label_path, "EIS000XXX_2032116T234928_0000C35F-WAC-PUSHB-IMG_RAW_V1.cub")


			# Build kernel set
			kernels = ['base/kernels/pck/pck00009.tpc',
			'base/kernels/spk/de430.bsp',
			'base/kernels/spk/mar097.bsp',
			'base/kernels/lsk/naif0012.tls']
			kernels = [os.path.join(isisdata_root, kernel) for kernel in kernels]

			eis_path = "/Path/to/eis_data/kernels"
			eis_kernels = ["sclk/europaclipper_00000.tsc",
			"EUROPA_SCLKSCET.00001.tsc",
			"ik/clipper_eis_v06.ti",
			"fk/clipper_v09.tf",
			"iak/eis_data.iak"]
			eis_kernels = [os.path.join(eis_path, kernel) for kernel in eis_kernels]
			eis_kernels = eis_kernels + kernels
			eis_kernels
			```

			%% Cell type:code id:f7adad53-2e5a-40c2-8932-ae3fd3d19de6 tags:

			``` python
			def normalize(a):
			a = np.array(a)
			return a/np.linalg.norm(a)

			# Degree coords to center the approach on
			latitude = -5
			longitude = 75

			scale = mars_radii[0]/europa_radii[0]

			# Height in km(?)
			height = 50 * scale

			# Direction of flight use "east" and "west"
			# "east" is flying west toward east
			# "west" is flygin east toward west
			direction = "east"

			# spacecraft speed in km/s
			speed = 4.5 * scale

			# +/- time from closest approach.
			# Total time = start + (time_offset * 2)
			time_offset = 500
			# time_offset = time_offset * 2

			# Radius definition to use
			radii = mars_radii

			# Compute and normalize nadir vector at closest approach
			y, x, z = reproject([longitude, latitude, height], radii[0], radii[2], "latlong", "geocent")
			spacecraft_point = np.array([x, y, z])
			nadir_vector = normalize(spacecraft_point)

			# Create a parallel trajectory that is perpendicular with nadir and parallel
			# to the x, y axis
			normal = np.array([0, 0, 1])
			perpendicular_nadir = np.cross(nadir_vector, normal)

			# compute fly_distance in km given speed and expected time from
			# closest approach
			fly_distance = speed * time_offset

			# Compute the ending X,Y,Z given by extending the trajectory vector by the fly_distance
			if direction == "west":
			start_point = spacecraft_point + (perpendicular_nadir * fly_distance)
			end_point = spacecraft_point - (perpendicular_nadir * fly_distance)
			elif direction == "east":
			start_point = spacecraft_point - (perpendicular_nadir * fly_distance)
			end_point = spacecraft_point + (perpendicular_nadir * fly_distance)
			else:
			print("Direction " + direction + " not defined")

			# Generate X, Y, Z some number of points from start and end
			num_points = 10000
			x_line_space = np.linspace(start_point[0], end_point[0], num_points)
			y_line_space = np.linspace(start_point[1], end_point[1], num_points)
			z_line_space = np.linspace(start_point[2], end_point[2], num_points)


			# Zip up X, Y, Z line spaces and compute velocities
			positions = np.array(list(zip(x_line_space, y_line_space, z_line_space)))
			velocities = np.array([positions[1] - positions[0]] * num_points)
			velocities = np.array([normalize(velocity) * speed for velocity in velocities])
			times = np.linspace(-time_offset, time_offset, num_points)
			times[int((num_points/2) - 1)], positions[int((num_points/2) - 1)], velocities[int((num_points/2) - 1)]
			```

			%% Cell type:code id:79b66a6c-f6f2-4a9c-9e7f-5b69485fcacc tags:

			``` python
			# Cell examines computed trajectory and look vector against the radii/wireframe
			%matplotlib qt

			def WireframeSphere(centre=[0.,0.,0.], radius=[1., 1., 1.],
			n_meridians=20, n_circles_latitude=None):
			"""
			Create the arrays of values to plot the wireframe of a sphere.

			Parameters
			----------
			centre: array like
			A point, defined as an iterable of three numerical values.
			radius: number
			The radius of the sphere.
			n_meridians: int
			The number of meridians to display (circles that pass on both poles).
			n_circles_latitude: int
			The number of horizontal circles (akin to the Equator) to display.
			Notice this includes one for each pole, and defaults to 4 or half
			of the n_meridians if the latter is larger.

			Returns
			-------
			sphere_x, sphere_y, sphere_z: arrays
			The arrays with the coordinates of the points to make the wireframe.
			Their shape is (n_meridians, n_circles_latitude).

			Examples
			--------
			>>> fig = plt.figure()
			>>> ax = fig.gca(projection='3d')
			>>> ax.set_aspect("equal")
			>>> sphere = ax.plot_wireframe(*WireframeSphere(), color="r", alpha=0.5)
			>>> fig.show()

			>>> fig = plt.figure()
			>>> ax = fig.gca(projection='3d')
			>>> ax.set_aspect("equal")
			>>> frame_xs, frame_ys, frame_zs = WireframeSphere()
			>>> sphere = ax.plot_wireframe(frame_xs, frame_ys, frame_zs, color="r", alpha=0.5)
			>>> fig.show()
			"""
			if n_circles_latitude is None:
			n_circles_latitude = max(n_meridians/2, 4)
			u, v = np.mgrid[0:2np.pi:n_meridians1j, 0:np.pi:n_circles_latitude*1j]
			sphere_x = centre[0] + radius[0] * np.cos(u) * np.sin(v)
			sphere_y = centre[1] + radius[1] * np.sin(u) * np.sin(v)
			sphere_z = centre[2] + radius[2] * np.cos(v)
			return sphere_x, sphere_y, sphere_z

			fig = plt.figure()
			ax = fig.add_subplot(projection = '3d')
			ax.set_aspect("auto")

			x_line_space = positions[:, 0]
			y_line_space = positions[:, 1]
			z_line_space = positions[:, 2]

			look_vectors = np.array([normalize(position) * 1000 for position in positions])

			# Plot sphere as x, y, z
			frame_xs, frame_ys, frame_zs = WireframeSphere(radius=mars_radii)
			sphere = ax.plot_wireframe(frame_xs, frame_ys, frame_zs, color="r", alpha=0.5)
			ax.plot3D(x_line_space, y_line_space, z_line_space, 'gray', marker='o')
			for i, lv in enumerate(look_vectors):
			ax.plot3D([lv[0], positions[i][0]], [lv[1], positions[i][1]], [lv[2], positions[i][2]], 'orange', marker='o')

			# Plot radii and center point in ECEF
			ax.plot3D(mars_radii[0], 0, 0, 'blue', marker='o')
			ax.plot3D(0, mars_radii[1], 0, 'orange', marker='o')
			ax.plot3D(0, 0, mars_radii[2], 'green', marker='o')
			ax.plot3D(0, 0, 0, 'black', marker='o')

			fig.show()
			```

			%% Cell type:code id:e7dedf0d-46da-4590-aedd-9192bc52e8dc tags:

			``` python
			spk_file = "/Path/to/custom_clipper_path.spk"
			if Path(spk_file).is_file():
			os.remove(spk_file)

			driver = ClipperEISWACPBIsisLabelNaifSpiceDriver(clipper_file, props={"nadir": True, "kernels": eis_kernels})
			with driver as active_driver:
			ephem_start = active_driver.ephemeris_start_time
			ephemeris_time = [time + ephem_start for time in times]

			handle = spice.spkopn(spk_file, "SPK", 512)
			states = pyspiceql.concatStates(positions, velocities)

			spice.spkw13(handle,
			-159,
			499,
			"IAU_MARS",
			ephemeris_time[0],
			ephemeris_time[-1],
			"Custom Trajectory",
			3,
			len(ephemeris_time),
			states,
			ephemeris_time)

			spice.spkcls(handle)
			eis_kernels.append(spk_file)
			```

			%% Cell type:code id:0d1609e8-db68-4212-921b-62b200272afa tags:

			``` python
			def interp_3d(time, xyz_array, times):
			f1 = np.interp(time, times, xyz_array[:, 0])
			f2 = np.interp(time, times, xyz_array[:, 1])
			f3 = np.interp(time, times, xyz_array[:, 2])

			return np.array([f1, f2, f3])

			def compute_segments_with_trajectory(ts,
			radii,
			IFOV,
			block_size,
			smear_fraction,
			tset,
			hmax,
			positions,
			velocities,
			times):
			segment_data = []

			# Setup initial conditions
			position = interp_3d(ts, positions, times)
			velocity = interp_3d(ts, velocities, times)
			r = np.linalg.norm(position)
			h = r - radii

			while (h > hmax):
			ts = ts + 0.5
			position = interp_3d(ts, positions, times)
			velocity = interp_3d(ts, velocities, times)
			r = np.linalg.norm(position)
			h = r - radii

			while h < hmax:
			direction = math.copysign(1, ts)
			position = interp_3d(ts, positions, times)
			velocity = interp_3d(ts, velocities, times)
			r = np.linalg.norm(position)
			h = r - radii
			pixwid = IFOV * h
			v = np.linalg.norm(velocity)
			unitR = normalize(position)
			unitC = normalize(np.cross(position, velocity))
			unitI = normalize(np.cross(unitC, unitR))
			vg = np.dot(velocity, unitI) * (radii/r)
			lts = pixwid / vg
			ltm = lts * (1 + direction * smear_fraction)
			Nb = 0
			smear = smear_fraction

			while smear <= smear_fraction:
			Nb = Nb + 1
			te = ts + ((Nb * block_size) - 1) * ltm
			position = interp_3d(te, positions, times)
			velocity = interp_3d(te, velocities, times)
			r = np.linalg.norm(position)
			v = np.linalg.norm(velocity)
			h = r - radii
			pixwid = IFOV * h
			unitR = normalize(position)
			unitC = normalize(np.cross(position, velocity))
			unitI = normalize(np.cross(unitC, unitR))
			vg = np.dot(velocity, unitI) * (radii/r)

			lte = pixwid / vg
			smear = abs(lte - ltm) / ltm

			Nb = max(1, Nb - 1)
			te = ts + ((Nb * block_size) - 1) * ltm
			data_set = [ts, te, Nb]
			segment_data.append(data_set)
			ts = te + tset

			return segment_data

			def compute_segments_const_velocity(ts, radii, speed, IFOV, block_size, smear_fraction, tset, height, hmax):
			segment_data = []

			v = speed
			s = v * ts
			r0 = radii + height
			r = math.sqrt(r02 + s2)
			h = r - radii

			while (h > hmax):
			ts = ts + 0.5
			s = v * ts
			r = math.sqrt(r02 + s2)
			h = r - radii

			while h < hmax:
			direction = math.copysign(1, ts)
			s = v * ts
			r = math.sqrt(r02 + s2)
			h = r - radii
			pixwid = IFOV * h
			vg = v * (radii/r0) / (1 + math.pow((s/r0), 2))
			lts = pixwid / vg
			ltm = lts * (1 + direction * smear_fraction)
			Nb = 0
			smear = smear_fraction
			while smear <= smear_fraction:
			Nb = Nb + 1
			te = ts + ((Nb * block_size) - 1) * ltm
			s = v * te
			r = math.sqrt(r02 + s2)
			h = r - radii
			pixwid = IFOV * h
			vg = v * (radii/r0) / (1 + math.pow((s/r0), 2))
			lte = pixwid / vg
			smear = abs(lte - ltm) / ltm

			Nb = max(1, Nb - 1)
			te = ts + ((Nb * block_size) - 1) * ltm
			data_set = [ts, te, Nb]
			segment_data.append(data_set)
			ts = te + tset
			return segment_data

			R = mars_radii[0]
			IFOV = 2.18e-4
			block_size = 256
			f = 0.05
			tset = 0.3
			hmax = 1000 * scale

			segments_from_trajectory = compute_segments_with_trajectory(-time_offset, R, IFOV, block_size, f, tset, hmax, positions, velocities, times)
			segments_from_const_velocity = compute_segments_const_velocity(-time_offset, R, speed, IFOV, block_size, f, tset, height, hmax)
			```

			%% Cell type:code id:30aa43ca-490d-4689-bd9f-d659ad636c65 tags:

			``` python
			# Define the fore, nadir, and aft sensor lines
			detector_offsets = [0, 1023, 2047]
			```

			%% Cell type:code id:0a492226-c1d6-47c1-88bd-e8f20d9b836e tags:

			``` python
			path = "/Path/to/eis_segments/"
			spk_file = eis_kernels[-1]
			ik_file = eis_kernels[3]
			sclk_file = eis_kernels[0]
			iak_file = eis_kernels[4]
			pck_file = eis_kernels[5]
			shape_file = os.path.join(isisroot_dir, "base/dems/molaMarsPlanetaryRadius0005.cub")
			for detector_offset in detector_offsets:
			i = 0
			driver = ClipperEISWACPBIsisLabelNaifSpiceDriver(clipper_file, props={"nadir": True, "kernels": eis_kernels})
			with driver as active_driver:
			ephem_start = active_driver.ephemeris_start_time

			for segment_data in segments_from_trajectory:
			image_path = os.path.join(path, f"eis_segment_{detector_offset}_{i:02}.cub")
			isis.makecube(to=image_path,
			pixels="VALUE",
			value=0.0,
			samples=4096,
			lines=int(segment_data[-1] * block_size),
			bands=1)
			isis.copylabel(from_=image_path,
			source=clipper_file,
			instrument=False,
			bandbin=True,
			kernels=True,
			mapping=False,
			radiometry=False,
			polygon=False,
			camstats=False)
			try:
			isis.editlab(from_=image_path, options="DELG", grpname="AlphaCube")
			except:
			pass
			isis.editlab(from_=image_path, options="ADDG", grpname="Instrument")

			# Most of the following PVL keyword edits are hard coded and may require changes
			# to the "value" argument to produce the expected/correct label
			isis.editlab(from_=image_path,
			options="SETKEY",
			grpname="Instrument",
			keyword="SpacecraftName",
			value="Clipper")
			isis.editlab(from_=image_path,
			options="SETKEY",
			grpname="Instrument",
			keyword="InstrumentId",
			value="WAC-PUSHBROOM")
			isis.editlab(from_=image_path,
			options="SETKEY",
			grpname="Instrument",
			keyword="TargetName",
			value="Mars")
			start_time = spice.et2utc(ephem_start + segment_data[0], 'C', 14)
			isis.editlab(from_=image_path,
			options="SETKEY",
			grpname="Instrument",
			keyword="StartTime",
			value=start_time)
			exposure_duration = segment_data[1] - segment_data[0]
			isis.editlab(from_=image_path,
			options="SETKEY",
			grpname="Instrument",
			keyword="ExposureDuration",
			value=exposure_duration, units="s")
			isis.editlab(from_=image_path,
			options="SETKEY",
			grpname="Instrument",
			keyword="DetectorOffset",
			value=detector_offset)

			# Create CSV for time table and write using csv2table
			csv_file_name = "temp.csv"
			with open(csv_file_name, 'w', newline='\n') as csvfile:
			fieldnames = ["EphemerisTime", "ExposureTime", "LineStart"]
			writer = csv.DictWriter(csvfile, fieldnames=fieldnames)

			writer.writeheader()
			writer.writerow({"EphemerisTime": segment_data[0] + ephem_start,
			"ExposureTime": exposure_duration/(segment_data[-1] * block_size),
			"LineStart": 1})
			isis.csv2table(csv=csv_file_name,
			tableName="LineScanTimes",
			to=image_path,
			coltypes="(Double, Double, Integer)")
			try:
			isis.spiceinit(from_=image_path,
			spk=spk_file,
			CKNADIR=True,
			CKRECON=False,
			ik=ik_file,
			sclk=sclk_file,
			iak=iak_file,
			pck=pck_file,
			shape="USER",
			model=shape_file)
			except subprocess.CalledProcessError as err:
			print('Had an ISIS error:')
			print(' '.join(err.cmd))
			print(err.stdout)
			print(err.stderr)
			raise err

			data_image_path = os.path.join(path, f"eis_segment_{detector_offset}_{i:02}.data.cub")
			pixel_data = os.path.join(path, "molaCenterShade.resamp.cub")
			if i <= 19:
			pixel_data = os.path.join(path, "molaRightShade.resamp.cub")
			elif i >= 60:
			pixel_data = os.path.join(path, "molaLeftShade.resamp.cub")
			isis.map2cam(from_=pixel_data,
			match=image_path,
			to=data_image_path)
			i += 1
			```

			%% Cell type:code id:30dd955e-5c28-4bee-b5df-98c2f062a676 tags:

			``` python
			for detector_offset in detector_offsets:
			files = glob(path, os.path.append(f"eis_segment_{detector_offset}_*.data.cub"))
			files.sort()

			for file in files:
			out_json = os.path.splitext(file)[0] + ".json"
			isd_generate.file_to_isd(file, out_json, log_level=logging.INFO, kernels = eis_kernels, only_naif_spice=True, nadir=True)
			```

			%% Cell type:code id:08086add-a50f-41ff-aa2c-ccb918bb97ea tags:

			``` python
			```