florent/ld2daps
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Global griddata and histogram voxel rasterisation

Browse Source
This commit is contained in:
Florent Guiotte 2018-07-02 18:19:45 +02:00
parent c27b6ac2fe
commit 16ca49cdb5
6 changed files with 1944 additions and 49 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

61
Notebooks/Custom rasters from LiDAR data.ipynb
View File

@ -69,8 +69,8 @@
" print('file: {}, point count: {}'.format(file, infile.header.get_count()))\n",
" xdata.append(infile.x)\n",
" ydata.append(infile.y)\n",
" #idata.append(infile.intensity)\n",
" idata.append(infile.num_returns)\n",
" idata.append(infile.intensity)\n",
" #idata.append(infile.num_returns)\n",
" \n",
" \n",
"data = np.array((np.concatenate(xdata), np.concatenate(ydata), np.concatenate(idata))).T\n",
@ -83,7 +83,7 @@
"metadata": {},
"outputs": [],
"source": [
"plt.hist(data[:,2], 4)\n",
"plt.hist(data[:,2], 1000)\n",
"plt.show()"
]
},
@ -101,7 +101,7 @@
"filtered_data[data[:,2] < tresholds[0]] = tresholds[0]\n",
"filtered_data[data[:,2] > tresholds[1]] = tresholds[1]\n",
"\n",
"plt.hist(filtered_data, 100)\n",
"plt.hist(filtered_data, 1000)\n",
"plt.show()"
]
},
@ -409,12 +409,12 @@
"sampled_data = np.array([[0, 0, 5], [9, 9, 5], [0, 9, 10], [20, 0, 10]])\n",
"\n",
"def rasterize(coords, values, resolution, method='cubic'):\n",
" #if isinstance(coords, np.array):\n",
" # xmin, xmax = coords[:,0].min(), coords[:,0].max()\n",
" # ymin, ymax = coords[:,1].min(), coords[:,1].max()\n",
" #else:\n",
" xmin, xmax = coords[0].min(), coords[0].max()\n",
" ymin, ymax = coords[1].min(), coords[1].max()\n",
" if isinstance(coords, np.ndarray):\n",
" xmin, xmax = coords[:,0].min(), coords[:,0].max()\n",
" ymin, ymax = coords[:,1].min(), coords[:,1].max()\n",
" else:\n",
" xmin, xmax = coords[0].min(), coords[0].max()\n",
" ymin, ymax = coords[1].min(), coords[1].max()\n",
"\n",
" gridx, gridy = np.mgrid[xmin:xmax:resolution, ymin:ymax:resolution]\n",
" \n",
@ -471,7 +471,7 @@
"outputs": [],
"source": [
"#raster = rasterize(data[:treshold,:2], should_have_keep_idata[:treshold], 0.5, 'nearest')\n",
"raster = rasterize((infile.x, infile.y), infile.num_returns, 1, 'nearest')\n",
"raster = rasterize((data[:,0], data[:,1]), data[:,2], 0.5, 'nearest')\n",
"\n",
"plt.figure(figsize=figsize)\n",
"plt.imshow(raster.T, origin='lower')\n",
@ -484,44 +484,7 @@
"metadata": {},
"outputs": [],
"source": [
"plt.imsave('../Res/my_raster_cubic.png', raster.T, origin='lower')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"idisp_complete = griddata(data[:,:2], data[:,2], test)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"X, Y = np.meshgrid(*([np.linspace(-1,1,200)] * 2))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"test = np.mgrid[0:bins[1], 0:bins[0]]\n",
"test"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"data[:,:2]"
"plt.imsave('../Res/raster.png', raster.T, origin='lower')"
]
},
{

461
Notebooks/Global Griddata Rasterisation Strategies.ipynb Normal file
View File

@ -0,0 +1,461 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import sys\n",
"from pathlib import Path\n",
"import numpy as np\n",
"from scipy.interpolate import griddata\n",
"import matplotlib.pyplot as plt\n",
"import laspy\n",
"\n",
"triskele_path = Path('../triskele/python/')\n",
"sys.path.append(str(triskele_path.resolve()))\n",
"import triskele\n",
"\n",
"figsize = np.array((16, 3)) * 1.5"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Setup"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load LiDAR data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def find_las(path):\n",
" path = Path(path)\n",
" las = list()\n",
" \n",
" if path.is_dir():\n",
" for child in path.iterdir():\n",
" las.extend(find_las(child))\n",
" \n",
" if path.is_file() and path.suffix == '.las':\n",
" las.append(path)\n",
" \n",
" return las"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def auto_filter(data, treshold=.5):\n",
" tresholds = np.percentile(data, [treshold, 100 - treshold])\n",
"\n",
" data[data < tresholds[0]] = tresholds[0]\n",
" data[data > tresholds[1]] = tresholds[1]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def bulk_load(path, name=None):\n",
" data = {'file': path, 'name': name}\n",
"\n",
" attributes = ['x', 'y', 'z', 'intensity', 'num_returns', 'scan_angle_rank', 'pt_src_id', 'gps_time']\n",
" for a in attributes:\n",
" data[a] = list()\n",
" \n",
" print('Load data...')\n",
" for f in find_las(path):\n",
" print('{}: '.format(f), end='')\n",
" infile = laspy.file.File(f)\n",
" for i, a in enumerate(attributes):\n",
" print('\\r {}: [{:3d}%]'.format(f, int(i/len(attributes) * 100)), end='')\n",
" data[a].extend(getattr(infile, a))\n",
" infile.close()\n",
" print('\\r {}: [Done]'.format(f))\n",
"\n",
" \n",
" print('Create matrices...', end='')\n",
" for i, a in enumerate(attributes):\n",
" print('\\rCreate matrices: [{:3d}%]'.format(int(i/len(attributes) * 100)), end='')\n",
" data[a] = np.array(data[a])\n",
" print('\\rCreate matrices: [Done]')\n",
" \n",
" print('Filter data...', end='')\n",
" for i, a in enumerate(attributes):\n",
" print('\\rFilter data: [{:3d}%]'.format(int(i/len(attributes) * 100)), end='')\n",
" if a == 'x' or a == 'y':\n",
" continue\n",
" auto_filter(data[a])\n",
" print('\\rFilter data: [Done]')\n",
"\n",
" return data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"C1 = bulk_load('../Data/lidar/C1/', 'C1')\n",
"C2 = bulk_load('../Data/lidar/C2/', 'C2')\n",
"C3 = bulk_load('../Data/lidar/C3/', 'C3')\n",
"C123 = bulk_load('../Data/lidar/', 'C123') # Todo: Copy C1 C2 C3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define raster function"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def rasterize(coords, values, resolution=1, method='nearest'):\n",
" if isinstance(coords, np.ndarray):\n",
" xmin, xmax = coords[:,0].min(), coords[:,0].max()\n",
" ymin, ymax = coords[:,1].min(), coords[:,1].max()\n",
" else:\n",
" xmin, xmax = coords[0].min(), coords[0].max()\n",
" ymin, ymax = coords[1].min(), coords[1].max()\n",
"\n",
" gridx, gridy = np.mgrid[xmin:xmax:resolution, ymin:ymax:resolution]\n",
" \n",
" return griddata(coords, values, (gridx, gridy), method=method).T # linear, nearest, cubic"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def raster_show_and_save(coords, values, resolution=1, method='nearest', name=None, save_dir='../Res'):\n",
" \"\"\"Display and save rasters. If name is provided, will try to cache result to reload it later.\"\"\"\n",
" \n",
" if name is not None:\n",
" base_name = '{}/{}'.format(save_dir, name)\n",
" tiff_file = Path(base_name + '.tiff')\n",
" png_file = Path(base_name + '.png')\n",
" \n",
" if name is not None and tiff_file.exists() :\n",
" print ('WARNING: Loading cached result {}'.format(tiff_file))\n",
" raster = triskele.read(tiff_file)\n",
" else: \n",
" raster = rasterize(coords, values, resolution, method)\n",
" \n",
" plt.figure(figsize=figsize)\n",
" plt.imshow(raster, origin='lower')\n",
" plt.colorbar()\n",
" plt.title(name)\n",
" plt.show()\n",
" \n",
" if name is not None:\n",
" plt.imsave(png_file, raster, origin='lower')\n",
" triskele.write(tiff_file, raster)\n",
" \n",
" return raster"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def raster_assistant(data, field='z', resolution='1', method='nearest', weight_field=None):\n",
" if weight_field is None:\n",
" name = '{}_{}_{}_{}'.format(data['name'], field, str(resolution).replace('.', '_'), method)\n",
" raster_show_and_save((data['x'], data['y']), data[field], resolution, method, name)\n",
" else:\n",
" name = '{}_{}_{}_{}_w{}'.format(data['name'], field, str(resolution).replace('.', '_'), method, weight_field)\n",
" raster_show_and_save((data['x'], data['y']), data[field] * data[weight_field], resolution, method, name)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def raster_assistant_mp(data, field='z', resolution='1', method='nearest', weight_field=None):\n",
" if data == 'C1':\n",
" raster_assistant(C1, field, resolution, method)\n",
" if data == 'C2':\n",
" raster_assistant(C2, field, resolution, method)\n",
" if data == 'C3':\n",
" raster_assistant(C3, field, resolution, method)\n",
" if data == 'C123':\n",
" raster_assistant(C123, field, resolution, method)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rl = raster_assistant(C1, resolution=10)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Interpolation Strategies for Global Griddata Rasterisation\n",
"\n",
"## Fine textured (visible) rasters\n",
"\n",
"### Nearest"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"null = raster_assistant(C1, 'intensity', 0.5, 'nearest')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"No interpolation, accurate values. On the downside, there is a lot of aliasing. Maybe that's because take all the points without selection on the z (altitude) axis."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Linear"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"null = raster_assistant(C1, 'intensity', 0.5, 'linear')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Cubic"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The results are pretty bad, same remark than before."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"null = raster_assistant(C1, 'intensity', 0.5, 'cubic')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The cubic create new values. Still, seems pretty useless with the aliasing. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Flat area (DSM) rasters\n",
"\n",
"### Nearest"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"null = raster_assistant(C1, 'z', 0.5, 'nearest')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"No interpolation, accurate values. On the downside, there is a lot of aliasing. Maybe that's because take all the points without selection on the z (altitude) axis."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Linear"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"null = raster_assistant(C1, 'z', 0.5, 'linear')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Cubic"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The results are pretty bad, same remark than before."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"null = raster_assistant(C1, 'z', 0.5, 'cubic')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Weight with altitude\n",
"\n",
"## Visible"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"null = raster_assistant(C1, 'intensity', 0.5, 'nearest', 'z')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"null = raster_assistant(C1, 'intensity', 0.5, 'linear', 'z')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## DSM"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"null = raster_assistant(C1, 'z', 0.5, 'nearest', 'z')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"null = raster_assistant(C1, 'z', 0.5, 'linear', 'z')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Other"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from multiprocessing import Pool, Process, Queue\n",
"import multiprocessing as mp\n",
"#mp.set_start_method('spawn')\n",
"\n",
"pool = Pool(processes=20)\n",
"\n",
"job_args = list()\n",
"for field in ('z', 'intensity', 'num_returns', 'gps_time'):\n",
" for data in ('C1', 'C2', 'C3', 'C123'):\n",
" job_args.append([data, field, 1, 'nearest'])\n",
" \n",
"for i in pool.starmap(raster_assistant_mp, job_args):\n",
" pass "
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

484
Notebooks/Histogram Voxels Rasterisation.ipynb Normal file
View File

@ -0,0 +1,484 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import sys\n",
"from pathlib import Path\n",
"import numpy as np\n",
"from scipy.interpolate import griddata\n",
"import matplotlib.pyplot as plt\n",
"import laspy\n",
"\n",
"triskele_path = Path('../triskele/python/')\n",
"sys.path.append(str(triskele_path.resolve()))\n",
"import triskele\n",
"\n",
"figsize = np.array((16, 3)) * 1.5"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Setup"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load LiDAR data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def find_las(path):\n",
" path = Path(path)\n",
" las = list()\n",
" \n",
" if path.is_dir():\n",
" for child in path.iterdir():\n",
" las.extend(find_las(child))\n",
" \n",
" if path.is_file() and path.suffix == '.las':\n",
" las.append(path)\n",
" \n",
" return las"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def auto_filter(data, treshold=.5):\n",
" tresholds = np.percentile(data, [treshold, 100 - treshold])\n",
"\n",
" data[data < tresholds[0]] = tresholds[0]\n",
" data[data > tresholds[1]] = tresholds[1]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def bulk_load(path, name=None):\n",
" data = {'file': path, 'name': name}\n",
"\n",
" attributes = ['x', 'y', 'z', 'intensity', 'num_returns']#, 'scan_angle_rank', 'pt_src_id', 'gps_time']\n",
" for a in attributes:\n",
" data[a] = list()\n",
" \n",
" print('Load data...')\n",
" for f in find_las(path):\n",
" print('{}: '.format(f), end='')\n",
" infile = laspy.file.File(f)\n",
" for i, a in enumerate(attributes):\n",
" print('\\r {}: [{:3d}%]'.format(f, int(i/len(attributes) * 100)), end='')\n",
" data[a].extend(getattr(infile, a))\n",
" infile.close()\n",
" print('\\r {}: [Done]'.format(f))\n",
"\n",
" \n",
" print('Create matrices...', end='')\n",
" for i, a in enumerate(attributes):\n",
" print('\\rCreate matrices: [{:3d}%]'.format(int(i/len(attributes) * 100)), end='')\n",
" data[a] = np.array(data[a])\n",
" print('\\rCreate matrices: [Done]')\n",
" \n",
" print('Filter data...', end='')\n",
" for i, a in enumerate(attributes):\n",
" print('\\rFilter data: [{:3d}%]'.format(int(i/len(attributes) * 100)), end='')\n",
" if a == 'x' or a == 'y':\n",
" continue\n",
" auto_filter(data[a])\n",
" print('\\rFilter data: [Done]')\n",
"\n",
" return data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"C1 = bulk_load('../Data/lidar/C1', 'C1')\n",
"#C2 = bulk_load('../Data/lidar/C2/', 'C2')\n",
"#C3 = bulk_load('../Data/lidar/C3/', 'C3')\n",
"#C123 = bulk_load('../Data/lidar/', 'C123') # Todo: Copy C1 C2 C3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define raster method"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"resolution = 0.5\n",
"\n",
"spatial_data = np.array((C1['x'], C1['y'], C1['z'])).T\n",
"\n",
"bins = np.rint((spatial_data.max(axis=0) - spatial_data.min(axis=0)) / resolution).astype(int)\n",
"display(bins)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"voxels_histo, edges = np.histogramdd(spatial_data, bins=bins)\n",
"voxels_histo.shape"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"voxels_histo_i, edges = np.histogramdd(spatial_data, bins=bins, weights=C1['intensity'])\n",
"voxels_histo_i.shape"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"layers = voxels_histo.shape[2]\n",
"\n",
"plt.figure(figsize=figsize * np.array((1, layers)))\n",
"\n",
"for i in range(layers):\n",
" plt.subplot(layers, 1, i+1)\n",
" plt.title('Layer {}'.format(i))\n",
" plt.imshow(voxels_histo_i.T[i] / voxels_histo.T[i], origin='lower')\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.imsave('../Res/voxel_layer', voxels_histo.T[0], origin='lower')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create voxels and flatten raster"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"layers = voxels_histo.shape[2]\n",
"\n",
"raster = np.zeros_like(voxels_histo.T[0])\n",
"raster[:] = np.nan\n",
"\n",
"voxels = (voxels_histo_i / voxels_histo).T\n",
"\n",
"for i in reversed(range(layers)):\n",
" nan_filter = np.isnan(raster)\n",
" raster[nan_filter] = voxels[i][nan_filter]\n",
"\n",
"raster"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=figsize)\n",
"plt.imshow(raster, origin='lower')\n",
"plt.show()\n",
"plt.imsave('../Res/voxel_raster_i_0_5.png', raster, origin='lower')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Interpolate NoData"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"nan_filter = np.isnan(raster)\n",
"val_filter = np.logical_not(nan_filter)\n",
"coords = np.argwhere(val_filter)\n",
"values = raster[val_filter]\n",
"grid = np.argwhere(nan_filter)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"interpol_values_linear = griddata(coords, values, grid, method='linear')\n",
"interpol_values_cubic = griddata(coords, values, grid, method='cubic')\n",
"interpol_values_nearest = griddata(coords, values, grid, method='nearest')\n",
"\n",
"raster_interpol_linear = raster.copy()\n",
"raster_interpol_cubic = raster.copy()\n",
"raster_interpol_nearest = raster.copy()\n",
"\n",
"raster_interpol_linear[nan_filter] = interpol_values_linear\n",
"raster_interpol_cubic[nan_filter] = interpol_values_cubic\n",
"raster_interpol_nearest[nan_filter] = interpol_values_nearest"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=figsize)\n",
"plt.imshow(raster_interpol_linear, origin='lower')\n",
"plt.show()\n",
"plt.imsave('../Res/voxel_raster_interpol_tmp_i_0_5.png', raster_interpol_linear, origin='lower')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Why do we have high values there?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raster_interpol_nearest.max(), raster_interpol_linear.max()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raster_interpol_nearest.min(), raster_interpol_linear.min()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raster_interpol_linear[raster_interpol_linear > 1e6] = 0"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=figsize)\n",
"plt.imshow((raster_interpol_linear - raster_interpol_nearest))\n",
"plt.colorbar()\n",
"plt.show()\n",
"plt.imsave('../Res/diff_interpol.png', raster_interpol_linear - raster_interpol_nearest, origin='lower')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(raster_interpol_linear - raster_interpol_nearest).max()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raster_interpol_cubic[np.isnan(raster_interpol_cubic)] = np.nanargmin(raster_interpol_cubic)\n",
"raster_interpol_linear[np.isnan(raster_interpol_linear)] = np.nanargmin(raster_interpol_linear)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#plt.hist(raster_interpol_cubic.reshape(-1), bins=100, label='cubic',log=True)\n",
"plt.hist(raster_interpol_linear[raster_interpol_linear < 1e6].reshape(-1), 100, label='linear',log=False)\n",
"plt.hist(raster_interpol_nearest.reshape(-1), 100, label='nearest',log=False, alpha=.8)\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"np.isnan(raster_interpol_nearest).sum()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Function "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def rasterize(spatial, values, resolution=1., method='nearest'):\n",
" '''\n",
" Rasterize spatialised data.\n",
"\n",
" :param spatial: Coordinates of the points (x, y, z)\n",
" :type spatial: ndarray\n",
" :param values: Values of the points\n",
" :type values: ndarray\n",
" :param resolution: Resolution of the raster in meters\n",
" :type resolution: float\n",
" :param method: Method of interpolation for NoData values {'linear', 'nearest', 'cubic'}\n",
" :type method: string\n",
" :return: 2D raster of the values along z axis\n",
" :rtype: ndarray\n",
"\n",
" .. warning:: WIP\n",
" '''\n",
" \n",
" resolution = 0.5\n",
"\n",
"spatial_data = np.array((C1['x'], C1['y'], C1['z'])).T\n",
"\n",
"bins = np.rint((spatial_data.max(axis=0) - spatial_data.min(axis=0)) / resolution).astype(int)\n",
"display(bins)\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rasterize griddata"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"infile = laspy.file.File('../Data/lidar/C1/272056_3289689.las')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"infile.reader.get_dimension('X')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# TMP"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dfc_c1 = triskele.read('../Data/phase1_rasters/Intensity_C1/UH17_GI1F051_TR.tif')\n",
"auto_filter(dfc_c1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.imshow(dfc_c1)\n",
"plt.show()\n",
"plt.imsave('../Res/dfc_c1.png', dfc_c1)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

