Fix issue 951 on accessing DEM raster of DEMInterpolator (#963)

rad-eng-59 · GitHub Enterprise · commit 70ff9d17672d · 2022-03-23T09:35:26.000-07:00
- Add missing pybind11 Eigen header to pybind script "DEMInterpolator.cpp".

- Add an overloaded constructor to pybind "DEMInterpolator.cpp" to
  construct DEM directly from a raster object.

- Augment the respective DEM python test suite "dem.py" to test
  new feature.

- Add a very short DEM raster file over Himalayas for testing DEM.

- Add a description for DEM raster file in "tests/data/README.md"

- Rename alias "DI" to "DEMInterp" in pybind "DEMInterpolator.cpp"

- Replace raw pointers with smart ones for members in class "DEMInterpolator"

- Replace raw pointer members "_proj" and "_interp" by shared smart ones

- Get rid of dynamic pointers for return values in the pybind code and
  instead build DEMInterpolator object on the stack

- Run clang-format only on pybind script
diff --git a/cxx/isce3/geometry/DEMInterpolator.cpp b/cxx/isce3/geometry/DEMInterpolator.cpp
@@ -10,20 +10,10 @@
 #include <isce3/core/Projections.h>
 #include <isce3/io/Raster.h>
 
-isce3::geometry::DEMInterpolator::
-~DEMInterpolator() {
-    if (_interp) {
-        delete _interp;
-    }
-    if (_proj) {
-        delete _proj;
-    }
-}
-
 /** Set EPSG code for input DEM */
 void isce3::geometry::DEMInterpolator::epsgCode(int epsgcode) {
     _epsgcode = epsgcode;
-    _proj = isce3::core::createProj(epsgcode);
+    _proj = isce3::core::makeProjection(epsgcode);
 }
 
 // Load DEM subset into memory
@@ -58,7 +48,7 @@ isce3::error::ErrorCode isce3::geometry::DEMInterpolator::loadDEM(
     //Initialize projection
     int epsgcode = demRaster.getEPSG();
     _epsgcode = epsgcode;
-    _proj = isce3::core::createProj(epsgcode);
+    _proj = isce3::core::makeProjection(epsgcode);
 
     // Validate requested geographic bounds with input DEM raster
     if (minX < firstX) {
@@ -112,7 +102,7 @@ isce3::error::ErrorCode isce3::geometry::DEMInterpolator::loadDEM(
     demRaster.getBlock(_dem.data(), xstart, ystart, width, length);
 
     // Initialize internal interpolator
-    _interp = isce3::core::createInterpolator<float>(_interpMethod);
+    _interp = std::unique_ptr<isce3::core::Interpolator<float>>(isce3::core::createInterpolator<float>(_interpMethod));
 
     // Indicate we have loaded a valid raster
     _haveRaster = true;
@@ -144,7 +134,7 @@ loadDEM(isce3::io::Raster & demRaster) {
     //Initialize projection
     int epsgcode = demRaster.getEPSG();
     _epsgcode = epsgcode;
-    _proj = isce3::core::createProj(epsgcode);
+    _proj = isce3::core::makeProjection(epsgcode);
 
     // Store actual starting lat/lon for raster subset
     _xstart = firstX;
@@ -159,7 +149,7 @@ loadDEM(isce3::io::Raster & demRaster) {
     demRaster.getBlock(_dem.data(), 0, 0, width, length);
 
     // Initialize internal interpolator
-    _interp = isce3::core::createInterpolator<float>(_interpMethod);
+    _interp = std::unique_ptr<isce3::core::Interpolator<float>>(isce3::core::createInterpolator<float>(_interpMethod));
 
     // Indicate we have loaded a valid raster
     _haveRaster = true;
diff --git a/cxx/isce3/geometry/DEMInterpolator.h b/cxx/isce3/geometry/DEMInterpolator.h
@@ -6,6 +6,7 @@
 
 #pragma once
 
+#include <memory>
 #include "forward.h"
 
 // pyre
@@ -52,9 +53,7 @@ class isce3::geometry::DEMInterpolator {
             _maxValue{height},
             _epsgcode{epsg},
             _interpMethod{method} {}
-
-        /** Destructor */
-        ~DEMInterpolator();
+       
 
         /** Read in subset of data from a DEM with a supported projection */
         isce3::error::ErrorCode loadDEM(isce3::io::Raster& demRaster,
@@ -138,7 +137,7 @@ class isce3::geometry::DEMInterpolator {
         void epsgCode(int epsgcode);
 
         /** Get Pointer to a ProjectionBase */
-        inline isce3::core::ProjectionBase* proj() const { return _proj; }
+        inline isce3::core::ProjectionBase* proj() const {return _proj.get(); }
 
         /** Get interpolator method enum */
         inline isce3::core::dataInterpMethod interpMethod() const {
@@ -159,10 +158,10 @@ class isce3::geometry::DEMInterpolator {
         float _maxValue;
         // Pointer to a ProjectionBase
         int _epsgcode;
-        isce3::core::ProjectionBase * _proj = nullptr;
+        std::shared_ptr<isce3::core::ProjectionBase> _proj;
         // Pointer to an Interpolator
         isce3::core::dataInterpMethod _interpMethod;
-        isce3::core::Interpolator<float> * _interp = nullptr;
+        std::shared_ptr<isce3::core::Interpolator<float>> _interp;
         // 2D array for storing DEM subset
         isce3::core::Matrix<float> _dem;
         // Starting x/y for DEM subset and spacing
diff --git a/python/extensions/pybind_isce3/geometry/DEMInterpolator.cpp b/python/extensions/pybind_isce3/geometry/DEMInterpolator.cpp
@@ -1,68 +1,94 @@
 #include "DEMInterpolator.h"
 
-#include <isce3/core/Constants.h>
-#include <isce3/io/Raster.h>
-#include <Eigen/Dense>
+#include <memory>
+#include <pybind11/eigen.h>
 #include <stdexcept>
 #include <string>
 
+#include <Eigen/Dense>
+
+#include <isce3/core/Constants.h>
+#include <isce3/io/Raster.h>
+
 namespace py = pybind11;
 
-using DI = isce3::geometry::DEMInterpolator;
+using DEMInterp = isce3::geometry::DEMInterpolator;
 
-void addbinding(pybind11::class_<DI> & pyDEMInterpolator)
+void addbinding(pybind11::class_<DEMInterp>& pyDEMInterpolator)
 {
     pyDEMInterpolator
-        .def(py::init<double, isce3::core::dataInterpMethod, int>(),
-            py::arg("height") = 0.0,
-            py::arg("method") = isce3::core::BILINEAR_METHOD,
-            py::arg("epsg") = 4326)
-        // For convenience allow a string, too.
-        .def(py::init([](double h, const std::string & method, int epsg) {
+            .def(py::init<double, isce3::core::dataInterpMethod, int>(),
+                    py::arg("height") = 0.0,
+                    py::arg("method") = isce3::core::BILINEAR_METHOD,
+                    py::arg("epsg") = 4326)
+            // For convenience allow a string, too.
+            .def(py::init([](double h, const std::string& method, int epsg) {
                 auto m = parseDataInterpMethod(method);
-                return new DI(h, m, epsg);
+                return DEMInterp(h, m, epsg);
+            }),
+                    py::arg("height") = 0.0, py::arg("method") = "bilinear",
+                    py::arg("epsg") = 4326)
+
+            // This constructor is similar to method "loadDEM" but is more
+            // convenient!
+            .def(py::init([](isce3::io::Raster& raster_obj) {
+                DEMInterp dem {};
+                dem.loadDEM(raster_obj);
+                return dem;
             }),
-            py::arg("height") = 0.0,
-            py::arg("method") = "bilinear",
-            py::arg("epsg") = 4326)
+                    "Construct DEM from ISCE3 Raster object",
+                    py::arg("raster_obj"))
 
-        .def("load_dem",
-            py::overload_cast<isce3::io::Raster&>(&DI::loadDEM))
-        .def("load_dem",
-            py::overload_cast<isce3::io::Raster&, double, double, double, double>
-                (&DI::loadDEM),
-            py::arg("raster"), py::arg("min_x"), py::arg("max_x"),
-                py::arg("min_y"), py::arg("max_y"))
+            .def("load_dem",
+                    py::overload_cast<isce3::io::Raster&>(&DEMInterp::loadDEM))
+            .def("load_dem",
+                    py::overload_cast<isce3::io::Raster&, double, double,
+                            double, double>(&DEMInterp::loadDEM),
+                    py::arg("raster"), py::arg("min_x"), py::arg("max_x"),
+                    py::arg("min_y"), py::arg("max_y"))
 
-        .def("interpolate_lonlat", &DI::interpolateLonLat)
-        .def("interpolate_xy", &DI::interpolateXY)
+            .def("interpolate_lonlat", &DEMInterp::interpolateLonLat)
+            .def("interpolate_xy", &DEMInterp::interpolateXY)
 
-        .def_property("ref_height",
-            py::overload_cast<>(&DI::refHeight, py::const_),
-            py::overload_cast<double>(&DI::refHeight))
-        .def_property_readonly("have_raster", &DI::haveRaster)
-        .def_property("interp_method",
-            py::overload_cast<>(&DI::interpMethod, py::const_),
-            py::overload_cast<isce3::core::dataInterpMethod>(&DI::interpMethod))
+            .def_property("ref_height",
+                    py::overload_cast<>(&DEMInterp::refHeight, py::const_),
+                    py::overload_cast<double>(&DEMInterp::refHeight))
+            .def_property_readonly("have_raster", &DEMInterp::haveRaster)
+            .def_property("interp_method",
+                    py::overload_cast<>(&DEMInterp::interpMethod, py::const_),
+                    py::overload_cast<isce3::core::dataInterpMethod>(
+                            &DEMInterp::interpMethod))
 
-        // Define all these as readonly even though writable in C++ API.
-        // Probably better to just convert your data to a GDAL format than try
-        // to build a DEM on the fly.
-        .def_property_readonly("data", [](DI & self) { // .data() isn't const
-            if (!self.haveRaster()) {
-                throw std::out_of_range("Tried to access DEM data but size=0");
-            }
-            using namespace Eigen;
-            using MatF = Eigen::Matrix<float, Dynamic, Dynamic, RowMajor>;
-            Map<const MatF> mat(self.data(), self.length(), self.width());
-            return mat;
-        }, py::return_value_policy::reference_internal)
-        .def_property_readonly("x_start", py::overload_cast<>(&DI::xStart, py::const_))
-        .def_property_readonly("y_start", py::overload_cast<>(&DI::yStart, py::const_))
-        .def_property_readonly("delta_x", py::overload_cast<>(&DI::deltaX, py::const_))
-        .def_property_readonly("delta_y", py::overload_cast<>(&DI::deltaY, py::const_))
-        .def_property_readonly("width", py::overload_cast<>(&DI::width, py::const_))
-        .def_property_readonly("length", py::overload_cast<>(&DI::length, py::const_))
-        .def_property_readonly("epsg_code", py::overload_cast<>(&DI::epsgCode, py::const_))
-        ;
+            // Define all these as readonly even though writable in C++ API.
+            // Probably better to just convert your data to a GDAL format than
+            // try to build a DEM on the fly.
+            .def_property_readonly(
+                    "data",
+                    [](DEMInterp& self) { // .data() isn't const
+                        if (!self.haveRaster()) {
+                            throw std::out_of_range(
+                                    "Tried to access DEM data but size=0");
+                        }
+                        using namespace Eigen;
+                        using MatF = Eigen::Matrix<float, Dynamic, Dynamic,
+                                RowMajor>;
+                        Map<const MatF> mat(
+                                self.data(), self.length(), self.width());
+                        return mat;
+                    },
+                    py::return_value_policy::reference_internal)
+            .def_property_readonly("x_start",
+                    py::overload_cast<>(&DEMInterp::xStart, py::const_))
+            .def_property_readonly("y_start",
+                    py::overload_cast<>(&DEMInterp::yStart, py::const_))
+            .def_property_readonly("delta_x",
+                    py::overload_cast<>(&DEMInterp::deltaX, py::const_))
+            .def_property_readonly("delta_y",
+                    py::overload_cast<>(&DEMInterp::deltaY, py::const_))
+            .def_property_readonly(
+                    "width", py::overload_cast<>(&DEMInterp::width, py::const_))
+            .def_property_readonly("length",
+                    py::overload_cast<>(&DEMInterp::length, py::const_))
+            .def_property_readonly("epsg_code",
+                    py::overload_cast<>(&DEMInterp::epsgCode, py::const_));
 }
diff --git a/tests/data/README.md b/tests/data/README.md
@@ -117,4 +117,13 @@ and then post-processed and converted into L0B product via its python utility to
   The Tx range lines types are of HPA, LNA, and BYPASS. BYPASS range line interval is 20.
 
 [1]: M. Shimada et al., "PALSAR Radiometric and Geometric Calibration",
-*IEEE Trans. Geosci. Remote Sens.*, pp. 3915-3932, December 2009.     
+*IEEE Trans. Geosci. Remote Sens.*, pp. 3915-3932, December 2009.
+
+## DEM
+- **dem_himalayas_E81p5_N28p3_short.tiff**
+
+  DEM raster file over a very small area of Himalayas downloaded from NISAR-DEM AWS S3
+  bucket as follows:
+  ```
+  $ stage_dem.py -b 81.45 28.29 81.5 28.3 -o dem_himalayas_E81p5_N28p3_short.vrt
+  ```
diff --git a/tests/data/dem_himalayas_E81p5_N28p3_short.tiff b/tests/data/dem_himalayas_E81p5_N28p3_short.tiff
diff --git a/tests/python/extensions/pybind/geometry/dem.py b/tests/python/extensions/pybind/geometry/dem.py
@@ -1,42 +1,79 @@
 #!/usr/bin/env python3
 import pytest
-import pybind_isce3.geometry as m
+from isce3.geometry import DEMInterpolator
+from isce3.io import Raster
 from pybind_isce3.core import DataInterpMethod
+import iscetest
 
-def test_const():
+import os
+import collections as cl
+import numpy.testing as npt
+
+from osgeo import gdal
+
+
+def dem_info_from_gdal(file_raster: str) -> cl.namedtuple:
+    """Get shape, min, max, mean of dem"""
+    dset = gdal.Open(file_raster, gdal.GA_ReadOnly)
+    band = dset.GetRasterBand(1)
+    dem = band.ReadAsArray()
+    return cl.namedtuple('dem_info', 'shape min max mean')(
+        dem.shape, dem.min(), dem.max(), dem.mean())
+
+
+def test_constructor_ref_height():
     href = 10.
 
-    dem = m.DEMInterpolator()
+    dem = DEMInterpolator()
     dem.ref_height = href
     assert dem.ref_height == href
 
-    dem = m.DEMInterpolator(href)
+    dem = DEMInterpolator(href)
     assert dem.ref_height == href
 
     assert dem.interpolate_xy(0, 0) == href
     assert dem.interpolate_lonlat(0, 0) == href
 
-    assert dem.have_raster == False
+    npt.assert_equal(dem.have_raster, False)
+
+
+def test_constructor_raster_obj():
+    # filename of the DEM ratster
+    filename_dem = 'dem_himalayas_E81p5_N28p3_short.tiff'
+    file_raster = os.path.join(iscetest.data, filename_dem)
+    # get some DEM info via gdal to be used as a reference for V&V
+    dem_info = dem_info_from_gdal(file_raster)
+    # build DEM object
+    raster_obj = Raster(file_raster)
+    dem_obj = DEMInterpolator(raster_obj)
+    # validate existence and details of DEM data
+    npt.assert_equal(dem_obj.have_raster, True, err_msg='No DEM ratser data')
+    npt.assert_equal(dem_obj.data.shape, dem_info.shape,
+                     err_msg='Wrong shape of DEM ratser data')
+    npt.assert_allclose(dem_obj.data.min(), dem_info.min,
+                        err_msg='Wrong min DEM height')
+    npt.assert_allclose(dem_obj.data.max(), dem_info.max,
+                        err_msg='Wrong max DEM height')
+    npt.assert_allclose(dem_obj.data.mean(), dem_info.mean,
+                        err_msg='Wrong mean DEM height')
 
 
 def test_methods():
     # pybind11::enum_ is not iterable
     for name in "SINC BILINEAR BICUBIC NEAREST BIQUINTIC".split():
         # enum constructor
         method = getattr(DataInterpMethod, name)
-        dem = m.DEMInterpolator(method=method)
+        dem = DEMInterpolator(method=method)
         assert dem.interp_method == method
         # string constructor
-        dem = m.DEMInterpolator(method=name)
+        dem = DEMInterpolator(method=name)
         assert dem.interp_method == method
 
-    dem = m.DEMInterpolator(method="bicubic")
+    dem = DEMInterpolator(method="bicubic")
     assert dem.interp_method == DataInterpMethod.BICUBIC
 
-    dem = m.DEMInterpolator(method="biCUBic")
+    dem = DEMInterpolator(method="biCUBic")
     assert dem.interp_method == DataInterpMethod.BICUBIC
 
     with pytest.raises(ValueError):
-        dem = m.DEMInterpolator(method="TigerKing")
-
-    # TODO Test other methods once we have isce::io::Raster bindings.
+        dem = DEMInterpolator(method="TigerKing")
