Merge pull request #125 from jlaura/libload

Adds surface

.github/workflows/ci.yml

@@ -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

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

environment.yml

@@ -19,9 +19,10 @@ dependencies:
   - plotly-orca 
   - psutil 
   - pvl
+  - pytest
+  - pytest-xdist
   - pyproj
   - kalasiris
-  - pytest
   - python>=3
   - requests
   - scipy

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))
+        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

knoten/surface.py

+"""
+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

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:

tests/test_csm.py

+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')

tests/test_surface.py

+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