Add helpers module with dsm() and interpolate()

idefix/__init__.py

@@ -9,6 +9,6 @@ Utils and production pipelines for processing LiDAR point clouds.
 
 """
 
-__all__ = ['utils', 'io', 'vxl']
+__all__ = ['utils', 'io', 'vxl', 'helpers']
 
-from . import utils, io, vxl
+from . import utils, io, vxl, helpers

idefix/helpers.py

@@ -0,0 +1,89 @@
+#!/usr/bin/env python
+# file helpers.py
+# author Florent Guiotte <florent.guiotte@irisa.fr>
+# version 0.0
+# date 24 août 2020
+"""High-level helper functions.
+
+This module contains high-level helper functions. This module shows many
+exemple on the use of idefix package and other packages (sap, rasterio,
+...) to process point clouds.
+
+"""
+
+import numpy as np
+from scipy.interpolate import griddata
+from rasterio import fill
+import sap
+
+from .vxl import get_grid, bin, squash
+
+def interpolate(raster, method='linear'):
+    """Interpolate masked raster.
+
+    Wrapper function to interpolate missing values in masked raster.
+    The 'linear', 'nearest' and 'cubic' implementation are from `Scipy`_
+    while the 'idw' (inverse distance weighting) is provided by
+    `rasterio`_.
+
+    .. _Scipy: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html
+    .. _rasterio: https://rasterio.readthedocs.io/en/latest/api/rasterio.fill.html
+
+    Parameters
+    ----------
+    raster : masked ndarray
+        The raster with missing values masked.
+    method : str
+        Can be 'linear', 'nearest', 'cubic' or 'idw'.
+
+    Returns
+    -------
+    out : ndarray
+        The raster with filled missing values.
+
+    """
+    if method == 'idw':
+        raster = fill.fillnodata(raster)
+    else:
+        coords = np.argwhere(~raster.mask)
+        values = raster.compressed()
+        grid   = np.argwhere(raster.mask)
+
+        raster[raster.mask] = griddata(coords, values, grid, method=method)
+
+        if method != 'nearest':
+            raster.mask = np.isnan(raster)
+            raster = interpolate(raster, 'nearest')
+
+        raster = np.array(raster)
+
+    assert not np.isnan(raster).any()
+
+    return raster
+
+def dsm(pcloud, cell_size=1., last=False):
+    """Create the digital surface model (DSM) of the point cloud.
+
+    Parameters
+    ----------
+    pcloud : recarray
+        A point cloud loaded with :mod:`idefix.io`.
+    cell_size : scalar
+        The size of the cells in meter. Cells are square. Default is 1
+        meter.
+    last : bool
+        Specifies whether the first echo (`False`) or the last echo
+        (`True`) should be taken into account. Default is `False`.
+
+    Returns
+    -------
+    dsm : ndarray
+        The DSM of the point cloud.
+
+    """
+    grid = get_grid(pcloud.spatial, cell_size)
+    vxlg = bin(grid, pcloud.spatial, pcloud.spatial[:,2], 'mean')
+    rstr = squash(vxlg, 'bottom' if last else 'top')
+    rstr = interpolate(rstr, 'idw')
+
+    return rstr

idefix/utils.py

@@ -10,14 +10,12 @@ This module contains common utils for basic point cloud management and dataviz.
 Notes
 -----
 
-Everything should be highly tested there.
+Everything is well tested there.
 
 """
 
 import numpy as np
 import logging
-from scipy.interpolate import griddata
-from rasterio import fill
 
 log = logging.getLogger(__name__)
 
