Merge pull request #121 from jlaura/120

Adds valid height check for DEM intersection

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.

knoten/csm.py

@@ -161,20 +161,19 @@ def _(dem, image_pt, camera, max_its = 20, tolerance = 0.001):
     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

tests/test_csm.py

+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