idefix/utils.py

@@ -16,8 +16,10 @@ Everything should be highly tested there.
 
 import numpy as np
 import logging
+from scipy.interpolate import griddata
+from rasterio import fill
 
-log = logging.getLogger(__name__) 
+log = logging.getLogger(__name__)
 
 def first(a):
     """Returns the inverse of the parameter.
@@ -39,7 +41,7 @@ def first(a):
 
 def bbox(data):
     """Returns bounding box of data.
-    
+
     This function returns the lower and the upper points describing the
     bounding box of the points contained in data.
 
@@ -57,3 +59,47 @@ 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

setup.py

@@ -8,14 +8,34 @@
 #
 # TODO details
 
-from distutils.core import setup
+import setuptools
 
-setup(name='idefix',
+with open('README.md', 'r') as fh:
-      version='1.4',
+    long_description = fh.read()
-      description='Utils and processing pipelines for LiDAR point clouds',
+
-      author='Florent Guiotte',
+setuptools.setup(
-      author_email='florent.guiotte@uhb.fr',
+    name='idefix',
-      url='https://git.guiotte.fr/Florent/Idefix',
+    version='0.1.4',
-      packages=['idefix', 'idefix.tools'],
+    description='Utils and processing pipelines for LiDAR point clouds',
-      entry_points = {'console_scripts':['txt2npz = idefix.tools.txt_to_npz:main',]},
+    author='Florent Guiotte',
+    author_email='florent.guiotte@uhb.fr',
+    url='https://git.guiotte.fr/Florent/Idefix',
+    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',]},
+    classifiers=[
+        'Programming Language :: Python :: 3',
+        'License :: OSI Approved',
+        'Operating System :: OS Independent',
+    ],
+    python_requires='>=3.6',
+    install_requires=[
+        'numpy',
+        'sap',
+        'tqdm',
+        'matplotlib',
+        'pathlib',
+        'rasterio'
+    ],
       )

test/test_utils.py

@@ -11,6 +11,13 @@ 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),
@@ -30,3 +37,9 @@ 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)
+