addressed PR feedback (2257afbe) · Commits · aflab / astrogeology / Knoten

examples/sensor_utils.ipynb

+22 −19

Original line number	Diff line number	Diff line
		%% Cell type:markdown id: tags:

		# Sensor Utils

		%% Cell type:code id: tags:

		``` python
		import os

		os.environ["ISISROOT"] = "/Users/astamile/ISIS3/build"
		os.environ["ISISDATA"] = "/Volumes/isis_data1/isis_data/"

		from csmapi import csmapi
		from knoten import csm, sensor_utils

		from knoten.shape import Ellipsoid
		from knoten.illuminator import Illuminator

		import ale
		import json
		```

		%% Cell type:markdown id: tags:

		## Create a usgscsm sensor model

		%% Cell type:code id: tags:

		``` python
		fileName = "data/N1573082850_1.cub"

		kernels = ale.util.generate_kernels_from_cube(fileName, expand=True)
		isd_string = ale.loads(fileName, props={'kernels': kernels})
		csm_isd = os.path.splitext(fileName)[0] + '.json'

		with open(csm_isd, 'w') as isd_file:
		isd_file.write(isd_string)
		```

		%% Output

		/Users/astamile/opt/anaconda3/envs/knoten/lib/python3.9/site-packages/osgeo/gdal.py:312: FutureWarning: Neither gdal.UseExceptions() nor gdal.DontUseExceptions() has been explicitly called. In GDAL 4.0, exceptions will be enabled by default.
		warnings.warn(

		%% Cell type:markdown id: tags:

		## Run Sensor Utils with usgscsm sensor model and image point

		%% Cell type:code id: tags:

		``` python
		camera = csm.create_csm(csm_isd)
		image_pt = csmapi.ImageCoord(511.5, 511.5)
		shape = Ellipsoid.from_csm_sensor(camera)
		illuminator = Illuminator()
		```

		%% Cell type:code id: tags:

		``` python
		phaseAngle = sensor_utils.phase_angle(image_pt, camera)
		phaseAngle = sensor_utils.phase_angle(image_pt, camera, shape, illuminator)

		phaseAngle
		```

		%% Output

		38.87212509629893
		38.87212509629895

		%% Cell type:code id: tags:

		``` python
		emissionAngle = sensor_utils.emission_angle(image_pt, camera)
		emissionAngle = sensor_utils.emission_angle(image_pt, camera, shape)

		emissionAngle
		```

		%% Output

		49.60309924896046
		49.60309924893989

		%% Cell type:code id: tags:

		``` python
		slantDistance = sensor_utils.slant_distance(image_pt, camera)
		slantDistance = sensor_utils.slant_distance(image_pt, camera, shape)

		slantDistance
		```

		%% Output

		2903512972.146144
		2903512972.146115

		%% Cell type:code id: tags:

		``` python
		targetCenterDistance = sensor_utils.target_center_distance(image_pt, camera)

		targetCenterDistance
		```

		%% Output

		2943536048.858226

		%% Cell type:code id: tags:

		``` python
		subSpacecraftPoint = sensor_utils.sub_spacecraft_point(image_pt, camera)

		subSpacecraftPoint
		```

		%% Output

		LatLon(lat=3.2229625890973583, lon=258.6197326526089)

		%% Cell type:code id: tags:

		``` python
		localRadius = sensor_utils.local_radius(image_pt, camera)
		localRadius = sensor_utils.local_radius(image_pt, camera, shape)

		localRadius
		```

		%% Output

		59096282.02424558
		59096282.024265066

		%% Cell type:code id: tags:

		``` python
		rightAscDec = sensor_utils.right_ascension_declination(image_pt, camera)

		rightAscDec
		```

		%% Output

		(79.34815579474038, -2.7790780986459485)

		%% Cell type:code id: tags:

		``` python
		lineResolution = sensor_utils.line_resolution(image_pt, camera)
		lineResolution = sensor_utils.line_resolution(image_pt, camera, shape)

		lineResolution
		```

		%% Output

		17397.960941876583
		17397.96094194587

		%% Cell type:code id: tags:

		``` python
		sampleResolution = sensor_utils.sample_resolution(image_pt, camera)
		sampleResolution = sensor_utils.sample_resolution(image_pt, camera, shape)

		sampleResolution
		```

		%% Output

		17397.933700407997
		17397.93370038153

		%% Cell type:code id: tags:

		``` python
		pixelResolution = sensor_utils.pixel_resolution(image_pt, camera)
		pixelResolution = sensor_utils.pixel_resolution(image_pt, camera, shape)

		pixelResolution
		```

		%% Output

		17397.94732114229
		17397.9473211637

knoten/csm.py

+3 −23

Original line number	Diff line number	Diff line
		@@ -9,7 +9,6 @@ import numpy as np
		import requests
		import scipy.stats
		from functools import singledispatch
		import spiceypy as spice

		from knoten.surface import EllipsoidDem

		@@ -42,28 +41,6 @@ def get_radii(camera):
		semi_minor = ellipsoid.getSemiMinorRadius()
		return semi_major, semi_minor

		def get_surface_normal(sensor, ground_pt):
		"""
		Given a sensor model and ground point, calculate the surface normal.

		Parameters
		----------
		sensor : object
		A CSM compliant sensor model object

		ground_pt: tuple
		The ground point as an (x, y, z) tuple

		Returns
		-------
		: tuple
		in the form (x, y, z)
		"""
		semi_major, semi_minor = get_radii(sensor)

		normal = spice.surfnm(semi_major, semi_major, semi_minor, np.array([ground_pt.x, ground_pt.y, ground_pt.z]))
		return utils.Point(normal[0], normal[1], normal[2])

		def create_camera(label, url='http://pfeffer.wr.usgs.gov/api/1.0/pds/'):
		"""
		Given an ALE supported label, create a CSM compliant ISD file. This func
		@@ -114,6 +91,9 @@ def get_state(sensor, image_pt):
		: dict
		Dictionary containing lookVec, sensorPos, sensorTime, and imagePoint
		"""
		if not isinstance(sensor, csmapi.RasterGM):
		raise TypeError("inputted sensor not a csm.RasterGM object")

		sensor_time = sensor.getImageTime(image_pt)
		locus = sensor.imageToRemoteImagingLocus(image_pt)
		sensor_position = sensor.getSensorPosition(image_pt)

knoten/illuminator.py

0 → 100644

+23 −0

Original line number	Diff line number	Diff line
		from knoten import csm, utils
		import spiceypy as spice
		import numpy as np

		class Illuminator:

		def __init__(self, position=None, velocity=None):
		self.position = position
		self.velocity = velocity

		def get_position_from_csm_sensor(self, sensor, ground_pt):
		sunEcefVec = sensor.getIlluminationDirection(ground_pt)
		self.position = utils.Point(ground_pt.x - sunEcefVec.x, ground_pt.y - sunEcefVec.y, ground_pt.z - sunEcefVec.z)
		return self.position

		# def get_position_from_time(self, sensor_state):
		# return 0

		# def get_velocity(self, sensor_state):
		# return 0



		No newline at end of file

knoten/sensor_utils.py

+111 −19

Original line number	Diff line number	Diff line
		@@ -2,7 +2,7 @@ from knoten import csm, utils, bundle
		from csmapi import csmapi
		import numpy as np

		def phase_angle(image_pt, sensor):
		def phase_angle(image_pt, sensor, shape, illuminator):
		"""
		Computes and returns phase angle, in degrees for a given image point.

		@@ -18,6 +18,12 @@ def phase_angle(image_pt, sensor):
		sensor : object
		A CSM compliant sensor model object

		shape : object
		A shape model object

		illuminator: object
		An illuminator object

		Returns
		-------
		: np.ndarray
		@@ -27,10 +33,9 @@ def phase_angle(image_pt, sensor):
		image_pt = csmapi.ImageCoord(*image_pt)

		sensor_state = csm.get_state(sensor, image_pt)
		ground_pt = csm.generate_ground_point(0.0, image_pt, sensor)
		ground_pt = shape.intersect_surface(sensor_state["sensorPos"], sensor_state["lookVec"])

		sunEcefVec = sensor.getIlluminationDirection(ground_pt)
		illum_pos = utils.Point(ground_pt.x - sunEcefVec.x, ground_pt.y - sunEcefVec.y, ground_pt.z - sunEcefVec.z)
		illum_pos = illuminator.get_position_from_csm_sensor(sensor, ground_pt)

		vec_a = utils.Point(sensor_state["sensorPos"].x - ground_pt.x,
		sensor_state["sensorPos"].y - ground_pt.y,
		@@ -42,7 +47,7 @@ def phase_angle(image_pt, sensor):

		return np.rad2deg(utils.sep_angle(vec_a, vec_b))

		def emission_angle(image_pt, sensor):
		def emission_angle(image_pt, sensor, shape):
		"""
		Computes and returns emission angle, in degrees, for a given image point.

		@@ -62,6 +67,9 @@ def emission_angle(image_pt, sensor):
		sensor : object
		A CSM compliant sensor model object

		shape : object
		A shape model object

		Returns
		-------
		: np.ndarray
		@@ -71,9 +79,9 @@ def emission_angle(image_pt, sensor):
		image_pt = csmapi.ImageCoord(*image_pt)

		sensor_state = csm.get_state(sensor, image_pt)
		ground_pt = csm.generate_ground_point(0.0, image_pt, sensor)
		ground_pt = shape.intersect_surface(sensor_state["sensorPos"], sensor_state["lookVec"])

		normal = csm.get_surface_normal(sensor, ground_pt)
		normal = shape.get_surface_normal(ground_pt)

		sensor_diff = utils.Point(sensor_state["sensorPos"].x - ground_pt.x,
		sensor_state["sensorPos"].y - ground_pt.y,
		@@ -81,7 +89,7 @@ def emission_angle(image_pt, sensor):

		return np.rad2deg(utils.sep_angle(normal, sensor_diff))

		def slant_distance(image_pt, sensor):
		def slant_distance(image_pt, sensor, shape):
		"""
		Computes the slant distance from the sensor to the ground point.

		@@ -93,6 +101,9 @@ def slant_distance(image_pt, sensor):
		sensor : object
		A CSM compliant sensor model object

		shape : object
		A shape model object

		Returns
		-------
		: np.ndarray
		@@ -102,7 +113,7 @@ def slant_distance(image_pt, sensor):
		image_pt = csmapi.ImageCoord(*image_pt)

		sensor_state = csm.get_state(sensor, image_pt)
		ground_pt = csm.generate_ground_point(0.0, image_pt, sensor)
		ground_pt = shape.intersect_surface(sensor_state["sensorPos"], sensor_state["lookVec"])

		return utils.distance(sensor_state["sensorPos"], ground_pt)

		@@ -154,7 +165,7 @@ def sub_spacecraft_point(image_pt, sensor):

		return utils.radians_to_degrees(lat_lon_rad)

		def local_radius(image_pt, sensor):
		def local_radius(image_pt, sensor, shape):
		"""
		Gets the local radius for a given image point.

		@@ -166,6 +177,9 @@ def local_radius(image_pt, sensor):
		sensor : object
		A CSM compliant sensor model object

		shape : object
		A shape model object

		Returns
		-------
		: np.ndarray
		@@ -174,7 +188,9 @@ def local_radius(image_pt, sensor):
		if not isinstance(image_pt, csmapi.ImageCoord):
		image_pt = csmapi.ImageCoord(*image_pt)

		ground_pt = csm.generate_ground_point(0.0, image_pt, sensor)
		sensor_state = csm.get_state(sensor, image_pt)
		ground_pt = shape.intersect_surface(sensor_state["sensorPos"], sensor_state["lookVec"])

		return utils.magnitude(ground_pt)

		def right_ascension_declination(image_pt, sensor):
		@@ -204,7 +220,8 @@ def right_ascension_declination(image_pt, sensor):

		return (ra_dec.lon, ra_dec.lat)

		def line_resolution(image_pt, sensor):

		def line_resolution(image_pt, sensor, shape):
		"""
		Compute the line resolution in meters per pixel for the current set point.

		@@ -224,6 +241,9 @@ def line_resolution(image_pt, sensor):
		sensor : object
		A CSM compliant sensor model object

		shape : object
		A shape model object

		Returns
		-------
		: np.ndarray
		@@ -232,14 +252,16 @@ def line_resolution(image_pt, sensor):
		if not isinstance(image_pt, csmapi.ImageCoord):
		image_pt = csmapi.ImageCoord(*image_pt)

		ground_pt = csm.generate_ground_point(0.0, image_pt, sensor)
		sensor_state = csm.get_state(sensor, image_pt)
		ground_pt = shape.intersect_surface(sensor_state["sensorPos"], sensor_state["lookVec"])

		image_partials = bundle.compute_image_partials(sensor, ground_pt)

		return np.sqrt(image_partials[0] * image_partials[0] +
		image_partials[2] * image_partials[2] +
		image_partials[4] * image_partials[4])

		def sample_resolution(image_pt, sensor):
		def sample_resolution(image_pt, sensor, shape):
		"""
		Compute the sample resolution in meters per pixel for the current set point.

		@@ -255,6 +277,9 @@ def sample_resolution(image_pt, sensor):
		sensor : object
		A CSM compliant sensor model object

		shape : object
		A shape model object

		Returns
		-------
		: np.ndarray
		@@ -263,14 +288,16 @@ def sample_resolution(image_pt, sensor):
		if not isinstance(image_pt, csmapi.ImageCoord):
		image_pt = csmapi.ImageCoord(*image_pt)

		ground_pt = csm.generate_ground_point(0.0, image_pt, sensor)
		sensor_state = csm.get_state(sensor, image_pt)
		ground_pt = shape.intersect_surface(sensor_state["sensorPos"], sensor_state["lookVec"])

		image_partials = bundle.compute_image_partials(sensor, ground_pt)

		return np.sqrt(image_partials[1] * image_partials[1] +
		image_partials[3] * image_partials[3] +
		image_partials[5] * image_partials[5])

		def pixel_resolution(image_pt, sensor):
		def pixel_resolution(image_pt, sensor, shape):
		"""
		Returns the pixel resolution at the current position in meters/pixel.

		@@ -282,15 +309,80 @@ def pixel_resolution(image_pt, sensor):
		sensor : object
		A CSM compliant sensor model object

		shape : object
		A shape model object

		Returns
		-------
		: np.ndarray
		pixel resolution in meters/pixel
		"""
		line_res = line_resolution(image_pt, sensor)
		samp_res = sample_resolution(image_pt, sensor)
		line_res = line_resolution(image_pt, sensor, shape)
		samp_res = sample_resolution(image_pt, sensor, shape)
		if (line_res < 0.0):
		return None
		if (samp_res < 0.0):
		return None
		return (line_res + samp_res) / 2.0


		# def sub_solar_point(image_pt, sensor, illuminator, shape):
		# sensor_state = csm.get_state(sensor, image_pt)

		# illum_pos = illuminator.get_position_from_time(sensor_state)
		# lookVec = utils.Point(-illum_pos.x, -illum_pos.y, -illum_pos.z)

		# pt = shape.intersect_surface(illum_pos, lookVec)
		# return utils.Point(pt.x, pt.y, pt.z)


		# def local_solar_time(image_pt, sensor, illuminator, shape):
		# if not isinstance(image_pt, csmapi.ImageCoord):
		# image_pt = csmapi.ImageCoord(*image_pt)

		# sensor_state = csm.get_state(sensor, image_pt)

		# sub_solar_pt = sub_solar_point(image_pt, sensor, illuminator, shape)

		# sub_solar_pt_degrees = utils.radians_to_degrees(utils.rect_to_spherical(sub_solar_pt))

		# ground_pt = shape.intersect_surface(sensor_state["sensorPos"], sensor_state["lookVec"])
		# spherical_pt = utils.rect_to_spherical(ground_pt)

		# spherical_pt_degrees = utils.radians_to_degrees(spherical_pt)

		# lst = spherical_pt_degrees.lon - sub_solar_pt_degrees.lon + 180.0
		# lst = lst / 15.0
		# if (lst < 0.0):
		# lst += 24.0
		# if (lst > 24.0):
		# lst -= 24.0
		# return lst

		# def solar_longitude(image_pt, sensor, illuminator, body):
		# sensor_state = csm.get_state(sensor, image_pt)

		# illum_pos = illuminator.get_position_from_time(sensor_state)
		# illum_vel = illuminator.get_velocity(sensor_state)

		# body_rot = body.rotation(sensor_state)

		# sun_av = utils.unit_vector(utils.cross_product(illum_pos, illum_vel))

		# npole = [body_rot[6], body_rot[7], body_rot[8]]

		# z = sun_av
		# x = utils.unit_vector(utils.cross_product(npole, z))
		# y = utils.unit_vector(utils.cross_product(z, x))

		# trans = np.matrix([[x.x, x.y, x.z], [y.x, y.y, y.z], [z.x, z.y, z.z]])

		# pos = np.matmul(trans, illum_pos)
		# spherical_pos = utils.rect_to_spherical(pos)

		# longitude360 = np.rad2deg(spherical_pos.lon)

		# if (longitude360 != 360.0):
		# longitude360 -= 360 * np.floor(longitude360 / 360)

		# return longitude360
		No newline at end of file

knoten/shape.py

0 → 100644

+58 −0

Original line number	Diff line number	Diff line
		from knoten import csm, utils
		import spiceypy as spice
		import numpy as np
		import csmapi

		class Ellipsoid:
		"""
		A biaxial ellipsoid shape model.
		"""

		def __init__(self, semi_major, semi_minor = None):
		"""
		Create an ellipsoid shape model from a set of radii

		Parameters
		----------
		semi_major : float
		The equatorial semi-major radius of the ellipsoid.
		semi_minor : float
		The polar semi-minor radius of the ellipsoid.
		"""
		self.a = semi_major
		self.b = semi_major
		self.c = semi_major

		if semi_minor is not None:
		self.c = semi_minor

		@classmethod
		def from_csm_sensor(cls, sensor):
		semi_major, semi_minor = csm.get_radii(sensor)
		return cls(semi_major, semi_minor)

		def get_surface_normal(self, ground_pt):
		"""
		Given a ground point, calculate the surface normal.

		Parameters
		----------
		ground_pt: tuple
		The ground point as an (x, y, z) tuple

		Returns
		-------
		: tuple
		in the form (x, y, z)
		"""
		normal = spice.surfnm(self.a, self.b, self.c, np.array([ground_pt.x, ground_pt.y, ground_pt.z]))
		return utils.Point(normal[0], normal[1], normal[2])


		def intersect_surface(self, sensor_pos, look_vec):
		sensor_pos = np.array([sensor_pos.x, sensor_pos.y, sensor_pos.z])
		look_vec = np.array([look_vec.x, look_vec.y, look_vec.z])

		ground_pt = spice.surfpt(sensor_pos, look_vec, self.a, self.b, self.c)
		return csmapi.EcefCoord(ground_pt[0], ground_pt[1], ground_pt[2])