462 lines
12 KiB
Plaintext
462 lines
12 KiB
Plaintext
{
|
|
"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
|
|
}
|