diff --git a/idefix/__init__.py b/idefix/__init__.py index 96b59e1..7266644 100644 --- a/idefix/__init__.py +++ b/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 diff --git a/idefix/helpers.py b/idefix/helpers.py new file mode 100644 index 0000000..ab22eaf --- /dev/null +++ b/idefix/helpers.py @@ -0,0 +1,89 @@ +#!/usr/bin/env python +# file helpers.py +# author Florent Guiotte +# 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 diff --git a/idefix/utils.py b/idefix/utils.py index 3334c7c..6fde6e7 100644 --- a/idefix/utils.py +++ b/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 diff --git a/idefix/vxl.py b/idefix/vxl.py index 3f2d803..fcee31a 100644 --- a/idefix/vxl.py +++ b/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: diff --git a/setup.py b/setup.py index 79785af..fbc8e19 100644 --- a/setup.py +++ b/setup.py @@ -36,6 +36,9 @@ setuptools.setup( 'tqdm', 'matplotlib', 'pathlib', - 'rasterio' + 'rasterio', + 'laspy', + 'humanize', + #'mayavi', Optional, for vxl.plot() ], ) diff --git a/test/test_helpers.py b/test/test_helpers.py new file mode 100644 index 0000000..378f7bf --- /dev/null +++ b/test/test_helpers.py @@ -0,0 +1,41 @@ +#!/usr/bin/env python +# file test_helpers.py +# author Florent Guiotte +# 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' diff --git a/test/test_helpers/pc0.txt b/test/test_helpers/pc0.txt new file mode 100644 index 0000000..b5d64b8 --- /dev/null +++ b/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 diff --git a/test/test_utils.py b/test/test_utils.py index 472f53a..566c53f 100644 --- a/test/test_utils.py +++ b/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) -