diff --git a/examples/sensor_utils.ipynb b/examples/sensor_utils.ipynb
index 582920f5e5d1fe9b9198397569fdc75cc32af720..8a6786c1b116ac682c8f103afe8344e05dcf52ce 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 0a569d30556e7eecc699121ca82751b1502600db..3cb7f535914ae567b6fbfdbff9ef4508536a4eb0 100644
--- a/knoten/csm.py
+++ b/knoten/csm.py
@@ -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)
diff --git a/knoten/illuminator.py b/knoten/illuminator.py
new file mode 100644
index 0000000000000000000000000000000000000000..09292ae1234092c8160536308a4926e3671e5fd0
--- /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 812fecf32804cbab0b408e8bb68821f4611119de..cfadaee2aa1b6d37d50e5d6e9d44161310103681 100644
--- a/knoten/sensor_utils.py
+++ b/knoten/sensor_utils.py
@@ -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
\ 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 0000000000000000000000000000000000000000..509e15da1768058955f1c813758bbed561d08b8b
--- /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])
+