diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 02573ee85dc4fd98ff3730099438c9a815370f72..ea1f0f53eaed04abd8c9db70cba19ef833a56659 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -32,7 +32,7 @@ jobs:
environment-file: environment.yml
auto-activate-base: false
auto-update-conda: true
- python-version: ${{ matrix.python-version }}
+ python-version: ${{ matrix.python-version }}
- name: Check build environment
run: |
conda list
@@ -41,4 +41,4 @@ jobs:
python setup.py install
- name: Test Python Package
run: |
- pytest
\ No newline at end of file
+ pytest -n4
\ No newline at end of file
diff --git a/CHANGELOG.md b/CHANGELOG.md
index cb73010f4a27e46c785db94082a16898b34a4a23..51cb0dac8185c29662825119f2eb9335e7fa68a9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -36,6 +36,8 @@ release.
## Unreleased
### Added
+- `generate_image_coordinate` to `csm.py`. This provides a similar interface to `generate_ground_coordinate` and abstracts away the `csmapi` from the user.
+- A surface class (moved from AutoCNet; credit @jessemapel) with support for Ellipsoid DEMs and basic support for raster DEMs readable by the plio.io.io_gdal.GeoDataset. Support is basic because it uses a single pixel intersection and not an interpolated elevation like ISIS does.
- 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. Fixes [#120](https://github.com/DOI-USGS/knoten/issues/120)
### Changed
diff --git a/environment.yml b/environment.yml
index 938bf87af05263f99a19c529685d89d812064ab8..792cdf89b07b606594199553573e298dbccc8b20 100644
--- a/environment.yml
+++ b/environment.yml
@@ -19,9 +19,10 @@ dependencies:
- plotly-orca
- psutil
- pvl
+ - pytest
+ - pytest-xdist
- pyproj
- kalasiris
- - pytest
- python>=3
- requests
- scipy
diff --git a/knoten/csm.py b/knoten/csm.py
index f80753052a4b4d9c8ff55b02efd29a3e3102b348..da04ae9596101ad51ba9f5e30dabf51cf4fbfc43 100644
--- a/knoten/csm.py
+++ b/knoten/csm.py
@@ -10,7 +10,7 @@ import requests
import scipy.stats
from functools import singledispatch
-from plio.io.io_gdal import GeoDataset
+from knoten.surface import EllipsoidDem
from knoten import utils
class NumpyEncoder(json.JSONEncoder):
@@ -127,7 +127,7 @@ def generate_ground_point(dem, image_pt, camera):
GeoDataset object generated from Plio off of a Digital Elevation
Model (DEM)
image_pt : tuple
- Pair of x, y coordinates in pixel space
+ Pair of x, y (sample, line) coordinates in pixel space
camera : object
USGSCSM camera model object
max_its : int, optional
@@ -145,11 +145,10 @@ def generate_ground_point(dem, image_pt, camera):
if not isinstance(image_pt, csmapi.ImageCoord):
# Support a call where image_pt is in the form (x,y)
image_pt = csmapi.ImageCoord(*image_pt)
-
return camera.imageToGround(image_pt, dem)
-@generate_ground_point.register(GeoDataset)
-def _(dem, image_pt, camera, max_its = 20, tolerance = 0.001):
+@generate_ground_point.register
+def _(dem: EllipsoidDem, image_pt, camera, max_its = 20, tolerance = 0.0001, dem_type='radius'):
if not isinstance(image_pt, csmapi.ImageCoord):
# Support a call where image_pt is in the form (x,y)
image_pt = csmapi.ImageCoord(*image_pt)
@@ -166,21 +165,23 @@ def _(dem, image_pt, camera, max_its = 20, tolerance = 0.001):
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:
+ height = dem.get_height(lat, lon)
+ if height is None:
raise ValueError(f'No DEM height at {lat}, {lon}')
-
- next_intersection = camera.imageToGround(image_pt, float(height))
- dist = _compute_intersection_distance(intersection, next_intersection)
+ next_intersection = generate_ground_point(float(height), image_pt, camera)
+ dist = _compute_intersection_distance(intersection, next_intersection)
intersection = next_intersection
iterations += 1
if dist < tolerance:
break
return intersection
+def generate_image_coordinate(ground_pt, camera):
+ if not isinstance(ground_pt, csmapi.EcefCoord):
+ ground_pt = csmapi.EcefCoord(*ground_pt)
+ return camera.groundToImage(ground_pt)
+
def _compute_intersection_distance(intersection, next_intersection):
"""
Private func that takes two csmapi Ecef objects or other objects with
diff --git a/knoten/surface.py b/knoten/surface.py
new file mode 100644
index 0000000000000000000000000000000000000000..2b85452b3de3fb4281f9a9f720e6680a91ec7aae
--- /dev/null
+++ b/knoten/surface.py
@@ -0,0 +1,145 @@
+"""
+A set of classes that represent the target surface. Each class implements the
+get_height and get_radius functions for computing the height and radius respectively
+at a given ground location (geocentric latitude and longitude).
+"""
+
+import numpy as np
+from plio.io.io_gdal import GeoDataset
+
+class EllipsoidDem:
+ """
+ A biaxial ellipsoid surface model.
+ """
+
+ def __init__(self, semi_major, semi_minor = None):
+ """
+ Create an ellipsoid DEM 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
+
+ def get_height(self, lat, lon):
+ """
+ Get the height above the ellipsoid at a ground location
+
+ Parameters
+ ----------
+ lat : float
+ The geocentric latitude in degrees
+ lon : float
+ The longitude in degrees
+ """
+ return 0
+
+ def get_radius(self, lat, lon):
+ """
+ Get the radius at a ground location
+
+ Parameters
+ ----------
+ lat : float
+ The geocentric latitude in degrees
+ lon : float
+ The longitude in degrees
+ """
+ cos_lon = np.cos(np.deg2rad(lon))
+ sin_lon = np.sin(np.deg2rad(lon))
+ cos_lat = np.cos(np.deg2rad(lat))
+ sin_lat = np.sin(np.deg2rad(lat))
+
+ denom = self.b * self.b * cos_lon * cos_lon
+ denom += self.a * self.a * sin_lon * sin_lon
+ denom *= self.c * self.c * cos_lat * cos_lat
+ denom += self.a * self.a * self.b * self.b * sin_lat * sin_lat
+ radius = (self.a * self.b * self.c) / np.sqrt(denom)
+ return radius
+
+class GdalDem(EllipsoidDem):
+ """
+ A raster DEM surface model.
+ """
+
+ def __init__(self, dem, semi_major, semi_minor = None, dem_type=None):
+ """
+ Create a GDAL dem from a dem file
+
+ Parameters
+ ----------
+ dem : str
+ The DEM file
+ semi_major : float
+ The equatorial semi-major radius of the reference ellipsoid.
+ semi_minor : float
+ The polar semi-minor radius of the reference ellipsoid.
+ dem_type : str
+ The type of DEM, either height above reference ellipsoid or radius.
+ """
+ super().__init__(semi_major, semi_minor)
+ dem_types = ('height', 'radius')
+ if dem_type is None:
+ dem_type = dem_types[0]
+ if dem_type not in dem_types:
+ raise ValueError(f'DEM type {dem_type} is not a valid option.')
+ self.dem = GeoDataset(dem)
+ self.dem_type = dem_type
+
+ def get_raster_value(self, lat, lon):
+ """
+ Get the value of the dem raster at a ground location
+
+ Parameters
+ ----------
+ lat : float
+ The geocentric latitude in degrees
+ lon : float
+ The longitude in degrees
+ """
+ px, py = self.dem.latlon_to_pixel(lat, lon)
+ value = self.dem.read_array(1, [px, py, 1, 1])[0][0]
+ if value == self.dem.no_data_value:
+ return None
+ return value
+
+ def get_height(self, lat, lon):
+ """
+ Get the height above the ellipsoid at a ground location
+
+ Parameters
+ ----------
+ lat : float
+ The geocentric latitude in degrees
+ lon : float
+ The longitude in degrees
+ """
+ height = self.get_raster_value(lat, lon)
+ if self.dem_type == 'radius' and height is not None:
+ height -= super().get_radius(lat, lon)
+ return height
+
+ def get_radius(self, lat, lon):
+ """
+ Get the radius at a ground location
+
+ Parameters
+ ----------
+ lat : float
+ The geocentric latitude in degrees
+ lon : float
+ The longitude in degrees
+ """
+ radius = self.get_raster_value(lat, lon)
+ if self.dem_type == 'height':
+ radius += super().get_radius(lat, lon)
+ return radius
diff --git a/recipe/meta.yaml b/recipe/meta.yaml
index 215af236191ed6de50acbfc4c74c47993dd5393f..50d5571872f77561ab4e1db558257debe843212f 100644
--- a/recipe/meta.yaml
+++ b/recipe/meta.yaml
@@ -23,7 +23,7 @@ requirements:
- python
- setuptools
run:
- - ale>=0.3.2
+ - ale>=0.10.0
- python
- csmapi>=1.0
- gdal>=3.0.0
@@ -34,7 +34,7 @@ requirements:
- requests
- scipy
- shapely
- - usgscsm>=1.3.1
+ - usgscsm>=2.0
test:
imports:
diff --git a/tests/test_csm.py b/tests/test_csm.py
index 94fb524ec5daeb83b5b779625c31e40932892f3a..d10320788fca4bb8e7b6985b49eea95e02d6b534 100644
--- a/tests/test_csm.py
+++ b/tests/test_csm.py
@@ -1,18 +1,21 @@
+from collections import namedtuple
from unittest import mock
import pytest
from plio.io.io_gdal import GeoDataset
import csmapi
-from knoten import csm
+from knoten import csm, surface
@pytest.fixture
def mock_dem():
- mock_dem = mock.MagicMock(spec_set=GeoDataset)
+ mock_surface = mock.MagicMock(spec=surface.GdalDem)
+ mock_dem = mock.MagicMock(spec=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
+ mock_surface.dem = mock_dem
+ return mock_surface
@pytest.fixture
def mock_sensor():
@@ -39,7 +42,7 @@ def test_generate_ground_point_with_imagecoord(mock_sensor, pt):
@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_compute_intersection,
mock_pyproj_transformer,
mock_sensor, pt, mock_dem):
csm.generate_ground_point(mock_dem, pt, mock_sensor)
@@ -48,16 +51,15 @@ def test_generate_ground_point_with_dtm(mock_get_radii,
# 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
+ mock_dem.get_height.return_value = None
with pytest.raises(ValueError):
- csm.generate_ground_point(mock_dem, pt, mock_sensor)
+ res = csm.generate_ground_point(mock_dem, pt, mock_sensor)
def test__compute_intersection_distance():
Point = namedtuple("Point", 'x, y, z')
diff --git a/tests/test_surface.py b/tests/test_surface.py
new file mode 100644
index 0000000000000000000000000000000000000000..dcc733a18fd8973bb57f517ceedb3256dc70237e
--- /dev/null
+++ b/tests/test_surface.py
@@ -0,0 +1,68 @@
+import unittest
+
+from unittest import mock
+
+from knoten import surface
+
+class TestEllipsoidDem(unittest.TestCase):
+
+ def test_height(self):
+ test_dem = surface.EllipsoidDem(3396190, 3376200)
+ self.assertEqual(test_dem.get_height(0, 0), 0)
+ self.assertEqual(test_dem.get_height(0, 180), 0)
+ self.assertEqual(test_dem.get_height(90, 100), 0)
+
+ def test_radius(self):
+ test_dem = surface.EllipsoidDem(3396190, 3376200)
+ self.assertEqual(test_dem.get_radius(0, 0), 3396190)
+ self.assertEqual(test_dem.get_radius(0, 180), 3396190)
+ self.assertEqual(test_dem.get_radius(90, 300), 3376200)
+
+ def tearDown(self):
+ pass
+
+
+class TestGdalDem(unittest.TestCase):
+
+ def test_height(self):
+ with mock.patch('knoten.surface.GeoDataset') as mockDataset:
+ mockInstance = mockDataset.return_value
+ mockInstance.latlon_to_pixel.return_value = (1,2)
+ mockInstance.read_array.return_value = [[100]]
+ test_dem = surface.GdalDem('TestDem.cub', 3396190, 3376200)
+ self.assertEqual(test_dem.get_height(0, 0), 100)
+ self.assertEqual(test_dem.get_height(0, 180), 100)
+ self.assertEqual(test_dem.get_height(90, 300), 100)
+
+ def test_height_from_radius(self):
+ with mock.patch('knoten.surface.GeoDataset') as mockDataset:
+ mockInstance = mockDataset.return_value
+ mockInstance.latlon_to_pixel.return_value = (1,2)
+ mockInstance.read_array.return_value = [[3396190]]
+ test_dem = surface.GdalDem('TestDem.cub', 3396190, 3376200, 'radius')
+ self.assertEqual(test_dem.get_height(0, 0), 0)
+ self.assertEqual(test_dem.get_height(0, 180), 0)
+ self.assertEqual(test_dem.get_height(90, 300), 19990)
+
+ def test_radius(self):
+ with mock.patch('knoten.surface.GeoDataset') as mockDataset:
+ mockInstance = mockDataset.return_value
+ mockInstance.latlon_to_pixel.return_value = (1,2)
+ mockInstance.read_array.return_value = [[3396190]]
+ test_dem = surface.GdalDem('TestDem.cub', 3396190, 3376200, 'radius')
+ self.assertEqual(test_dem.get_radius(0, 0), 3396190)
+ self.assertEqual(test_dem.get_radius(0, 180), 3396190)
+ self.assertEqual(test_dem.get_radius(90, 300), 3396190)
+
+ def test_radius_from_height(self):
+ with mock.patch('knoten.surface.GeoDataset') as mockDataset:
+ mockInstance = mockDataset.return_value
+ mockInstance.latlon_to_pixel.return_value = (1,2)
+ mockInstance.read_array.return_value = [[100]]
+ test_dem = surface.GdalDem(mockInstance, 3396190, 3376200)
+ self.assertEqual(test_dem.get_radius(0, 0), 3396290)
+ self.assertEqual(test_dem.get_radius(0, 180), 3396290)
+ self.assertEqual(test_dem.get_radius(90, 300), 3376300)
+
+ def tearDown(self):
+ pass