added sensor util functions (b233cdd3) · Commits · aflab / astrogeology / Knoten

examples/sensor_utils.ipynb

0 → 100644

+168 −0

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

		# Sensor Utils

		%% Cell type:code id: tags:

		``` python
		import os

		from csmapi import csmapi
		from knoten import csm, sensor_utils

		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)
		```

		%% Cell type:code id: tags:

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

		phaseAngle
		```

		%% Output

		38.87212509629893

		%% Cell type:code id: tags:

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

		emissionAngle
		```

		%% Output

		49.60309924896046

		%% Cell type:code id: tags:

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

		slantDistance
		```

		%% Output

		2903512972.146144

		%% 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
		```

		%% Output

		59096282.02424558

		%% 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
		```

		%% Output

		17397.960941876583

		%% Cell type:code id: tags:

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

		sampleResolution
		```

		%% Output

		17397.933700407997

		%% Cell type:code id: tags:

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

		pixelResolution
		```

		%% Output

		17397.94732114229

knoten/bundle.py

+29 −0

Original line number	Diff line number	Diff line
		@@ -212,6 +212,35 @@ def compute_ground_partials(sensor, ground_pt):
		partials = np.array(sensor.computeGroundPartials(csm_ground))
		return np.reshape(partials, (2, 3))

		def compute_image_partials(sensor, ground_pt):
		"""
		Compute the partial derivatives of the ground point with respect to
		the line and sample at a ground point.

		These are not normally available from the CSM model, so we use
		csm::RasterGM::computeGroundPartials to get the Jacobian of the ground to
		image transformation. Then we use the pseudoinverse of that to get the
		Jacobian of the image to ground transformation.

		Parameters
		----------
		sensor : CSM sensor
		The CSM sensor model
		ground_pt : array
		The (x, y, z) ground point to compute the partial derivatives W.R.T.

		Returns
		-------
		: array
		The partial derivatives of the image to ground transformation
		"""
		if isinstance(ground_pt, csmapi.EcefCoord):
		ground_pt = [ground_pt.x, ground_pt.y, ground_pt.z]
		ground_matrix = compute_ground_partials(sensor, ground_pt)
		image_matrix = np.linalg.pinv(ground_matrix)

		return image_matrix.flatten()

		def compute_coefficient_columns(network, sensors, parameters):
		"""
		Compute the columns for different coefficients

knoten/csm.py

+52 −0

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

		from knoten.surface import EllipsoidDem

		@@ -41,6 +42,28 @@ 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
		@@ -74,6 +97,35 @@ def create_camera(label, url='http://pfeffer.wr.usgs.gov/api/1.0/pds/'):
		model = plugin.constructModelFromISD(isd, model_name)
		return model

		def get_state(sensor, image_pt):
		"""
		Get the state of the sensor model at a given image point.

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

		image_pt : tuple
		Pair of x, y (sample, line) coordinates in pixel space

		Returns
		-------
		: dict
		Dictionary containing lookVec, sensorPos, sensorTime, and imagePoint
		"""
		sensor_time = sensor.getImageTime(image_pt)
		locus = sensor.imageToRemoteImagingLocus(image_pt)
		sensor_position = sensor.getSensorPosition(image_pt)

		sensor_state = {
		"lookVec": locus.direction,
		"sensorPos": sensor_position,
		"sensorTime": sensor_time,
		"imagePoint": image_pt
		}
		return sensor_state

		def _from_state(state, verbose):
		with open(state, 'r') as stream:
		model_name = stream.readline().rstrip()

knoten/sensor_utils.py

0 → 100644

+296 −0

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

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

		Phase Angle: The angle between the vector from the intersection point to
		the observer (usually the spacecraft) and the vector from the intersection
		point to the illuminator (usually the sun).

		Parameters
		----------
		image_pt : tuple
		Pair of x, y (sample, line) coordinates in pixel space

		sensor : object
		A CSM compliant sensor model object

		Returns
		-------
		: np.ndarray
		phase angle in degrees
		"""
		if not isinstance(image_pt, csmapi.ImageCoord):
		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)

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

		vec_a = utils.Point(sensor_state["sensorPos"].x - ground_pt.x,
		sensor_state["sensorPos"].y - ground_pt.y,
		sensor_state["sensorPos"].z - ground_pt.z)

		vec_b = utils.Point(illum_pos.x - ground_pt.x,
		illum_pos.y - ground_pt.y,
		illum_pos.z - ground_pt.z)

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

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

		Emission Angle: The angle between the surface normal vector at the
		intersection point and the vector from the intersection point to the
		observer (usually the spacecraft). The emission angle varies from 0 degrees
		when the observer is viewing the sub-spacecraft point (nadir viewing) to 90
		degrees when the intercept is tangent to the surface of the target body.
		Thus, higher values of emission angle indicate more oblique viewing of the
		target.

		Parameters
		----------
		image_pt : tuple
		Pair of x, y (sample, line) coordinates in pixel space

		sensor : object
		A CSM compliant sensor model object

		Returns
		-------
		: np.ndarray
		emission angle in degrees
		"""
		if not isinstance(image_pt, csmapi.ImageCoord):
		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)

		normal = csm.get_surface_normal(sensor, ground_pt)

		sensor_diff = utils.Point(sensor_state["sensorPos"].x - ground_pt.x,
		sensor_state["sensorPos"].y - ground_pt.y,
		sensor_state["sensorPos"].z - ground_pt.z)

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

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

		Parameters
		----------
		image_pt : tuple
		Pair of x, y (sample, line) coordinates in pixel space

		sensor : object
		A CSM compliant sensor model object

		Returns
		-------
		: np.ndarray
		slant distance in meters
		"""
		if not isinstance(image_pt, csmapi.ImageCoord):
		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)

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

		def target_center_distance(image_pt, sensor):
		"""
		Calculates and returns the distance from the spacecraft to the target center.

		Parameters
		----------
		image_pt : tuple
		Pair of x, y (sample, line) coordinates in pixel space

		sensor : object
		A CSM compliant sensor model object

		Returns
		-------
		: np.ndarray
		target center distance in meters
		"""
		if not isinstance(image_pt, csmapi.ImageCoord):
		image_pt = csmapi.ImageCoord(*image_pt)

		sensor_state = csm.get_state(sensor, image_pt)
		return utils.distance(sensor_state["sensorPos"], utils.Point(0,0,0))

		def sub_spacecraft_point(image_pt, sensor):
		"""
		Get the latitude and longitude of the sub-spacecraft point.

		Parameters
		----------
		image_pt : tuple
		Pair of x, y (sample, line) coordinates in pixel space

		sensor : object
		A CSM compliant sensor model object

		Returns
		-------
		: np.ndarray
		sub spacecraft point in degrees
		"""
		if not isinstance(image_pt, csmapi.ImageCoord):
		image_pt = csmapi.ImageCoord(*image_pt)

		sensor_state = csm.get_state(sensor, image_pt)
		lat_lon_rad = utils.rect_to_spherical(sensor_state["sensorPos"])

		return utils.radians_to_degrees(lat_lon_rad)

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

		Parameters
		----------
		image_pt : tuple
		Pair of x, y (sample, line) coordinates in pixel space

		sensor : object
		A CSM compliant sensor model object

		Returns
		-------
		: np.ndarray
		local radius in meters
		"""
		if not isinstance(image_pt, csmapi.ImageCoord):
		image_pt = csmapi.ImageCoord(*image_pt)

		ground_pt = csm.generate_ground_point(0.0, image_pt, sensor)
		return utils.magnitude(ground_pt)

		def right_ascension_declination(image_pt, sensor):
		"""
		Computes the right ascension and declination for a given image point.

		Parameters
		----------
		image_pt : tuple
		Pair of x, y (sample, line) coordinates in pixel space

		sensor : object
		A CSM compliant sensor model object

		Returns
		-------
		: tuple
		in the form (ra, dec) in degrees
		"""
		if not isinstance(image_pt, csmapi.ImageCoord):
		image_pt = csmapi.ImageCoord(*image_pt)

		sensor_state = csm.get_state(sensor, image_pt)
		spherical_pt = utils.rect_to_spherical(sensor_state["lookVec"])

		ra_dec = utils.radians_to_degrees(spherical_pt)

		return (ra_dec.lon, ra_dec.lat)

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

		CSM sensor models do not expose all of the necessary parameters to do the
		same calculation as ISIS sensor models, so this uses a more time consuming but
		more accurate method and thus is equivalent to the oblique line resolution.

		For time dependent sensor models, this may also be the line-to-line resolution
		and not the resolution within a line or framelet. This is determined by the
		CSM model's ground computeGroundPartials method.

		Parameters
		----------
		image_pt : tuple
		Pair of x, y (sample, line) coordinates in pixel space

		sensor : object
		A CSM compliant sensor model object

		Returns
		-------
		: np.ndarray
		line resolution in meters/pixel
		"""
		if not isinstance(image_pt, csmapi.ImageCoord):
		image_pt = csmapi.ImageCoord(*image_pt)

		ground_pt = csm.generate_ground_point(0.0, image_pt, sensor)
		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):
		"""
		Compute the sample resolution in meters per pixel for the current set point.

		CSM sensor models do not expose all of the necessary parameters to do the
		same calculation as ISIS sensor models, so this uses a more time consuming but
		more accurate method and thus is equivalent to the oblique sample resolution.

		Parameters
		----------
		image_pt : tuple
		Pair of x, y (sample, line) coordinates in pixel space

		sensor : object
		A CSM compliant sensor model object

		Returns
		-------
		: np.ndarray
		sample resolution in meters/pixel
		"""
		if not isinstance(image_pt, csmapi.ImageCoord):
		image_pt = csmapi.ImageCoord(*image_pt)

		ground_pt = csm.generate_ground_point(0.0, image_pt, sensor)
		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):
		"""
		Returns the pixel resolution at the current position in meters/pixel.

		Parameters
		----------
		image_pt : tuple
		Pair of x, y (sample, line) coordinates in pixel space

		sensor : object
		A CSM compliant sensor model object

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

tests/test_bundle.py

+9 −0

Original line number	Diff line number	Diff line
		@@ -121,6 +121,15 @@ def test_compute_ground_partials():
		partials = bundle.compute_ground_partials(sensor, ground_pt)
		np.testing.assert_array_equal(partials, [[1, 2, 3], [4, 5, 6]])

		def test_compute_image_partials():
		ground_pt = [9, 8, 10]
		sensor = mock.MagicMock(spec=csmapi.RasterGM)
		sensor.computeGroundPartials.return_value = (1, 2, 3, 4, 5, 6)
		partials = bundle.compute_image_partials(sensor, ground_pt)
		np.testing.assert_allclose(partials, np.array([-0.94444444, 0.44444444,
		-0.11111111, 0.11111111,
		0.72222222, -0.22222222]))

		def test_compute_jacobian(control_network, sensors):
		parameters = {sn: [mock.MagicMock()]*2 for sn in sensors}
		sensor_partials = [(i+1) * np.ones((2, 2)) for i in range(9)]