@@ -59,47 +57,3 @@ def bbox(data):
          [xmax, ymax, zmax]]
     """
     return np.array((np.min(data, axis=0), np.max(data, axis=0)))
-
-
-def interpolate(raster, method='linear'):
-    """Interpolate masked raster.
-
-    Wrapper function to interpolate missing values in masked raster.
-    The 'linear', 'nearest' and 'cubic' implementation are from `Scipy`_
-    while the 'idw' (inverse distance weighting) is provided by
-    `rasterio`_.
-
-    .. _Scipy: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html
-    .. _rasterio: https://rasterio.readthedocs.io/en/latest/api/rasterio.fill.html
-
-    Parameters
-    ----------
-    raster : masked ndarray
-        The raster with missing values masked.
-    method : str
-        Can be 'linear', 'nearest', 'cubic' or 'idw'.
-
-    Returns
-    -------
-    out : ndarray
-        The raster with filled missing values.
-
-    """
-    if method == 'idw':
-        raster = fill.fillnodata(raster)
-    else:
-        coords = np.argwhere(~raster.mask)
-        values = raster.compressed()
-        grid   = np.argwhere(raster.mask)
-
-        raster[raster.mask] = griddata(coords, values, grid, method=method)
-
-        if method != 'nearest':
-            raster.mask = np.isnan(raster)
-            raster = interpolate(raster, 'nearest')
-
-        raster = np.array(raster)
-
-    assert not np.isnan(raster).any()
-
-    return raster

idefix/vxl.py

@@ -12,7 +12,6 @@ import logging
 import numpy as np
 import humanize
 from .utils import bbox
-import mayavi.mlab as mlab
 
 log = logging.getLogger(__name__)
 
@@ -91,7 +90,7 @@ def bin(grid, spatial, feature=None, method='density'):
         - 'mean': The mean of feature value in each cell.
         - 'mode': The modal (most common) in each cell. Designed for labels on
           point cloud, can be long with rich spectral data. If there is an
-          equal number of elements, then the smallest is returned.
+          equal number of elements, the smallest is returned.
         The default is 'density'.
 
     Returns
@@ -371,6 +370,8 @@ def plot(voxel_grid, vmin=None, vmax=None):
     >>> mlab.savefig(fname, magnification=4)
     >>> mlab.show()
     """
+    import mayavi.mlab as mlab
+
     points = np.where(~voxel_grid.mask)
 
     if vmin or vmax:

setup.py

@@ -36,6 +36,9 @@ setuptools.setup(
         'tqdm',
         'matplotlib',
         'pathlib',
-        'rasterio'
+        'rasterio',
+        'laspy',
+        'humanize',
+        #'mayavi', Optional, for vxl.plot()
     ],
       )

test/test_helpers.py

@@ -0,0 +1,41 @@
+#!/usr/bin/env python
+# file test_helpers.py
+# author Florent Guiotte <florent.guiotte@irisa.fr>
+# version 0.0
+# date 24 août 2020
+"""Abstract
+
+doc.
+"""
+
+import numpy as np
+import pytest
+from idefix import helpers, io
+
+@pytest.fixture
+def ma_raster():
+    rs = np.random.RandomState(42)
+    raster = rs.random((10,10))
+    raster = np.ma.array(raster, mask=raster<.1)
+    return raster
+
+@pytest.mark.parametrize('method', 
+        ['nearest', 'linear', 'cubic', 'idw'])
+def test_interpolate(ma_raster, method):
+    helpers.interpolate(ma_raster, method)
+
+def _data_pc(datadir, set_id):
+    path = datadir.join('pc{}.txt'.format(set_id))
+    data = io.load_txt(path, 'x y z i'.split())
+    return data
+
+@pytest.mark.parametrize('params', [
+    {},
+    {'cell_size': 2.},
+    {'last': True}])
+def test_dsm(datadir, params):
+    pc = _data_pc(datadir, 0)
+    dsm = helpers.dsm(pc, **params)
+
+    assert dsm is not None, 'Did not return anything...'
+    assert not np.isnan(dsm).any(), 'Some missing values in DSM'

test/test_helpers/pc0.txt

@@ -0,0 +1,8 @@
+# x y z feature
+1 1 1 2
+1 3 2 5
+1 3 2 5
+1 3 2 10
+1 3 2 20
+10 10 10 1
+5 5 5 0

test/test_utils.py

@@ -11,13 +11,6 @@ import numpy as np
 import pytest
 from idefix import utils
 
-@pytest.fixture
-def ma_raster():
-    rs = np.random.RandomState(42)
-    raster = rs.random((10,10))
-    raster = np.ma.array(raster, mask=raster<.1)
-    return raster
-
 @pytest.mark.parametrize("first_input,first_expected", [
     (1, -1),
     (-4, 4),
@@ -37,9 +30,3 @@ def test_bbox(fix_data):
 def test_read(datadir):
     with open(datadir.join('first.txt')) as f:
         assert f.read() == 'hullo\n'
-
-@pytest.mark.parametrize('method', 
-        ['nearest', 'linear', 'cubic', 'idw'])
-def test_interpolate(ma_raster, method):
-    utils.interpolate(ma_raster, method)
-