Merge branch 'develop' into 'master'

Merge v1.5

See merge request Florent/idefix!2

idefix/helpers.py

@@ -14,13 +14,10 @@ exemple on the use of idefix package and other packages (sap, rasterio,
 import numpy as np
 from scipy.interpolate import griddata
 from rasterio import fill
-import rasterio as rio
 import sap
 import higra as hg
-from pathlib import Path
 
-from .vxl import get_grid, bin, squash, fit_bbox
+from .vxl import get_grid, bin, squash
-from .io import load_pc
 
 def interpolate(raster, method='linear'):
     """Interpolate masked raster.
@@ -102,7 +99,7 @@ def dtm_dh_filter(dsm, sigma=.5, epsilon=20000, alpha=2):
     dsm : ndarray
         The DSM.
     sigma : scalar
-        The height theshold to trigger object detection. Default is
+        The height theshold to trigger object detection. Default is 
         0.5 m.
     epsilon : scalar
         The area theshold for ground objects. All objects with surface
@@ -120,12 +117,12 @@ def dtm_dh_filter(dsm, sigma=.5, epsilon=20000, alpha=2):
     """
     mt = sap.MaxTree(dsm)
     area = mt.get_attribute('area')
-
+    
     area_child = hg.accumulate_parallel(mt._tree, area, hg.Accumulators.max)
     pruned = (area - area_child) <= alpha
 
     pruned_tree, pruned_map = hg.simplify_tree(mt._tree, pruned)
-
+    
     dh = mt._alt[pruned_map] - mt._alt[pruned_map][pruned_tree.parents()]
     remove = dh > sigma
 
@@ -133,98 +130,8 @@ def dtm_dh_filter(dsm, sigma=.5, epsilon=20000, alpha=2):
     original_map[pruned_map] = np.arange(pruned_map.size)
     original_map = hg.accumulate_and_max_sequential(mt._tree, original_map, np.arange(mt._tree.num_leaves()), hg.Accumulators.max)
     original_remove = remove[original_map] & (area < epsilon)
-
+    
     dtm = mt.reconstruct(original_remove, filtering='min')
-
+    
     return dtm
 
-def write_raster(raster, bbox, crs, fname):
-    """Write a GeoTiff
-    """
-    west, east = bbox[:,0]
-    south, north = bbox[:,1]
-    width, height = raster.shape
-
-    west, east, south, north, width, height
-
-    transform = rio.transform.from_bounds(west, south, east, north, width, height)
-    dtype = raster.dtype
-
-    with rio.open(fname,
-                     mode='w',
-                     driver='GTiff',
-                     width=raster.shape[0],
-                     height=raster.shape[1],
-                     count=1,
-                     crs=crs,
-                     dtype=dtype,
-                     transform=transform,
-                     compress='lzw') as out_tif:
-        out_tif.write(raster[None])
-
-def rasterize(pc_file, feature='elevation', resolution=1.,
-              bin_structure='voxel', bin_method='mean',
-              squash_method='top', interpolation='idw',
-              out_dir=None, crs=None):
-    """Rasterize a point cloud file.
-
-    Parameters
-    ----------
-    pc_file : Path
-        Path of the point cloud file.
-    feature : str
-        'elevation' or LiDAR feature name contained in the point cloud
-        file.
-    resolution : float
-        The spatial resolution of the raster.
-    bin_structure : str
-        'voxel' or 'pixel'.
-    bin_method : str
-        'mean', 'density' or 'mode'.
-    squash_method : str or tuple of str
-        'top', 'center', 'bottom', 'min', 'max', 'mean', 'median', 'std'
-    interpolation : str
-        'idw', 'nearest', 'linear', 'cubic'
-    out_dir : Path or str or None
-        Path of the directory for the result files. Optional.
-    crs : str or None
-        The coordinate reference system of the point cloud. Required if
-        out_dir is defined to write the GeoTiff.
-
-    Returns
-    -------
-    raster : ndarray
-        Raster.
-    """
-    pc = load_pc(pc_file)
-
-    steps = resolution if bin_structure == 'voxel' else (resolution, resolution, None)
-
-    bbox = fit_bbox(pc.spatial)
-    grid = get_grid(bbox, steps)
-    #ix.vxl.insight(grid, method=bin_method, verbose=True)
-
-    fval = pc.spatial[:,2] if feature == 'elevation' else getattr(pc.feature, feature)
-
-    vxl = bin(grid, pc.spatial, fval, bin_method)
-
-    squash_method = squash_method if not isinstance(squash_method, str) else (squash_method,)
-
-    rasters = []
-    for s in squash_method:
-        raster = squash(vxl, s)
-        raster = interpolate(raster, interpolation)
-        if out_dir:
-            format_dict = {'name': Path(pc_file).stem,
-                           'feature': feature,
-                           'resolution': resolution,
-                           'bin_structure': bin_structure,
-                           'bin_method': bin_method,
-                           'squash_method': s,
-                           'interpolation': interpolation}
-            out_tif_name = Path(out_dir) / \
-            '{name}_{feature}_{resolution}_{bin_structure}_{bin_method}_{squash_method}_{interpolation}.tif'.format(**format_dict)
-            write_raster(raster, bbox, crs, out_tif_name)
-        rasters += [raster]
-
-    return np.array(rasters)

idefix/tools/rasterize.py

@@ -1,101 +0,0 @@
-#!/usr/bin/env python
-# file rasterize-lidar.py
-# author Florent Guiotte <florent.guiotte@irisa.fr>
-# version 0.0
-# date 18 févr. 2021
-"""Create raster of LiDAR features.
-
-Use npz point clouds from Idefix.
-"""
-
-import json
-import argparse
-from multiprocessing import Pool
-from tqdm.auto import tqdm
-from pathlib import Path
-
-import idefix as ix
-
-parser = argparse.ArgumentParser(description='Compute features rasters from .npz point cloud.',
-                                 formatter_class=argparse.RawDescriptionHelpFormatter,
-                                 epilog="""
-The config file can contain any parameters of the
-idefix.helpers.rasterize function in a json file.
-
-You can define 'global' parameters (for all the rasters) and raster
-specific parameters in a list 'rasters'.
-
-See the following `config.json` example file:
-
-{
- "global": {
-  "resolution": 5,
-  "interpolation": "idw",
-  "out_dir": "./results"
- },
- "rasters": [
-  {
-   "feature": "elevation",
-   "bin_structure": "voxel",
-   "bin_method": "mean",
-   "squash_method": ["top", "bottom", "std"]
-  },
-  {
-   "feature": "elevation",
-   "bin_structure": "pixel",
-   "bin_method": "mean"
-  },
-  {
-   "bin_structure": "pixel",
-   "bin_method": "density"
-  },
-  {
-   "feature": "num_returns",
-   "bin_structure": "pixel",
-   "bin_method": "mode"
-  }
- ]
-}
-""")
-
-parser.add_argument('-c', '--config', type=str, help='json file to setup the rasterization processes', dest='c')
-parser.add_argument('-n', '--nprocess', type=int, help='number of child processes to use', default=1, dest='n')
-parser.add_argument('in_dir', type=str, help='the path to the point cloud directory')
-parser.add_argument('out_dir', type=str, help='path to output raster results')
-
-args = parser.parse_args()
-
-def _map_rasterize(kwargs):
-    return ix.helpers.rasterize(**kwargs)
-
-def main():
-    in_dir = Path(args.in_dir)
-    out_dir = Path(args.out_dir)
-    out_dir.mkdir(exist_ok=True)
-
-
-    pc_files = list(in_dir.glob('*.npz'))
-
-
-    config = {'global': {'out_dir': out_dir}}
-    config_json = json.load(open(args.c)) if args.c else {}
-    config['global'].update(config_json['global'] if 'global' in config_json else {})
-    config['rasters'] = config_json['rasters'] if 'rasters' in config_json else {}
-
-    globalc = config['global']
-    rasters = config['rasters'] if 'rasters' in config else [{}]
-    
-    queue = []
-    for pc_file in pc_files:
-        for raster in rasters:
-            job = globalc.copy()
-            job.update(raster)
-            job.update({'pc_file': pc_file})
-            queue += [job]
-
-    pool = Pool(processes=args.n)
-    for _ in tqdm(pool.imap_unordered(_map_rasterize, queue), total=len(queue)):
-        pass
-
-if __name__ == '__main__':
-    main()

idefix/utils.py

@@ -55,42 +55,5 @@ def bbox(data):
 
         [[xmin, ymin, zmin],
          [xmax, ymax, zmax]]
-
-    See Also
-    --------
-    fit_bbox : Return a bounding box fit on fixed coordinates.
     """
     return np.array((np.min(data, axis=0), np.max(data, axis=0)))
-
-def fit_bbox(data, decimals=0):
-    """Return a bounding box fit on fixed coordinates.
-
-    - Round $x$ and $y$ coordinates to match most orthoimagery tiles.
-    - Ceil and floor $z$ coordinates to include all the point on the vertical axis.
-
-    Parameters
-    ----------
-    data : ndarray (N, 3)
-        Bbox or point cloud data of shape (N, 3), i.e. (x,y,z).
-    decimals : int
-        The precision for the rounding, ceiling and flooring operations.
-
-    Returns
-    -------
-    bbox : ndarray
-        Lower and upper points describing the bounding box such as::
-
-        [[xmin, ymin, zmin],
-         [xmax, ymax, zmax]]
-
-    See Also
-    --------
-    bbox : Returns a raw bounding box on the data.
-    """
-    raw_bbox = bbox(data)
-
-    nbbox = np.round(raw_bbox, decimals)
-    nbbox[0,2] = np.floor(raw_bbox * 10 ** decimals)[0,2] / 10 ** decimals
-    nbbox[1,2] = np.ceil(raw_bbox * 10 ** decimals)[1,2] / 10 ** decimals
-
-    return nbbox

idefix/vxl.py

@@ -11,7 +11,7 @@ General functions to transform point clouds to voxels compatible with numpy.
 import logging
 import numpy as np
 import humanize
-from .utils import bbox, fit_bbox
+from .utils import bbox
 
 log = logging.getLogger(__name__)
 
@@ -38,8 +38,8 @@ def _ui_step(step, spatial):
 def get_grid(spatial, step):
     '''Return grid bins.
 
-    Compute the grid bins of a point cloud or the corresponding bounding
+    Compute the grid bins of a point cloud or the corresponding bounding box
-    box according to given step (or steps for anisotropic grid).
+    according to given step (or steps for anisotropic grid).
 
     Parameters
     ----------
@@ -54,26 +54,8 @@ def get_grid(spatial, step):
     Returns
     -------
     grid : array of array (n,)
-        Grid of spatial given step. Return three arrays (not necessarily
+        Grid of spatial given step. Return three arrays (not necessarily of the
-        of the same size) defining the bins of axis `x`, `y` and `z`.
+        same size) defining the bins of axis `x`, `y` and `z`.
-
-    Notes
-    -----
-    The grid is built considering the half-open interval
-    $[min_of_the_axis, max_of_the_axis)$. If the positions of the points
-    are directly used as spatial parameter, the points at the upper
-    limits will be excluded from further processing.
-
-    You can define a more precise bounding box to take into account all
-    the points of your dataset (e.g. by adding a small distance on the
-    upper limits).
-
-    See Also
-    --------
-    bbox : Returns the raw bounding box of the point cloud (excluding
-        points on upper limit).
-    fit_bbox : Returns a bounding box on rounded coordinates (can
-        include all the points).
     '''
     spatial = np.array(spatial)
     bb = bbox(spatial)
@@ -83,7 +65,7 @@ def get_grid(spatial, step):
     for a_min, a_max, a_s in zip(bb[0], bb[1], step):
         # Beware of float underflow
         if a_s:
-            bins = np.trunc((a_max - a_min) / a_s).astype(int)
+            bins = np.trunc((a_max - a_min) / a_s).astype(int) + 1
             grid += [np.linspace(a_min, a_min + bins * a_s, bins + 1)]
         else:
             grid += [np.array((a_min, a_max + 1))]
@@ -365,7 +347,7 @@ def squash(voxel_grid, method='top', axis=-1):
         return voxel_grid.max(axis=axis)
     elif method == 'std':
         return voxel_grid.std(axis=axis)
-
+    
     raise NotImplementedError('Method \'{}\' does not exist.'.format(method))
 
 def plot(voxel_grid, vmin=None, vmax=None):

setup.py

@@ -15,7 +15,7 @@ with open('README.md', 'r') as fh:
 
 setuptools.setup(
     name='idefix',
-    version='0.2.0',
+    version='0.1.5',
     description='Utils and processing pipelines for LiDAR point clouds',
     author='Florent Guiotte',
     author_email='florent.guiotte@uhb.fr',
@@ -23,10 +23,7 @@ setuptools.setup(
     long_description=long_description,
     long_description_content_type='text/markdown',
     packages=['idefix', 'idefix.tools'],
-    entry_points = {'console_scripts':[
+    entry_points = {'console_scripts':['txt2npz = idefix.tools.txt_to_npz:main',]},
-        'txt2npz=idefix.tools.txt_to_npz:main',
-        'rasterize=idefix.tools.rasterize:main',
-    ]},
     classifiers=[
         'Programming Language :: Python :: 3',
         'License :: OSI Approved',

test/test_helpers.py

@@ -44,16 +44,3 @@ def test_dtm(ma_raster):
     dtm = helpers.dtm_dh_filter(ma_raster)
 
     assert dtm is not None, 'Did not return anything...'
-
-@pytest.mark.parametrize('params', [
-    {},
-    {'bin_structure': 'pixel'},
-    {'out_dir': True, 'crs': 'EPSG:26910'}])
-def test_rasterize(datadir, params):
-    # Workaround for out_dir with pytest
-    if 'out_dir' in params:
-        params['out_dir'] = datadir
-
-    raster = helpers.rasterize(datadir.join('test.npz'), **params)
-
-    assert raster is not None, 'Did not return anything...'

test/test_utils.py

@@ -27,11 +27,6 @@ def test_bbox(fix_data):
     res = np.array([fix_data.min(axis=0), fix_data.max(axis=0)])
     assert (utils.bbox(fix_data) == res).all()
 
-def test_fit_bbox(fix_data):
-    res = np.array([[0, 0, 0],
-                    [1, 1, 1]])
-    assert (utils.fit_bbox(fix_data) == res).all()
-
 def test_read(datadir):
     with open(datadir.join('first.txt')) as f:
         assert f.read() == 'hullo\n'

test/test_vxl.py

@@ -86,14 +86,7 @@ def test_get_grid(datadir, set_id, step, grid_id):
     assert res is not None, 'Test data empty, test function is broken!'
     assert len(res) == 3, 'Test data malformed, test function is broken!'
 
-    bbox = vxl.bbox(spatial)
+    test = vxl.get_grid(spatial, step)
-    try:
-        bbox[1,:] += step * 1.01
-    except:
-        # Workaround for the None on last axis
-        bbox[1,:] += step[0]
-
-    test = vxl.get_grid(bbox, step)
 
     assert test is not None, 'Function did not return anything :('
     assert len(test) == 3, 'Function doesn\'t give right number of axis'
@@ -185,7 +178,7 @@ def test__geo_to_np_coordinate():
     raster_truth[0, -1] = 25
     raster_truth[-1, 2] = 7
 
-    assert (raster_truth == vxl._geo_to_np_coordinate(raster)).all(), 'Missmatch between 2D raters'
+    assert (raster_truth == vxl._geo_to_np_coordinate(raster)).all(), 'Missmatch between 2D raters' 
 
     raster = np.zeros((5, 5, 3), dtype=np.uint8)
     raster[0, 0, 0] = 42
@@ -197,7 +190,7 @@ def test__geo_to_np_coordinate():
     raster_truth[0, -1, 1] = 25
     raster_truth[-1, 2, 2] = 7
 
-    assert (raster_truth == vxl._geo_to_np_coordinate(raster)).all(), 'Missmatch between 3D raters'
+    assert (raster_truth == vxl._geo_to_np_coordinate(raster)).all(), 'Missmatch between 3D raters' 
 
 @pytest.mark.parametrize('set_id, grid_id, axis, method', [
     (1, 1, 2, 'top'),

test/test_vxl/pc0_grid_s1-1-n.txt

@@ -1,3 +1,3 @@
 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0
 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0
-1.0 12.0
+1.0 11.0