From b233cdd365cb7ee455124e6afe5623979d3a0654 Mon Sep 17 00:00:00 2001 From: Amy Stamile Date: Wed, 24 Apr 2024 12:39:18 -0700 Subject: [PATCH 1/5] added sensor util functions --- examples/sensor_utils.ipynb | 316 ++++++++++++++++++++++++++++++++++++ knoten/bundle.py | 29 ++++ knoten/csm.py | 52 ++++++ knoten/sensor_utils.py | 296 +++++++++++++++++++++++++++++++++ tests/test_bundle.py | 9 + tests/test_csm.py | 31 +++- tests/test_sensorutils.py | 103 ++++++++++++ 7 files changed, 834 insertions(+), 2 deletions(-) create mode 100644 examples/sensor_utils.ipynb create mode 100644 knoten/sensor_utils.py create mode 100644 tests/test_sensorutils.py diff --git a/examples/sensor_utils.ipynb b/examples/sensor_utils.ipynb new file mode 100644 index 0000000..582920f --- /dev/null +++ b/examples/sensor_utils.ipynb @@ -0,0 +1,316 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Sensor Utils\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "\n", + "from csmapi import csmapi\n", + "from knoten import csm, sensor_utils\n", + "\n", + "import ale\n", + "import json" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Create a usgscsm sensor model" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/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.\n", + " warnings.warn(\n" + ] + } + ], + "source": [ + "fileName = \"data/N1573082850_1.cub\"\n", + "\n", + "kernels = ale.util.generate_kernels_from_cube(fileName, expand=True)\n", + "isd_string = ale.loads(fileName, props={'kernels': kernels})\n", + "csm_isd = os.path.splitext(fileName)[0] + '.json'\n", + "\n", + "with open(csm_isd, 'w') as isd_file:\n", + " isd_file.write(isd_string)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run Sensor Utils with usgscsm sensor model and image point" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "camera = csm.create_csm(csm_isd)\n", + "image_pt = csmapi.ImageCoord(511.5, 511.5)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "38.87212509629893" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "phaseAngle = sensor_utils.phase_angle(image_pt, camera)\n", + "\n", + "phaseAngle" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "49.60309924896046" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "emissionAngle = sensor_utils.emission_angle(image_pt, camera)\n", + "\n", + "emissionAngle" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2903512972.146144" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "slantDistance = sensor_utils.slant_distance(image_pt, camera)\n", + "\n", + "slantDistance" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2943536048.858226" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "targetCenterDistance = sensor_utils.target_center_distance(image_pt, camera)\n", + "\n", + "targetCenterDistance" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "LatLon(lat=3.2229625890973583, lon=258.6197326526089)" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "subSpacecraftPoint = sensor_utils.sub_spacecraft_point(image_pt, camera)\n", + "\n", + "subSpacecraftPoint" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "59096282.02424558" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "localRadius = sensor_utils.local_radius(image_pt, camera)\n", + "\n", + "localRadius" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(79.34815579474038, -2.7790780986459485)" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "rightAscDec = sensor_utils.right_ascension_declination(image_pt, camera)\n", + "\n", + "rightAscDec" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "17397.960941876583" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lineResolution = sensor_utils.line_resolution(image_pt, camera)\n", + "\n", + "lineResolution" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "17397.933700407997" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sampleResolution = sensor_utils.sample_resolution(image_pt, camera)\n", + "\n", + "sampleResolution" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "17397.94732114229" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "pixelResolution = sensor_utils.pixel_resolution(image_pt, camera)\n", + "\n", + "pixelResolution" + ] + } + ], + "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.18" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/knoten/bundle.py b/knoten/bundle.py index d59250a..9d1389c 100644 --- a/knoten/bundle.py +++ b/knoten/bundle.py @@ -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 diff --git a/knoten/csm.py b/knoten/csm.py index 9912283..0a569d3 100644 --- a/knoten/csm.py +++ b/knoten/csm.py @@ -9,6 +9,7 @@ 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 @@ -73,6 +96,35 @@ def create_camera(label, url='http://pfeffer.wr.usgs.gov/api/1.0/pds/'): if plugin.canModelBeConstructedFromISD(isd, model_name): 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: diff --git a/knoten/sensor_utils.py b/knoten/sensor_utils.py new file mode 100644 index 0000000..812fecf --- /dev/null +++ b/knoten/sensor_utils.py @@ -0,0 +1,296 @@ +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 diff --git a/tests/test_bundle.py b/tests/test_bundle.py index 6b95d67..eac4c17 100644 --- a/tests/test_bundle.py +++ b/tests/test_bundle.py @@ -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)] diff --git a/tests/test_csm.py b/tests/test_csm.py index d103207..ede5308 100644 --- a/tests/test_csm.py +++ b/tests/test_csm.py @@ -5,7 +5,7 @@ from plio.io.io_gdal import GeoDataset import csmapi -from knoten import csm, surface +from knoten import csm, surface, utils @pytest.fixture def mock_dem(): @@ -66,4 +66,31 @@ def test__compute_intersection_distance(): pt1 = Point(0,0,0) pt2 = Point(1,1,1) dist = csm._compute_intersection_distance(pt1, pt2) - assert dist == 1 \ No newline at end of file + assert dist == 1 + + +def test_get_state(mock_sensor, pt): + Locus = namedtuple("Locus", 'point direction') + + mock_sensor.getImageTime.return_value = 0.0 + mock_sensor.imageToRemoteImagingLocus.return_value = Locus(utils.Point(0, 1, 2), utils.Point(0, 1, 2)) + mock_sensor.getSensorPosition.return_value = utils.Point(2, 2, 2) + + state = csm.get_state(mock_sensor, pt) + + expected = { + "lookVec": utils.Point(0, 1, 2), + "sensorPos": utils.Point(2, 2, 2), + "sensorTime": 0.0, + "imagePoint": pt + } + + assert state == expected + + +@mock.patch.object(csm, 'get_radii', return_value=(10, 10)) +def test_get_surface_normal(mock_sensor): + ground_pt = utils.Point(1, 0, 0) + normal = csm.get_surface_normal(mock_sensor, ground_pt) + + assert normal == (1.0, 0.0, 0.0) \ No newline at end of file diff --git a/tests/test_sensorutils.py b/tests/test_sensorutils.py new file mode 100644 index 0000000..c7be63a --- /dev/null +++ b/tests/test_sensorutils.py @@ -0,0 +1,103 @@ +from unittest import mock +import pytest + +import numpy as np +import csmapi +from knoten import csm, sensor_utils, utils, bundle + +@pytest.fixture +def mock_sensor(): + mock_sensor = mock.MagicMock(spec=csmapi.RasterGM) + return mock_sensor + +@pytest.fixture +def pt(): + return csmapi.ImageCoord(0.0, 0.0) + +@mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(100.0, 0.0, 0.0)}) +@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(0,0,0)) +def test_phase_angle(mock_sensor, pt): + mock_sensor.getIlluminationDirection.return_value = utils.Point(100.0, 100.0, 0.0) + phase_angle = sensor_utils.phase_angle(pt, mock_sensor) + + np.testing.assert_array_equal(phase_angle, 135.0) + + +@mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(100.0, 0.0, 0.0)}) +@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(0,0,0)) +def test_emission_angle(mock_sensor, pt): + with mock.patch.object(csm, 'get_surface_normal', return_value=utils.Point(-1,0,0)) as mock_normal: + emission_angle = sensor_utils.emission_angle(pt, mock_sensor) + + mock_normal.assert_called() + + np.testing.assert_array_equal(emission_angle, 180.0) + + +@mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(-100.0, 0.0, 0.0)}) +@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(0,0,0)) +def test_slant_distance(mock_sensor, pt): + slant_distance = sensor_utils.slant_distance(pt, mock_sensor) + + np.testing.assert_array_equal(slant_distance, 100.0) + + +@mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(-100.0, 0.0, 0.0)}) +def test_target_center_distance(mock_sensor, pt): + target_center_distance = sensor_utils.target_center_distance(pt, mock_sensor) + + np.testing.assert_array_equal(target_center_distance, 100.0) + + +@mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(0.0, 0.0, 100.0)}) +def test_sub_spacecraft_point(mock_sensor, pt): + sub_spacecraft_point = sensor_utils.sub_spacecraft_point(pt, mock_sensor) + + np.testing.assert_array_equal(sub_spacecraft_point, [90.0, 0.0]) + + +@mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(100.0, 0.0, 0.0)}) +@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(10.0, 0.0, 0.0)) +def test_local_radius_intersection(mock_sensor, pt): + local_radius = sensor_utils.local_radius(pt, mock_sensor) + + np.testing.assert_array_equal(local_radius, 10.0) + + +@mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(-1000.0, 0.0, 0.0)}) +@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(10.0, 0.0, 0.0)) +def test_local_radius_ground(mock_sensor, pt): + local_radius = sensor_utils.local_radius(pt, mock_sensor) + + np.testing.assert_array_equal(local_radius, 10.0) + + +@mock.patch.object(csm, 'get_state', return_value={'lookVec': utils.Point(-1.0, 0.0, 0.0)}) +def test_right_ascension_declination(mock_sensor, pt): + right_ascension_declination = sensor_utils.right_ascension_declination(pt, mock_sensor) + + np.testing.assert_array_equal(right_ascension_declination, [180.0, 0.0]) + + +@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(0.0, 0.0, 0.0)) +@mock.patch.object(bundle, 'compute_image_partials', return_value=np.array([2, 1, 4, 4, 4, 8])) +def test_line_resolution(mock_sensor, pt): + line_resolution = sensor_utils.line_resolution(pt, mock_sensor) + + np.testing.assert_array_equal(line_resolution, 6.0) + + +@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(0.0, 0.0, 0.0)) +@mock.patch.object(bundle, 'compute_image_partials', return_value=np.array([2, 1, 4, 4, 4, 8])) +def test_sample_resolution(mock_sensor, pt): + sample_resolution = sensor_utils.sample_resolution(pt, mock_sensor) + + np.testing.assert_array_equal(sample_resolution, 9.0) + + +@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(0.0, 0.0, 0.0)) +@mock.patch.object(bundle, 'compute_image_partials', return_value=np.array([2, 1, 4, 4, 4, 8])) +def test_pixel_resolution(mock_sensor, pt): + pixel_resolution = sensor_utils.pixel_resolution(pt, mock_sensor) + + np.testing.assert_array_equal(pixel_resolution, 7.5) \ No newline at end of file From 2257afbe8cc5212ef7d181463268a151305891cd Mon Sep 17 00:00:00 2001 From: Amy Stamile Date: Thu, 25 Apr 2024 13:06:12 -0700 Subject: [PATCH 2/5] addressed PR feedback --- examples/sensor_utils.ipynb | 51 +++++++------- knoten/csm.py | 26 +------- knoten/illuminator.py | 23 +++++++ knoten/sensor_utils.py | 130 ++++++++++++++++++++++++++++++------ knoten/shape.py | 58 ++++++++++++++++ 5 files changed, 220 insertions(+), 68 deletions(-) create mode 100644 knoten/illuminator.py create mode 100644 knoten/shape.py diff --git a/examples/sensor_utils.ipynb b/examples/sensor_utils.ipynb index 582920f..8a6786c 100644 --- a/examples/sensor_utils.ipynb +++ b/examples/sensor_utils.ipynb @@ -15,9 +15,15 @@ "source": [ "import os\n", "\n", + "os.environ[\"ISISROOT\"] = \"/Users/astamile/ISIS3/build\"\n", + "os.environ[\"ISISDATA\"] = \"/Volumes/isis_data1/isis_data/\"\n", + "\n", "from csmapi import csmapi\n", "from knoten import csm, sensor_utils\n", "\n", + "from knoten.shape import Ellipsoid\n", + "from knoten.illuminator import Illuminator\n", + "\n", "import ale\n", "import json" ] @@ -31,18 +37,9 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/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.\n", - " warnings.warn(\n" - ] - } - ], + "outputs": [], "source": [ "fileName = \"data/N1573082850_1.cub\"\n", "\n", @@ -68,7 +65,9 @@ "outputs": [], "source": [ "camera = csm.create_csm(csm_isd)\n", - "image_pt = csmapi.ImageCoord(511.5, 511.5)" + "image_pt = csmapi.ImageCoord(511.5, 511.5)\n", + "shape = Ellipsoid.from_csm_sensor(camera)\n", + "illuminator = Illuminator()" ] }, { @@ -79,7 +78,7 @@ { "data": { "text/plain": [ - "38.87212509629893" + "38.87212509629895" ] }, "execution_count": 4, @@ -88,7 +87,7 @@ } ], "source": [ - "phaseAngle = sensor_utils.phase_angle(image_pt, camera)\n", + "phaseAngle = sensor_utils.phase_angle(image_pt, camera, shape, illuminator)\n", "\n", "phaseAngle" ] @@ -101,7 +100,7 @@ { "data": { "text/plain": [ - "49.60309924896046" + "49.60309924893989" ] }, "execution_count": 5, @@ -110,7 +109,7 @@ } ], "source": [ - "emissionAngle = sensor_utils.emission_angle(image_pt, camera)\n", + "emissionAngle = sensor_utils.emission_angle(image_pt, camera, shape)\n", "\n", "emissionAngle" ] @@ -123,7 +122,7 @@ { "data": { "text/plain": [ - "2903512972.146144" + "2903512972.146115" ] }, "execution_count": 6, @@ -132,7 +131,7 @@ } ], "source": [ - "slantDistance = sensor_utils.slant_distance(image_pt, camera)\n", + "slantDistance = sensor_utils.slant_distance(image_pt, camera, shape)\n", "\n", "slantDistance" ] @@ -189,7 +188,7 @@ { "data": { "text/plain": [ - "59096282.02424558" + "59096282.024265066" ] }, "execution_count": 9, @@ -198,7 +197,7 @@ } ], "source": [ - "localRadius = sensor_utils.local_radius(image_pt, camera)\n", + "localRadius = sensor_utils.local_radius(image_pt, camera, shape)\n", "\n", "localRadius" ] @@ -233,7 +232,7 @@ { "data": { "text/plain": [ - "17397.960941876583" + "17397.96094194587" ] }, "execution_count": 11, @@ -242,7 +241,7 @@ } ], "source": [ - "lineResolution = sensor_utils.line_resolution(image_pt, camera)\n", + "lineResolution = sensor_utils.line_resolution(image_pt, camera, shape)\n", "\n", "lineResolution" ] @@ -255,7 +254,7 @@ { "data": { "text/plain": [ - "17397.933700407997" + "17397.93370038153" ] }, "execution_count": 12, @@ -264,7 +263,7 @@ } ], "source": [ - "sampleResolution = sensor_utils.sample_resolution(image_pt, camera)\n", + "sampleResolution = sensor_utils.sample_resolution(image_pt, camera, shape)\n", "\n", "sampleResolution" ] @@ -277,7 +276,7 @@ { "data": { "text/plain": [ - "17397.94732114229" + "17397.9473211637" ] }, "execution_count": 13, @@ -286,7 +285,7 @@ } ], "source": [ - "pixelResolution = sensor_utils.pixel_resolution(image_pt, camera)\n", + "pixelResolution = sensor_utils.pixel_resolution(image_pt, camera, shape)\n", "\n", "pixelResolution" ] diff --git a/knoten/csm.py b/knoten/csm.py index 0a569d3..3cb7f53 100644 --- a/knoten/csm.py +++ b/knoten/csm.py @@ -9,7 +9,6 @@ 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) diff --git a/knoten/illuminator.py b/knoten/illuminator.py new file mode 100644 index 0000000..09292ae --- /dev/null +++ b/knoten/illuminator.py @@ -0,0 +1,23 @@ +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 diff --git a/knoten/sensor_utils.py b/knoten/sensor_utils.py index 812fecf..cfadaee 100644 --- a/knoten/sensor_utils.py +++ b/knoten/sensor_utils.py @@ -2,7 +2,7 @@ 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 \ No newline at end of file + 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 diff --git a/knoten/shape.py b/knoten/shape.py new file mode 100644 index 0000000..509e15d --- /dev/null +++ b/knoten/shape.py @@ -0,0 +1,58 @@ +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]) + From 92116f14860980326f4fe4361f2f00163a6afcd2 Mon Sep 17 00:00:00 2001 From: Amy Stamile Date: Thu, 25 Apr 2024 19:14:33 -0700 Subject: [PATCH 3/5] fixed tests. --- knoten/sensor_utils.py | 4 +- knoten/shape.py | 6 ++- knoten/utils.py | 1 + tests/test_csm.py | 10 +--- tests/test_illuminator.py | 24 +++++++++ tests/test_sensorutils.py | 108 +++++++++++++++++++++++--------------- tests/test_shape.py | 28 ++++++++++ 7 files changed, 127 insertions(+), 54 deletions(-) create mode 100644 tests/test_illuminator.py create mode 100644 tests/test_shape.py diff --git a/knoten/sensor_utils.py b/knoten/sensor_utils.py index cfadaee..be985aa 100644 --- a/knoten/sensor_utils.py +++ b/knoten/sensor_utils.py @@ -34,7 +34,7 @@ def phase_angle(image_pt, sensor, shape, illuminator): sensor_state = csm.get_state(sensor, image_pt) ground_pt = shape.intersect_surface(sensor_state["sensorPos"], sensor_state["lookVec"]) - + illum_pos = illuminator.get_position_from_csm_sensor(sensor, ground_pt) vec_a = utils.Point(sensor_state["sensorPos"].x - ground_pt.x, @@ -44,7 +44,7 @@ def phase_angle(image_pt, sensor, shape, illuminator): 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, shape): diff --git a/knoten/shape.py b/knoten/shape.py index 509e15d..441ab06 100644 --- a/knoten/shape.py +++ b/knoten/shape.py @@ -8,7 +8,7 @@ class Ellipsoid: A biaxial ellipsoid shape model. """ - def __init__(self, semi_major, semi_minor = None): + def __init__(self, semi_major, semi_minor=None, median=None): """ Create an ellipsoid shape model from a set of radii @@ -23,9 +23,13 @@ def __init__(self, semi_major, semi_minor = None): self.b = semi_major self.c = semi_major + if median is not None: + self.b = median + 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) diff --git a/knoten/utils.py b/knoten/utils.py index cde4fe5..1165b85 100644 --- a/knoten/utils.py +++ b/knoten/utils.py @@ -29,6 +29,7 @@ def sep_angle(a_vec, b_vec): : np.ndarray """ dot_prod = a_vec.x * b_vec.x + a_vec.y * b_vec.y + a_vec.z * b_vec.z + print(dot_prod) dot_prod /= magnitude(a_vec) * magnitude(b_vec) if(dot_prod >= 1.0): return 0.0 diff --git a/tests/test_csm.py b/tests/test_csm.py index ede5308..b31e9b9 100644 --- a/tests/test_csm.py +++ b/tests/test_csm.py @@ -85,12 +85,4 @@ def test_get_state(mock_sensor, pt): "imagePoint": pt } - assert state == expected - - -@mock.patch.object(csm, 'get_radii', return_value=(10, 10)) -def test_get_surface_normal(mock_sensor): - ground_pt = utils.Point(1, 0, 0) - normal = csm.get_surface_normal(mock_sensor, ground_pt) - - assert normal == (1.0, 0.0, 0.0) \ No newline at end of file + assert state == expected \ No newline at end of file diff --git a/tests/test_illuminator.py b/tests/test_illuminator.py new file mode 100644 index 0000000..e7d8688 --- /dev/null +++ b/tests/test_illuminator.py @@ -0,0 +1,24 @@ +import unittest +from unittest import mock + +import csmapi + +from knoten import utils +from knoten.illuminator import Illuminator + +class TestIlluminator(unittest.TestCase): + + def test_get_position_from_csm_sensor(self): + mock_sensor = mock.MagicMock(spec=csmapi.RasterGM) + mock_sensor.getIlluminationDirection.return_value = utils.Point(1.0, 10.0, 0.0) + + test_illum = Illuminator() + + ground_pt = utils.Point(10.0, 10.0, 0.0) + + position = test_illum.get_position_from_csm_sensor(mock_sensor, ground_pt) + + self.assertEqual(position, (9.0, 0.0, 0.0)) + + def tearDown(self): + pass \ No newline at end of file diff --git a/tests/test_sensorutils.py b/tests/test_sensorutils.py index c7be63a..b6731c2 100644 --- a/tests/test_sensorutils.py +++ b/tests/test_sensorutils.py @@ -3,7 +3,7 @@ import numpy as np import csmapi -from knoten import csm, sensor_utils, utils, bundle +from knoten import csm, sensor_utils, utils, bundle, shape, illuminator @pytest.fixture def mock_sensor(): @@ -14,32 +14,46 @@ def mock_sensor(): def pt(): return csmapi.ImageCoord(0.0, 0.0) -@mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(100.0, 0.0, 0.0)}) -@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(0,0,0)) -def test_phase_angle(mock_sensor, pt): - mock_sensor.getIlluminationDirection.return_value = utils.Point(100.0, 100.0, 0.0) - phase_angle = sensor_utils.phase_angle(pt, mock_sensor) +@pytest.fixture +def test_shape(): + return shape.Ellipsoid(0, 0) - np.testing.assert_array_equal(phase_angle, 135.0) +@pytest.fixture +def test_illuminator(): + return illuminator.Illuminator(0, 0) -@mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(100.0, 0.0, 0.0)}) -@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(0,0,0)) -def test_emission_angle(mock_sensor, pt): - with mock.patch.object(csm, 'get_surface_normal', return_value=utils.Point(-1,0,0)) as mock_normal: - emission_angle = sensor_utils.emission_angle(pt, mock_sensor) +def test_phase_angle(mock_sensor, pt, test_shape, test_illuminator): + with ( + mock.patch.object(illuminator.Illuminator, 'get_position_from_csm_sensor', return_value=utils.Point(100.0, 100.0, 0.0)), + mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(100.0, 0.0, 0.0), 'lookVec': utils.Point(-1.0, 0.0, 0.0)}), + mock.patch.object(shape.Ellipsoid, 'intersect_surface', return_value=utils.Point(0,0,0)) + ): + phase_angle = sensor_utils.phase_angle(pt, mock_sensor, test_shape, test_illuminator) - mock_normal.assert_called() + np.testing.assert_allclose(phase_angle, 45.0) - np.testing.assert_array_equal(emission_angle, 180.0) +def test_emission_angle(mock_sensor, pt, test_shape): + with ( + mock.patch.object(shape.Ellipsoid, 'get_surface_normal', return_value=utils.Point(-1,0,0)), + mock.patch.object(shape.Ellipsoid, 'intersect_surface', return_value=utils.Point(0,0,0)), + mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(100.0, 0.0, 0.0), 'lookVec': utils.Point(-1.0, 0.0, 0.0)}), + mock.patch.object(illuminator.Illuminator, 'get_position_from_csm_sensor', return_value=utils.Point(100.0, 100.0, 0.0)) + ): + emission_angle = sensor_utils.emission_angle(pt, mock_sensor, test_shape) -@mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(-100.0, 0.0, 0.0)}) -@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(0,0,0)) -def test_slant_distance(mock_sensor, pt): - slant_distance = sensor_utils.slant_distance(pt, mock_sensor) + np.testing.assert_array_equal(emission_angle, 180.0) - np.testing.assert_array_equal(slant_distance, 100.0) + +def test_slant_distance(mock_sensor, pt, test_shape): + with ( + mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(100.0, 0.0, 0.0), 'lookVec': utils.Point(-1.0, 0.0, 0.0)}), + mock.patch.object(shape.Ellipsoid, 'intersect_surface', return_value=utils.Point(0,0,0)) + ): + slant_distance = sensor_utils.slant_distance(pt, mock_sensor, test_shape) + + np.testing.assert_array_equal(slant_distance, 100.0) @mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(-100.0, 0.0, 0.0)}) @@ -56,18 +70,18 @@ def test_sub_spacecraft_point(mock_sensor, pt): np.testing.assert_array_equal(sub_spacecraft_point, [90.0, 0.0]) -@mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(100.0, 0.0, 0.0)}) -@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(10.0, 0.0, 0.0)) -def test_local_radius_intersection(mock_sensor, pt): - local_radius = sensor_utils.local_radius(pt, mock_sensor) +@mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(100.0, 0.0, 0.0), 'lookVec': utils.Point(-1.0, 0.0, 0.0)}) +@mock.patch.object(shape.Ellipsoid, 'intersect_surface', return_value=utils.Point(10.0, 0.0, 0.0)) +def test_local_radius_intersection(mock_sensor, pt, test_shape): + local_radius = sensor_utils.local_radius(pt, mock_sensor, test_shape) np.testing.assert_array_equal(local_radius, 10.0) -@mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(-1000.0, 0.0, 0.0)}) -@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(10.0, 0.0, 0.0)) -def test_local_radius_ground(mock_sensor, pt): - local_radius = sensor_utils.local_radius(pt, mock_sensor) +@mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(1000.0, 0.0, 0.0), 'lookVec': utils.Point(-1000.0, 0.0, 0.0)}) +@mock.patch.object(shape.Ellipsoid, 'intersect_surface', return_value=utils.Point(10.0, 0.0, 0.0)) +def test_local_radius_ground(mock_sensor, pt, test_shape): + local_radius = sensor_utils.local_radius(pt, mock_sensor, test_shape) np.testing.assert_array_equal(local_radius, 10.0) @@ -79,25 +93,35 @@ def test_right_ascension_declination(mock_sensor, pt): np.testing.assert_array_equal(right_ascension_declination, [180.0, 0.0]) -@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(0.0, 0.0, 0.0)) -@mock.patch.object(bundle, 'compute_image_partials', return_value=np.array([2, 1, 4, 4, 4, 8])) -def test_line_resolution(mock_sensor, pt): - line_resolution = sensor_utils.line_resolution(pt, mock_sensor) +def test_line_resolution(mock_sensor, pt, test_shape): + with ( + mock.patch.object(bundle, 'compute_image_partials', return_value=np.array([2, 1, 4, 4, 4, 8])), + mock.patch.object(shape.Ellipsoid, 'intersect_surface', return_value=utils.Point(0,0,0)), + mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(100.0, 0.0, 0.0), 'lookVec': utils.Point(-1.0, 0.0, 0.0)}) + ): + line_resolution = sensor_utils.line_resolution(pt, mock_sensor, test_shape) + + np.testing.assert_array_equal(line_resolution, 6.0) - np.testing.assert_array_equal(line_resolution, 6.0) -@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(0.0, 0.0, 0.0)) -@mock.patch.object(bundle, 'compute_image_partials', return_value=np.array([2, 1, 4, 4, 4, 8])) -def test_sample_resolution(mock_sensor, pt): - sample_resolution = sensor_utils.sample_resolution(pt, mock_sensor) +def test_sample_resolution(mock_sensor, pt, test_shape): + with ( + mock.patch.object(bundle, 'compute_image_partials', return_value=np.array([2, 1, 4, 4, 4, 8])), + mock.patch.object(shape.Ellipsoid, 'intersect_surface', return_value=utils.Point(0,0,0)), + mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(100.0, 0.0, 0.0), 'lookVec': utils.Point(-1.0, 0.0, 0.0)}) + ): + sample_resolution = sensor_utils.sample_resolution(pt, mock_sensor, test_shape) - np.testing.assert_array_equal(sample_resolution, 9.0) + np.testing.assert_array_equal(sample_resolution, 9.0) -@mock.patch.object(csm, 'generate_ground_point', return_value=utils.Point(0.0, 0.0, 0.0)) -@mock.patch.object(bundle, 'compute_image_partials', return_value=np.array([2, 1, 4, 4, 4, 8])) -def test_pixel_resolution(mock_sensor, pt): - pixel_resolution = sensor_utils.pixel_resolution(pt, mock_sensor) +def test_pixel_resolution(mock_sensor, pt, test_shape): + with ( + mock.patch.object(bundle, 'compute_image_partials', return_value=np.array([2, 1, 4, 4, 4, 8])), + mock.patch.object(shape.Ellipsoid, 'intersect_surface', return_value=utils.Point(0,0,0)), + mock.patch.object(csm, 'get_state', return_value={'sensorPos': utils.Point(100.0, 0.0, 0.0), 'lookVec': utils.Point(-1.0, 0.0, 0.0)}) + ): + pixel_resolution = sensor_utils.pixel_resolution(pt, mock_sensor, test_shape) - np.testing.assert_array_equal(pixel_resolution, 7.5) \ No newline at end of file + np.testing.assert_array_equal(pixel_resolution, 7.5) \ No newline at end of file diff --git a/tests/test_shape.py b/tests/test_shape.py new file mode 100644 index 0000000..07cb559 --- /dev/null +++ b/tests/test_shape.py @@ -0,0 +1,28 @@ +import unittest +from unittest import mock + +from knoten import shape, utils + +class TestEllipsoid(unittest.TestCase): + + def test_get_surface_normal(self): + test_shape = shape.Ellipsoid(10, 10) + + ground_pt = utils.Point(1, 0, 0) + normal = test_shape.get_surface_normal(ground_pt) + + self.assertEqual(normal, (1.0, 0.0, 0.0)) + + def test_intersect_surface(self): + test_shape = shape.Ellipsoid(10, 10) + sensor_pos = utils.Point(100, 0, 0) + look_vec = utils.Point(-1, 0, 0) + + intersection = test_shape.intersect_surface(sensor_pos, look_vec) + + self.assertEqual(intersection.x, 10.0) + self.assertEqual(intersection.y, 0.0) + self.assertEqual(intersection.z, 0.0) + + def tearDown(self): + pass \ No newline at end of file From e35fc471545ab48116ded6125ddf9daeccf7e9cd Mon Sep 17 00:00:00 2001 From: Amy Stamile Date: Fri, 26 Apr 2024 08:10:30 -0700 Subject: [PATCH 4/5] removed full local paths --- examples/sensor_utils.ipynb | 3 --- 1 file changed, 3 deletions(-) diff --git a/examples/sensor_utils.ipynb b/examples/sensor_utils.ipynb index 8a6786c..f2d3bc0 100644 --- a/examples/sensor_utils.ipynb +++ b/examples/sensor_utils.ipynb @@ -15,9 +15,6 @@ "source": [ "import os\n", "\n", - "os.environ[\"ISISROOT\"] = \"/Users/astamile/ISIS3/build\"\n", - "os.environ[\"ISISDATA\"] = \"/Volumes/isis_data1/isis_data/\"\n", - "\n", "from csmapi import csmapi\n", "from knoten import csm, sensor_utils\n", "\n", From 6e06809e431646912426ed0a653bd2c00819fa78 Mon Sep 17 00:00:00 2001 From: Amy Stamile Date: Fri, 26 Apr 2024 08:36:35 -0700 Subject: [PATCH 5/5] clean up docs --- knoten/shape.py | 11 +++++++---- knoten/utils.py | 1 - 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/knoten/shape.py b/knoten/shape.py index 441ab06..78a406a 100644 --- a/knoten/shape.py +++ b/knoten/shape.py @@ -5,7 +5,7 @@ class Ellipsoid: """ - A biaxial ellipsoid shape model. + A biaxial or triaxial ellipsoid shape model. """ def __init__(self, semi_major, semi_minor=None, median=None): @@ -15,9 +15,11 @@ def __init__(self, semi_major, semi_minor=None, median=None): Parameters ---------- semi_major : float - The equatorial semi-major radius of the ellipsoid. + Length of ellipsoid semi-axis along the x-axis. semi_minor : float - The polar semi-minor radius of the ellipsoid. + Length of ellipsoid semi-axis along the z-axis. + median : float + Length of ellipsoid semi-axis along the y-axis. """ self.a = semi_major self.b = semi_major @@ -34,7 +36,8 @@ def __init__(self, semi_major, semi_minor=None, median=None): 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. diff --git a/knoten/utils.py b/knoten/utils.py index 1165b85..cde4fe5 100644 --- a/knoten/utils.py +++ b/knoten/utils.py @@ -29,7 +29,6 @@ def sep_angle(a_vec, b_vec): : np.ndarray """ dot_prod = a_vec.x * b_vec.x + a_vec.y * b_vec.y + a_vec.z * b_vec.z - print(dot_prod) dot_prod /= magnitude(a_vec) * magnitude(b_vec) if(dot_prod >= 1.0): return 0.0