191
Notebooks/Triskel check.ipynb Normal file
View File

@ -0,0 +1,191 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import sys\n",
"from pathlib import Path\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import laspy\n",
"\n",
"triskele_path = Path('../triskele/python/')\n",
"sys.path.append(str(triskele_path.resolve()))\n",
"import triskele\n",
"\n",
"figsize = (16, 3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Read and Write"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"R = np.random.random((1000,4000))\n",
"\n",
"plt.figure(figsize=figsize)\n",
"plt.imshow(R)\n",
"plt.colorbar()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"R.dtype"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## F64 "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"triskele.write('../Res/rw_test.tiff', R)\n",
"R2 = triskele.read('../Res/rw_test.tiff')\n",
"R2.dtype"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=figsize)\n",
"plt.imshow(R2)\n",
"plt.colorbar()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## F32"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=figsize)\n",
"plt.imshow(R.astype(np.float32))\n",
"plt.colorbar()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"triskele.write('../Res/rw_test.tiff', R.astype(np.float32))\n",
"R2 = triskele.read('../Res/rw_test.tiff')\n",
"R2.dtype"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=figsize)\n",
"plt.imshow(R2)\n",
"plt.colorbar()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Int32"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"R = np.random.random_integers(0, 10000, (1000,4000)).astype(np.int32)\n",
"\n",
"plt.figure(figsize=figsize)\n",
"plt.imshow(R)\n",
"plt.colorbar()\n",
"plt.show()\n",
"\n",
"R.dtype"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"triskele.write('../Res/rw_test.tiff', R)\n",
"R2 = triskele.read('../Res/rw_test.tiff')\n",
"R2.dtype"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=figsize)\n",
"plt.imshow(R2)\n",
"plt.colorbar()\n",
"plt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

