Fix for Python < 3.9

idefix/helpers.py

@@ -14,10 +14,13 @@ 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
+from .vxl import get_grid, bin, squash, fit_bbox
+from .io import load_pc
 
 def interpolate(raster, method='linear'):
     """Interpolate masked raster.
@@ -99,7 +102,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
@@ -117,12 +120,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
 
@@ -130,8 +133,98 @@ 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

@@ -0,0 +1,101 @@
+#!/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,5 +55,42 @@ 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
+from .utils import bbox, fit_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 box
-    according to given step (or steps for anisotropic grid).
+    Compute the grid bins of a point cloud or the corresponding bounding
+    box according to given step (or steps for anisotropic grid).
 
     Parameters
     ----------
@@ -54,8 +54,26 @@ def get_grid(spatial, step):
     Returns
     -------
     grid : array of array (n,)
-        Grid of spatial given step. Return three arrays (not necessarily of the
-        same size) defining the bins of axis `x`, `y` and `z`.
+        Grid of spatial given step. Return three arrays (not necessarily
+        of the 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)
@@ -65,7 +83,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) + 1
+            bins = np.trunc((a_max - a_min) / a_s).astype(int)
             grid += [np.linspace(a_min, a_min + bins * a_s, bins + 1)]
         else:
             grid += [np.array((a_min, a_max + 1))]
@@ -347,7 +365,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.1.5',
+    version='0.2.0',
     description='Utils and processing pipelines for LiDAR point clouds',
     author='Florent Guiotte',
     author_email='florent.guiotte@uhb.fr',
@@ -23,7 +23,10 @@ setuptools.setup(
     long_description=long_description,
     long_description_content_type='text/markdown',
     packages=['idefix', 'idefix.tools'],
-    entry_points = {'console_scripts':['txt2npz = idefix.tools.txt_to_npz:main',]},
+    entry_points = {'console_scripts':[
+        '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,3 +44,16 @@ 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,6 +27,11 @@ 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,7 +86,14 @@ 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!'
 
-    test = vxl.get_grid(spatial, step)
+    bbox = vxl.bbox(spatial)
+    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'
@@ -178,7 +185,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
@@ -190,7 +197,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 11.0
+1.0 12.0