diff --git a/CHANGELOG.md b/CHANGELOG.md
index 07c497e74357af082d0142c43a2d37e2c039c8ea..8f7f78b07e86793ad4e952f14adf321be0bff3c1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -36,5 +36,7 @@ release.
## Unreleased
### Changed
-- Removed all `pyproj` calls from csm.py, abstracting them into the reprojection and pyproj.Transformer code inside utils.py. Updated the transformations to use the new pipeline style syntax to avoid deprecation warnings about old syntax.p
+- Removed all `pyproj` calls from csm.py, abstracting them into the reprojection and pyproj.Transformer code inside utils.py. Updated the transformations to use the new pipeline style syntax to avoid deprecation warnings about old syntax.
+### Fixed
+- Added a check to `generate_ground_point` when a GeoDataset is used to raise a `ValueError` if the algorithm intersects a no data value in the passed DEM. This ensures that valid heights are used in the intersection computation.
diff --git a/knoten/csm.py b/knoten/csm.py
index ee703142c56d4cb9eb3260429c0c932f7a8b83e6..f80753052a4b4d9c8ff55b02efd29a3e3102b348 100644
--- a/knoten/csm.py
+++ b/knoten/csm.py
@@ -157,24 +157,23 @@ def _(dem, image_pt, camera, max_its = 20, tolerance = 0.001):
intersection = generate_ground_point(0.0, image_pt, camera)
iterations = 0
semi_major, semi_minor = get_radii(camera)
-
+
source_proj = f'+proj=cart +a={semi_major} +b={semi_minor}'
dest_proj = f'+proj=lonlat +a={semi_major} +b={semi_minor}'
transformer = utils.create_transformer(source_proj, dest_proj)
-
while iterations != max_its:
- lon, lat, alt = transformer.transform(intersection.x,
+ lon, lat, _ = transformer.transform(intersection.x,
intersection.y,
intersection.z,
errcheck=True)
px, py = dem.latlon_to_pixel(lat, lon)
height = dem.read_array(1, [px, py, 1, 1])[0][0]
-
+ if height == dem.no_data_value:
+ raise ValueError(f'No DEM height at {lat}, {lon}')
+
next_intersection = camera.imageToGround(image_pt, float(height))
- dist = max(abs(intersection.x - next_intersection.x),
- abs(intersection.y - next_intersection.y),
- abs(intersection.z - next_intersection.z))
+ dist = _compute_intersection_distance(intersection, next_intersection)
intersection = next_intersection
iterations += 1
@@ -182,6 +181,30 @@ def _(dem, image_pt, camera, max_its = 20, tolerance = 0.001):
break
return intersection
+def _compute_intersection_distance(intersection, next_intersection):
+ """
+ Private func that takes two csmapi Ecef objects or other objects with
+ x,y,z properties and computes the distance between them. This is the
+ maximum distance in 3D space.
+
+ Parameters
+ ----------
+ intersection : object
+ Any object with x,y, and z properties that are numeric
+
+ next_intersection : object
+ Any object with x,y, and z properties that are numeric
+
+ Returns
+ -------
+ dist : float
+ The maximum distance between intersection and next_intersection
+ in one of the three planes (x,y,z)
+ """
+ return max(abs(intersection.x - next_intersection.x),
+ abs(intersection.y - next_intersection.y),
+ abs(intersection.z - next_intersection.z))
+
def generate_boundary(isize, npoints=10):
'''
Generates a bounding box given a camera model
diff --git a/tests/test_csm.py b/tests/test_csm.py
new file mode 100644
index 0000000000000000000000000000000000000000..94fb524ec5daeb83b5b779625c31e40932892f3a
--- /dev/null
+++ b/tests/test_csm.py
@@ -0,0 +1,67 @@
+from unittest import mock
+import pytest
+
+from plio.io.io_gdal import GeoDataset
+
+import csmapi
+from knoten import csm
+
+@pytest.fixture
+def mock_dem():
+ mock_dem = mock.MagicMock(spec_set=GeoDataset)
+ mock_dem.no_data_value = 10
+ mock_dem.read_array.return_value = [[100]]
+ mock_dem.latlon_to_pixel.return_value = (0.5,0.5)
+ return mock_dem
+
+@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)
+
+def test_generate_ground_point_with_float(mock_sensor):
+ csm.generate_ground_point(0, (0.5, 0.5), mock_sensor)
+ # The internal conversion from tuple to csmapi.ImageCoord means
+ # assert_called_once_with fails due to different addresses of
+ # different objects.
+ mock_sensor.imageToGround.assert_called_once()
+
+def test_generate_ground_point_with_imagecoord(mock_sensor, pt):
+ height = 0.0
+ csm.generate_ground_point(height, pt, mock_sensor)
+ mock_sensor.imageToGround.assert_called_once_with(pt, height)
+
+@mock.patch.object(csm, 'get_radii', return_value=(10,10))
+@mock.patch.object(csm, '_compute_intersection_distance', return_value=0)
+@mock.patch('pyproj.transformer.Transformer.transform', return_value=(0,0,0))
+def test_generate_ground_point_with_dtm(mock_get_radii,
+ mock_compute_intsesection,
+ mock_pyproj_transformer,
+ mock_sensor, pt, mock_dem):
+ csm.generate_ground_point(mock_dem, pt, mock_sensor)
+ # This call is mocked so that the intitial intersection and
+ # one iteration should occur. Therefore, the call count
+ # should always be 2.
+ assert mock_sensor.imageToGround.call_count == 2
+
+from collections import namedtuple
+
+@mock.patch.object(csm, 'get_radii', return_value=(10,10))
+@mock.patch('pyproj.transformer.Transformer.transform', return_value=(0,0,0))
+def test_generate_ground_point_with_dtm_ndv(mock_get_radii,
+ mock_pyproj_transformer,
+ mock_sensor, pt, mock_dem):
+ mock_dem.no_data_value = 100
+ with pytest.raises(ValueError):
+ csm.generate_ground_point(mock_dem, pt, mock_sensor)
+
+def test__compute_intersection_distance():
+ Point = namedtuple("Point", 'x, y, z')
+ 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