224
Notebooks/Triskele Python Demo.ipynb Normal file
View File

@ -0,0 +1,224 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Init"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Import Triskele\n",
"\n",
"I did not install triskele-python on my system folders, so this is a trick to load the module outside the triskele project."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import sys\n",
"sys.path.append('../triskele/python/') # Change with your triskele install path\n",
"import triskele"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Simple test\n",
"\n",
"Just to test the python import and the triskele compilation. Create an array with random values."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raster = np.random.randint(0, 255, (1000, 1000))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create the instance of triskele object with the array to filter."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ap_random = triskele.Triskele(raster)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Filter the raster with this method :"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"filtered_raster = ap_random.filter()\n",
"print(filtered_raster)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If there is no error so far everything is okay.\n",
"\n",
"Theoretically the filtered output should be a copy of the input. (it's not, there is something going on along the edges, we should tell François !) \n",
"\n",
"## Attribute filter\n",
"\n",
"With `filter` method you can ask for specific attributes filtering. The results are stacked."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"filtered_rasters = ap_random.filter(area=[10, 100, 1e3])\n",
"print(filtered_rasters.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I don't know if you can see the signature with yout IDE : `filter(tree='max-tree', area=None, standard_deviation=None, moment_of_inertia=None)`\n",
"\n",
"The parameters are the same than apGenerator from François, I should write the doc..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Advanced test\n",
"\n",
"With a real example.\n",
"\n",
"Load a TIFF file with 4 bands :"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raster = triskele.read('../triskele/data/10m.tif')\n",
"print(raster.shape)\n",
"\n",
"plt.imshow(raster)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create the instance with some tweaks."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ap_10m = triskele.Triskele(raster, verbose=False) # Verbose False to hide apGenerator gibbering :p"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Filter as you wish :"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tree = 'tos-tree' # {'max-tree', 'min-tree', ...}\n",
"area = [100, 1e3, 1e4]\n",
"standard_deviation = None\n",
"moment_of_inertia = None\n",
"\n",
"filtered_rasters = ap_10m.filter(tree=tree, area=area, standard_deviation=standard_deviation, moment_of_inertia=moment_of_inertia)\n",
"print(filtered_rasters.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Display results"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for i in range(filtered_rasters.shape[2]):\n",
" plt.imshow(filtered_rasters[:,:,i])\n",
" plt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

572
Notebooks/notebook.tex Normal file
View File

@ -0,0 +1,572 @@
% Default to the notebook output style
% Inherit from the specified cell style.
\documentclass[11pt]{article}
\usepackage[T1]{fontenc}
% Nicer default font (+ math font) than Computer Modern for most use cases
\usepackage{mathpazo}
% Basic figure setup, for now with no caption control since it's done
% automatically by Pandoc (which extracts ![](path) syntax from Markdown).
\usepackage{graphicx}
% We will generate all images so they have a width \maxwidth. This means
% that they will get their normal width if they fit onto the page, but
% are scaled down if they would overflow the margins.
\makeatletter
\def\maxwidth{\ifdim\Gin@nat@width>\linewidth\linewidth
\else\Gin@nat@width\fi}
\makeatother
\let\Oldincludegraphics\includegraphics
% Set max figure width to be 80% of text width, for now hardcoded.
\renewcommand{\includegraphics}[1]{\Oldincludegraphics[width=.8\maxwidth]{#1}}
% Ensure that by default, figures have no caption (until we provide a
% proper Figure object with a Caption API and a way to capture that
% in the conversion process - todo).
\usepackage{caption}
\DeclareCaptionLabelFormat{nolabel}{}
\captionsetup{labelformat=nolabel}
\usepackage{adjustbox} % Used to constrain images to a maximum size
\usepackage{xcolor} % Allow colors to be defined
\usepackage{enumerate} % Needed for markdown enumerations to work
\usepackage{geometry} % Used to adjust the document margins
\usepackage{amsmath} % Equations
\usepackage{amssymb} % Equations
\usepackage{textcomp} % defines textquotesingle
% Hack from http://tex.stackexchange.com/a/47451/13684:
\AtBeginDocument{%
\def\PYZsq{\textquotesingle}% Upright quotes in Pygmentized code
}
\usepackage{upquote} % Upright quotes for verbatim code
\usepackage{eurosym} % defines \euro
\usepackage[mathletters]{ucs} % Extended unicode (utf-8) support
\usepackage[utf8x]{inputenc} % Allow utf-8 characters in the tex document
\usepackage{fancyvrb} % verbatim replacement that allows latex
\usepackage{grffile} % extends the file name processing of package graphics
% to support a larger range
% The hyperref package gives us a pdf with properly built
% internal navigation ('pdf bookmarks' for the table of contents,
% internal cross-reference links, web links for URLs, etc.)
\usepackage{hyperref}
\usepackage{longtable} % longtable support required by pandoc >1.10
\usepackage{booktabs} % table support for pandoc > 1.12.2
\usepackage[inline]{enumitem} % IRkernel/repr support (it uses the enumerate* environment)
\usepackage[normalem]{ulem} % ulem is needed to support strikethroughs (\sout)
% normalem makes italics be italics, not underlines
% Colors for the hyperref package
\definecolor{urlcolor}{rgb}{0,.145,.698}
\definecolor{linkcolor}{rgb}{.71,0.21,0.01}
\definecolor{citecolor}{rgb}{.12,.54,.11}
% ANSI colors
\definecolor{ansi-black}{HTML}{3E424D}
\definecolor{ansi-black-intense}{HTML}{282C36}
\definecolor{ansi-red}{HTML}{E75C58}
\definecolor{ansi-red-intense}{HTML}{B22B31}
\definecolor{ansi-green}{HTML}{00A250}
\definecolor{ansi-green-intense}{HTML}{007427}
\definecolor{ansi-yellow}{HTML}{DDB62B}
\definecolor{ansi-yellow-intense}{HTML}{B27D12}
\definecolor{ansi-blue}{HTML}{208FFB}
\definecolor{ansi-blue-intense}{HTML}{0065CA}
\definecolor{ansi-magenta}{HTML}{D160C4}
\definecolor{ansi-magenta-intense}{HTML}{A03196}
\definecolor{ansi-cyan}{HTML}{60C6C8}
\definecolor{ansi-cyan-intense}{HTML}{258F8F}
\definecolor{ansi-white}{HTML}{C5C1B4}
\definecolor{ansi-white-intense}{HTML}{A1A6B2}
% commands and environments needed by pandoc snippets
% extracted from the output of `pandoc -s`
\providecommand{\tightlist}{%
\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}
\DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\}}
% Add ',fontsize=\small' for more characters per line
\newenvironment{Shaded}{}{}
\newcommand{\KeywordTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{\textbf{{#1}}}}
\newcommand{\DataTypeTok}[1]{\textcolor[rgb]{0.56,0.13,0.00}{{#1}}}
\newcommand{\DecValTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{{#1}}}
\newcommand{\BaseNTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{{#1}}}
\newcommand{\FloatTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{{#1}}}
\newcommand{\CharTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{{#1}}}
\newcommand{\StringTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{{#1}}}
\newcommand{\CommentTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textit{{#1}}}}
\newcommand{\OtherTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{{#1}}}
\newcommand{\AlertTok}[1]{\textcolor[rgb]{1.00,0.00,0.00}{\textbf{{#1}}}}
\newcommand{\FunctionTok}[1]{\textcolor[rgb]{0.02,0.16,0.49}{{#1}}}
\newcommand{\RegionMarkerTok}[1]{{#1}}
\newcommand{\ErrorTok}[1]{\textcolor[rgb]{1.00,0.00,0.00}{\textbf{{#1}}}}
\newcommand{\NormalTok}[1]{{#1}}
% Additional commands for more recent versions of Pandoc
\newcommand{\ConstantTok}[1]{\textcolor[rgb]{0.53,0.00,0.00}{{#1}}}
\newcommand{\SpecialCharTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{{#1}}}
\newcommand{\VerbatimStringTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{{#1}}}
\newcommand{\SpecialStringTok}[1]{\textcolor[rgb]{0.73,0.40,0.53}{{#1}}}
\newcommand{\ImportTok}[1]{{#1}}
\newcommand{\DocumentationTok}[1]{\textcolor[rgb]{0.73,0.13,0.13}{\textit{{#1}}}}
\newcommand{\AnnotationTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{{#1}}}}}
\newcommand{\CommentVarTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{{#1}}}}}
\newcommand{\VariableTok}[1]{\textcolor[rgb]{0.10,0.09,0.49}{{#1}}}
\newcommand{\ControlFlowTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{\textbf{{#1}}}}
\newcommand{\OperatorTok}[1]{\textcolor[rgb]{0.40,0.40,0.40}{{#1}}}
\newcommand{\BuiltInTok}[1]{{#1}}
\newcommand{\ExtensionTok}[1]{{#1}}
\newcommand{\PreprocessorTok}[1]{\textcolor[rgb]{0.74,0.48,0.00}{{#1}}}
\newcommand{\AttributeTok}[1]{\textcolor[rgb]{0.49,0.56,0.16}{{#1}}}
\newcommand{\InformationTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{{#1}}}}}
\newcommand{\WarningTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{{#1}}}}}
% Define a nice break command that doesn't care if a line doesn't already
% exist.
\def\br{\hspace*{\fill} \\* }
% Math Jax compatability definitions
\def\gt{>}
\def\lt{<}
% Document parameters
\title{Triskele Python Demo}
% Pygments definitions
\makeatletter
\def\PY@reset{\let\PY@it=\relax \let\PY@bf=\relax%
\let\PY@ul=\relax \let\PY@tc=\relax%
\let\PY@bc=\relax \let\PY@ff=\relax}
\def\PY@tok#1{\csname PY@tok@#1\endcsname}
\def\PY@toks#1+{\ifx\relax#1\empty\else%
\PY@tok{#1}\expandafter\PY@toks\fi}
\def\PY@do#1{\PY@bc{\PY@tc{\PY@ul{%
\PY@it{\PY@bf{\PY@ff{#1}}}}}}}
\def\PY#1#2{\PY@reset\PY@toks#1+\relax+\PY@do{#2}}
\expandafter\def\csname PY@tok@w\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.73,0.73}{##1}}}
\expandafter\def\csname PY@tok@c\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}}
\expandafter\def\csname PY@tok@cp\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.74,0.48,0.00}{##1}}}
\expandafter\def\csname PY@tok@k\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
\expandafter\def\csname PY@tok@kp\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
\expandafter\def\csname PY@tok@kt\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.69,0.00,0.25}{##1}}}
\expandafter\def\csname PY@tok@o\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
\expandafter\def\csname PY@tok@ow\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.67,0.13,1.00}{##1}}}
\expandafter\def\csname PY@tok@nb\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
\expandafter\def\csname PY@tok@nf\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,1.00}{##1}}}
\expandafter\def\csname PY@tok@nc\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,1.00}{##1}}}
\expandafter\def\csname PY@tok@nn\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,1.00}{##1}}}
\expandafter\def\csname PY@tok@ne\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.82,0.25,0.23}{##1}}}
\expandafter\def\csname PY@tok@nv\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}}
\expandafter\def\csname PY@tok@no\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.53,0.00,0.00}{##1}}}
\expandafter\def\csname PY@tok@nl\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.63,0.63,0.00}{##1}}}
\expandafter\def\csname PY@tok@ni\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.60,0.60,0.60}{##1}}}
\expandafter\def\csname PY@tok@na\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.49,0.56,0.16}{##1}}}
\expandafter\def\csname PY@tok@nt\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
\expandafter\def\csname PY@tok@nd\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.67,0.13,1.00}{##1}}}
\expandafter\def\csname PY@tok@s\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
\expandafter\def\csname PY@tok@sd\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
\expandafter\def\csname PY@tok@si\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.73,0.40,0.53}{##1}}}
\expandafter\def\csname PY@tok@se\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.73,0.40,0.13}{##1}}}
\expandafter\def\csname PY@tok@sr\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.40,0.53}{##1}}}
\expandafter\def\csname PY@tok@ss\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}}
\expandafter\def\csname PY@tok@sx\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
\expandafter\def\csname PY@tok@m\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
\expandafter\def\csname PY@tok@gh\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,0.50}{##1}}}
\expandafter\def\csname PY@tok@gu\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.50,0.00,0.50}{##1}}}
\expandafter\def\csname PY@tok@gd\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.63,0.00,0.00}{##1}}}
\expandafter\def\csname PY@tok@gi\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.63,0.00}{##1}}}
\expandafter\def\csname PY@tok@gr\endcsname{\def\PY@tc##1{\textcolor[rgb]{1.00,0.00,0.00}{##1}}}
\expandafter\def\csname PY@tok@ge\endcsname{\let\PY@it=\textit}
\expandafter\def\csname PY@tok@gs\endcsname{\let\PY@bf=\textbf}
\expandafter\def\csname PY@tok@gp\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,0.50}{##1}}}
\expandafter\def\csname PY@tok@go\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.53,0.53,0.53}{##1}}}
\expandafter\def\csname PY@tok@gt\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.27,0.87}{##1}}}
\expandafter\def\csname PY@tok@err\endcsname{\def\PY@bc##1{\setlength{\fboxsep}{0pt}\fcolorbox[rgb]{1.00,0.00,0.00}{1,1,1}{\strut ##1}}}
\expandafter\def\csname PY@tok@kc\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
\expandafter\def\csname PY@tok@kd\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
\expandafter\def\csname PY@tok@kn\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
\expandafter\def\csname PY@tok@kr\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
\expandafter\def\csname PY@tok@bp\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
\expandafter\def\csname PY@tok@fm\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,1.00}{##1}}}
\expandafter\def\csname PY@tok@vc\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}}
\expandafter\def\csname PY@tok@vg\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}}
\expandafter\def\csname PY@tok@vi\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}}
\expandafter\def\csname PY@tok@vm\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}}
\expandafter\def\csname PY@tok@sa\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
\expandafter\def\csname PY@tok@sb\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
\expandafter\def\csname PY@tok@sc\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
\expandafter\def\csname PY@tok@dl\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
\expandafter\def\csname PY@tok@s2\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
\expandafter\def\csname PY@tok@sh\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
\expandafter\def\csname PY@tok@s1\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
\expandafter\def\csname PY@tok@mb\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
\expandafter\def\csname PY@tok@mf\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
\expandafter\def\csname PY@tok@mh\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
\expandafter\def\csname PY@tok@mi\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
\expandafter\def\csname PY@tok@il\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
\expandafter\def\csname PY@tok@mo\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
\expandafter\def\csname PY@tok@ch\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}}
\expandafter\def\csname PY@tok@cm\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}}
\expandafter\def\csname PY@tok@cpf\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}}
\expandafter\def\csname PY@tok@c1\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}}
\expandafter\def\csname PY@tok@cs\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}}
\def\PYZbs{\char`\\}
\def\PYZus{\char`\_}
\def\PYZob{\char`\{}
\def\PYZcb{\char`\}}
\def\PYZca{\char`\^}
\def\PYZam{\char`\&}
\def\PYZlt{\char`\<}
\def\PYZgt{\char`\>}
\def\PYZsh{\char`\#}
\def\PYZpc{\char`\%}
\def\PYZdl{\char`\$}
\def\PYZhy{\char`\-}
\def\PYZsq{\char`\'}
\def\PYZdq{\char`\"}
\def\PYZti{\char`\~}
% for compatibility with earlier versions
\def\PYZat{@}
\def\PYZlb{[}
\def\PYZrb{]}
\makeatother
% Exact colors from NB
\definecolor{incolor}{rgb}{0.0, 0.0, 0.5}
\definecolor{outcolor}{rgb}{0.545, 0.0, 0.0}
% Prevent overflowing lines due to hard-to-break entities
\sloppy
% Setup hyperref package
\hypersetup{
breaklinks=true, % so long urls are correctly broken across lines
colorlinks=true,
urlcolor=urlcolor,
linkcolor=linkcolor,
citecolor=citecolor,
}
% Slightly bigger margins than the latex defaults
\geometry{verbose,tmargin=1in,bmargin=1in,lmargin=1in,rmargin=1in}
\begin{document}
\maketitle
\hypertarget{init}{%
\section{Init}\label{init}}
\begin{Verbatim}[commandchars=\\\{\}]
{\color{incolor}In [{\color{incolor}1}]:} \PY{k+kn}{import} \PY{n+nn}{numpy} \PY{k}{as} \PY{n+nn}{np}
\PY{k+kn}{import} \PY{n+nn}{matplotlib}\PY{n+nn}{.}\PY{n+nn}{pyplot} \PY{k}{as} \PY{n+nn}{plt}
\end{Verbatim}
\hypertarget{import-triskele}{%
\subsection{Import Triskele}\label{import-triskele}}
I did not install triskele-python on my system folders, so this is a
trick to load the module outside the triskele project.
\begin{Verbatim}[commandchars=\\\{\}]
{\color{incolor}In [{\color{incolor}2}]:} \PY{k+kn}{import} \PY{n+nn}{sys}
\PY{n}{sys}\PY{o}{.}\PY{n}{path}\PY{o}{.}\PY{n}{append}\PY{p}{(}\PY{l+s+s1}{\PYZsq{}}\PY{l+s+s1}{../triskele/python/}\PY{l+s+s1}{\PYZsq{}}\PY{p}{)} \PY{c+c1}{\PYZsh{} Change with your triskele install path}
\PY{k+kn}{import} \PY{n+nn}{triskele}
\end{Verbatim}
\hypertarget{simple-test}{%
\section{Simple test}\label{simple-test}}
Just to test the python import and the triskele compilation. Create an
array with random values.
\begin{Verbatim}[commandchars=\\\{\}]
{\color{incolor}In [{\color{incolor}3}]:} \PY{n}{raster} \PY{o}{=} \PY{n}{np}\PY{o}{.}\PY{n}{random}\PY{o}{.}\PY{n}{randint}\PY{p}{(}\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{255}\PY{p}{,} \PY{p}{(}\PY{l+m+mi}{1000}\PY{p}{,} \PY{l+m+mi}{1000}\PY{p}{)}\PY{p}{)}
\end{Verbatim}
Create the instance of triskele object with the array to filter.
\begin{Verbatim}[commandchars=\\\{\}]
{\color{incolor}In [{\color{incolor}4}]:} \PY{n}{ap\PYZus{}random} \PY{o}{=} \PY{n}{triskele}\PY{o}{.}\PY{n}{Triskele}\PY{p}{(}\PY{n}{raster}\PY{p}{)}
\end{Verbatim}
Filter the raster with this method :
\begin{Verbatim}[commandchars=\\\{\}]
{\color{incolor}In [{\color{incolor}5}]:} \PY{n}{filtered\PYZus{}raster} \PY{o}{=} \PY{n}{ap\PYZus{}random}\PY{o}{.}\PY{n}{filter}\PY{p}{(}\PY{p}{)}
\PY{n+nb}{print}\PY{p}{(}\PY{n}{filtered\PYZus{}raster}\PY{p}{)}
\end{Verbatim}
\begin{Verbatim}[commandchars=\\\{\}]
STDOUT:
Input:/tmp/infile.tif [1000,1000] (1 chanels of Byte)
Crop topLeft:(0,0) size:[1000,1000] band:[0]
Output:/tmp/outfile.tif
core count:24
STDERR:
*** apGenerator done!
Leaf Count Mean Min Max
max 1 1000000 1000000 1000000
Comp Count Mean Min Max
max 1 347157 347157 347157
Time Sum Count Mean Min Max
build tree 00:00:00.086426 1 00:00:00.086426 00:00:00.086426 00:00:00.086426
setup 00:00:00.038795 1 00:00:00.038795 00:00:00.038795 00:00:00.038795
parents 00:00:00.026013 1 00:00:00.026013 00:00:00.026013 00:00:00.026013
merge 00:00:00.005112 1 00:00:00.005112 00:00:00.005112 00:00:00.005112
index 00:00:00.006034 1 00:00:00.006034 00:00:00.006034 00:00:00.006034
compress 00:00:00.008091 1 00:00:00.008091 00:00:00.008091 00:00:00.008091
children 00:00:00.108123 1 00:00:00.108123 00:00:00.108123 00:00:00.108123
area 00:00:00.004463 1 00:00:00.004463 00:00:00.004463 00:00:00.004463
[[194 125 164 {\ldots} 213 182 168]
[ 25 48 242 {\ldots} 229 16 143]
[ 7 207 212 {\ldots} 112 94 100]
{\ldots}
[227 180 185 {\ldots} 34 10 229]
[129 149 27 {\ldots} 217 66 195]
[145 203 251 {\ldots} 186 116 165]]
\end{Verbatim}
If there is no error so far everything is okay.
Theoretically the filtered output should be a copy of the input. (it's
not, there is something going on along the edges, we should tell
François !)
\hypertarget{attribute-filter}{%
\subsection{Attribute filter}\label{attribute-filter}}
With \texttt{filter} method you can ask for specific attributes
filtering. The results are stacked.
\begin{Verbatim}[commandchars=\\\{\}]
{\color{incolor}In [{\color{incolor}6}]:} \PY{n}{filtered\PYZus{}rasters} \PY{o}{=} \PY{n}{ap\PYZus{}random}\PY{o}{.}\PY{n}{filter}\PY{p}{(}\PY{n}{area}\PY{o}{=}\PY{p}{[}\PY{l+m+mi}{10}\PY{p}{,} \PY{l+m+mi}{100}\PY{p}{,} \PY{l+m+mf}{1e3}\PY{p}{]}\PY{p}{)}
\PY{n+nb}{print}\PY{p}{(}\PY{n}{filtered\PYZus{}rasters}\PY{o}{.}\PY{n}{shape}\PY{p}{)}
\end{Verbatim}
\begin{Verbatim}[commandchars=\\\{\}]
STDOUT:
Input:/tmp/infile.tif [1000,1000] (1 chanels of Byte)
Crop topLeft:(0,0) size:[1000,1000] band:[0]
Output:/tmp/outfile.tif
core count:24
STDERR:
*** apGenerator done!
Leaf Count Mean Min Max
max 1 1000000 1000000 1000000
Comp Count Mean Min Max
max 1 347157 347157 347157
Time Sum Count Mean Min Max
build tree 00:00:00.094336 1 00:00:00.094336 00:00:00.094336 00:00:00.094336
setup 00:00:00.040030 1 00:00:00.040030 00:00:00.040030 00:00:00.040030
parents 00:00:00.034337 1 00:00:00.034337 00:00:00.034337 00:00:00.034337
merge 00:00:00.004679 1 00:00:00.004679 00:00:00.004679 00:00:00.004679
index 00:00:00.005624 1 00:00:00.005624 00:00:00.005624 00:00:00.005624
compress 00:00:00.007514 1 00:00:00.007514 00:00:00.007514 00:00:00.007514
children 00:00:00.113179 1 00:00:00.113179 00:00:00.113179 00:00:00.113179
area 00:00:00.004390 1 00:00:00.004390 00:00:00.004390 00:00:00.004390
filtering 00:00:00.007614 1 00:00:00.007614 00:00:00.007614 00:00:00.007614
(1000, 1000, 4)
\end{Verbatim}
I don't know if you can see the signature with yout IDE :
\texttt{filter(tree=\textquotesingle{}max-tree\textquotesingle{},\ area=None,\ standard\_deviation=None,\ moment\_of\_inertia=None)}
The parameters are the same than apGenerator from François, I should
write the doc\ldots{}
\hypertarget{advanced-test}{%
\section{Advanced test}\label{advanced-test}}
With a real example.
Load a TIFF file with 4 bands :
\begin{Verbatim}[commandchars=\\\{\}]
{\color{incolor}In [{\color{incolor}7}]:} \PY{n}{raster} \PY{o}{=} \PY{n}{triskele}\PY{o}{.}\PY{n}{read}\PY{p}{(}\PY{l+s+s1}{\PYZsq{}}\PY{l+s+s1}{../triskele/data/10m.tif}\PY{l+s+s1}{\PYZsq{}}\PY{p}{)}
\PY{n+nb}{print}\PY{p}{(}\PY{n}{raster}\PY{o}{.}\PY{n}{shape}\PY{p}{)}
\PY{n}{plt}\PY{o}{.}\PY{n}{imshow}\PY{p}{(}\PY{n}{raster}\PY{p}{)}
\PY{n}{plt}\PY{o}{.}\PY{n}{show}\PY{p}{(}\PY{p}{)}
\end{Verbatim}
\begin{Verbatim}[commandchars=\\\{\}]
(291, 685, 4)
\end{Verbatim}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_14_1.png}
\end{center}
{ \hspace*{\fill} \\}
Create the instance with some tweaks.
\begin{Verbatim}[commandchars=\\\{\}]
{\color{incolor}In [{\color{incolor}8}]:} \PY{n}{ap\PYZus{}10m} \PY{o}{=} \PY{n}{triskele}\PY{o}{.}\PY{n}{Triskele}\PY{p}{(}\PY{n}{raster}\PY{p}{,} \PY{n}{verbose}\PY{o}{=}\PY{k+kc}{False}\PY{p}{)} \PY{c+c1}{\PYZsh{} Verbose False to hide apGenerator gibbering :p}
\end{Verbatim}
Filter as you wish :
\begin{Verbatim}[commandchars=\\\{\}]
{\color{incolor}In [{\color{incolor}9}]:} \PY{n}{tree} \PY{o}{=} \PY{l+s+s1}{\PYZsq{}}\PY{l+s+s1}{tos\PYZhy{}tree}\PY{l+s+s1}{\PYZsq{}} \PY{c+c1}{\PYZsh{} \PYZob{}\PYZsq{}max\PYZhy{}tree\PYZsq{}, \PYZsq{}min\PYZhy{}tree\PYZsq{}, ...\PYZcb{}}
\PY{n}{area} \PY{o}{=} \PY{p}{[}\PY{l+m+mi}{100}\PY{p}{,} \PY{l+m+mf}{1e3}\PY{p}{,} \PY{l+m+mf}{1e4}\PY{p}{]}
\PY{n}{standard\PYZus{}deviation} \PY{o}{=} \PY{k+kc}{None}
\PY{n}{moment\PYZus{}of\PYZus{}inertia} \PY{o}{=} \PY{k+kc}{None}
\PY{n}{filtered\PYZus{}rasters} \PY{o}{=} \PY{n}{ap\PYZus{}10m}\PY{o}{.}\PY{n}{filter}\PY{p}{(}\PY{n}{tree}\PY{o}{=}\PY{n}{tree}\PY{p}{,} \PY{n}{area}\PY{o}{=}\PY{n}{area}\PY{p}{,} \PY{n}{standard\PYZus{}deviation}\PY{o}{=}\PY{n}{standard\PYZus{}deviation}\PY{p}{,} \PY{n}{moment\PYZus{}of\PYZus{}inertia}\PY{o}{=}\PY{n}{moment\PYZus{}of\PYZus{}inertia}\PY{p}{)}
\PY{n+nb}{print}\PY{p}{(}\PY{n}{filtered\PYZus{}rasters}\PY{o}{.}\PY{n}{shape}\PY{p}{)}
\end{Verbatim}
\begin{Verbatim}[commandchars=\\\{\}]
(291, 685, 16)
\end{Verbatim}
\hypertarget{display-results}{%
\subsection{Display results}\label{display-results}}
\begin{Verbatim}[commandchars=\\\{\}]
{\color{incolor}In [{\color{incolor}10}]:} \PY{k}{for} \PY{n}{i} \PY{o+ow}{in} \PY{n+nb}{range}\PY{p}{(}\PY{n}{filtered\PYZus{}rasters}\PY{o}{.}\PY{n}{shape}\PY{p}{[}\PY{l+m+mi}{2}\PY{p}{]}\PY{p}{)}\PY{p}{:}
\PY{n}{plt}\PY{o}{.}\PY{n}{imshow}\PY{p}{(}\PY{n}{filtered\PYZus{}rasters}\PY{p}{[}\PY{p}{:}\PY{p}{,}\PY{p}{:}\PY{p}{,}\PY{n}{i}\PY{p}{]}\PY{p}{)}
\PY{n}{plt}\PY{o}{.}\PY{n}{show}\PY{p}{(}\PY{p}{)}
\end{Verbatim}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_0.png}
\end{center}
{ \hspace*{\fill} \\}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_1.png}
\end{center}
{ \hspace*{\fill} \\}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_2.png}
\end{center}
{ \hspace*{\fill} \\}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_3.png}
\end{center}
{ \hspace*{\fill} \\}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_4.png}
\end{center}
{ \hspace*{\fill} \\}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_5.png}
\end{center}
{ \hspace*{\fill} \\}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_6.png}
\end{center}
{ \hspace*{\fill} \\}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_7.png}
\end{center}
{ \hspace*{\fill} \\}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_8.png}
\end{center}
{ \hspace*{\fill} \\}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_9.png}
\end{center}
{ \hspace*{\fill} \\}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_10.png}
\end{center}
{ \hspace*{\fill} \\}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_11.png}
\end{center}
{ \hspace*{\fill} \\}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_12.png}
\end{center}
{ \hspace*{\fill} \\}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_13.png}
\end{center}
{ \hspace*{\fill} \\}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_14.png}
\end{center}
{ \hspace*{\fill} \\}
\begin{center}
\adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_20_15.png}
\end{center}
{ \hspace*{\fill} \\}
% Add a bibliography block to the postdoc
\end{document}