initial commit
This commit is contained in:
commit
3810964728
353
demo.ipynb
Normal file
353
demo.ipynb
Normal file
@ -0,0 +1,353 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Interactive pixel to Spectrum x DTM Analysis in Attribute Space\n",
|
||||
"\n",
|
||||
"Fonctionnel, mériterais une réécriture :p"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 0 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"import numpy as np\n",
|
||||
"import rasterio as rio\n",
|
||||
"import sap\n",
|
||||
"import matplotlib as mpl\n",
|
||||
"from matplotlib import pyplot as plt\n",
|
||||
"import inspect\n",
|
||||
"import ipywidgets as ipw\n",
|
||||
"import higra as hg\n",
|
||||
"\n",
|
||||
"plt.style.use('dark_background')\n",
|
||||
"plt.set_cmap('plasma')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Load DTM"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(2001, 2001)"
|
||||
]
|
||||
},
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"dtm = rio.open('data/dsm_vox_50cm_tile_-12_0.tif').read()[0]#[-500:,-500:].copy()\n",
|
||||
"dtm.shape"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Setup utils"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def ui8_clip(x, vmin=None, vmax=None):\n",
|
||||
" vmin = vmin if vmin else np.nanmin(x)\n",
|
||||
" vmax = vmax if vmax else np.nanmax(x)\n",
|
||||
" \n",
|
||||
" ui8 = ((np.clip(x, vmin, vmax) - vmin) / (vmax - vmin) * 255).astype(np.uint8)\n",
|
||||
" \n",
|
||||
" return np.ma.array(ui8, mask=np.isnan(x))\n",
|
||||
"\n",
|
||||
"def hillshades(x):\n",
|
||||
" hs = np.zeros_like(x) \n",
|
||||
" hs[:-1,:-1] = x[:-1,:-1] - x[1:,1:]\n",
|
||||
" return hs\n",
|
||||
"\n",
|
||||
"from sap.spectra import get_bins\n",
|
||||
"\n",
|
||||
"def spectrum2d(tree, x_attribute, y_attribute, x_count=100, y_count=100, \n",
|
||||
" x_log=False, y_log=False, weighted=True, normalized=True,\n",
|
||||
" node_mask=None):\n",
|
||||
" x = tree.get_attribute(x_attribute)\n",
|
||||
" y = tree.get_attribute(y_attribute)\n",
|
||||
"\n",
|
||||
" bins = (get_bins(x, x_count, 'geo' if x_log else 'lin'),\n",
|
||||
" get_bins(y, y_count, 'geo' if y_log else 'lin'))\n",
|
||||
"\n",
|
||||
" weights = tree.get_attribute('area') if weighted else None\n",
|
||||
" weights = weights / tree._image.size if normalized and weighted else weights\n",
|
||||
"\n",
|
||||
" s, xedges, yedges = np.histogram2d(x[node_mask].ravel(), y[node_mask].ravel(),\n",
|
||||
" bins=bins, density=None, weights=weights[node_mask].ravel())\n",
|
||||
"\n",
|
||||
" return s, xedges, yedges, x_log, y_log\n",
|
||||
"\n",
|
||||
"def pixel_to_node(tree, mask):\n",
|
||||
" \"\"\"Compute the node mask from the pixel mask.\"\"\"\n",
|
||||
" node_mask = hg.accumulate_and_min_sequential(tree._tree, \n",
|
||||
" np.ones(tree._tree.num_vertices(), dtype=np.uint8), \n",
|
||||
" mask.ravel(), \n",
|
||||
" hg.Accumulators.min).astype(np.bool)\n",
|
||||
" return node_mask\n",
|
||||
"\n",
|
||||
"dtm_hs = ui8_clip(hillshades(dtm), -1, 1)\n",
|
||||
"cm_hud = mpl.colors.LinearSegmentedColormap.from_list('GreenHUD', [(0.,0.,0.,0.), (0.,1.,0.,1.)], 256)\n",
|
||||
"\n",
|
||||
"model = {}"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Interactive GUI"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"application/vnd.jupyter.widget-view+json": {
|
||||
"model_id": "87b90d0f6d584b83af83965b32f15118",
|
||||
"version_major": 2,
|
||||
"version_minor": 0
|
||||
},
|
||||
"text/plain": [
|
||||
"interactive(children=(Dropdown(description='tree', index=1, options=(('AlphaTree', <class 'sap.trees.AlphaTree…"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<function __main__.spectrum_widget(tree=<class 'sap.trees.MaxTree'>, x='area', y='compactness', x_log=True, y_log=False, hillshade_overlay=False, highlight=True)>"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# \"widget\" or \"qt\", qt provide a snappier GUI\n",
|
||||
"#%matplotlib widget\n",
|
||||
"%matplotlib qt\n",
|
||||
"plt.style.use('dark_background')\n",
|
||||
"\n",
|
||||
"model.update({\n",
|
||||
" 'ax_spectrum': None,\n",
|
||||
" 'ax_dtm': None,\n",
|
||||
" 'x': None,\n",
|
||||
" 'y': None,\n",
|
||||
" 'x1': 0,\n",
|
||||
" 'x2': 0,\n",
|
||||
" 'y1': 0,\n",
|
||||
" 'y2': 0,\n",
|
||||
" 'dtm_x1': 0,\n",
|
||||
" 'dtm_x2': 0,\n",
|
||||
" 'dtm_y1': 0,\n",
|
||||
" 'dtm_y2': 0,\n",
|
||||
" 'overlay': False,\n",
|
||||
" 'pixel_mask': None,\n",
|
||||
"})\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"def line_select_callback(eclick, erelease):\n",
|
||||
" model['x1'], model['y1'] = eclick.xdata, eclick.ydata\n",
|
||||
" model['x2'], model['y2'] = erelease.xdata, erelease.ydata\n",
|
||||
" filter_dtm()\n",
|
||||
" \n",
|
||||
"def dtm_select_callback(eclick, erelease):\n",
|
||||
" model['dtm_y1'], model['dtm_x1'] = np.int(eclick.xdata), np.int(eclick.ydata)\n",
|
||||
" model['dtm_y2'], model['dtm_x2'] = np.int(erelease.xdata), np.int(erelease.ydata)\n",
|
||||
" filter_nodes()\n",
|
||||
" \n",
|
||||
"def filter_nodes():\n",
|
||||
" print('HI!')\n",
|
||||
" x1, x2 = model['dtm_x1'], model['dtm_x2']\n",
|
||||
" y1, y2 = model['dtm_y1'], model['dtm_y2']\n",
|
||||
" pmask = np.zeros_like(dtm, dtype=np.bool)\n",
|
||||
" pmask[x1:x2,y1:y2] = True\n",
|
||||
" model['pixel_mask'] = pmask\n",
|
||||
" show_spectrum()\n",
|
||||
" \n",
|
||||
"def show_spectrum():\n",
|
||||
" x, y = model['x'], model['y']\n",
|
||||
" x_log, y_log = model['x_log'], model['y_log']\n",
|
||||
" highlight = model['highlight']\n",
|
||||
"\n",
|
||||
" # Compute spectrum\n",
|
||||
" ps = sap.spectrum2d(model['tree_cache'], x, y, 200, 100, x_log, y_log)\n",
|
||||
" \n",
|
||||
" if highlight:\n",
|
||||
" node_mask = pixel_to_node(model['tree_cache'], model['pixel_mask']) if model['pixel_mask'] is not None else None\n",
|
||||
" ps_mask = spectrum2d(model['tree_cache'], x, y, 200, 100, x_log, y_log, node_mask=node_mask)\n",
|
||||
" ps[0][:] = ps_mask[0] / ps[0]\n",
|
||||
"\n",
|
||||
" \n",
|
||||
" # Display spectrum\n",
|
||||
" plt.sca(model['ax_spectrum'])\n",
|
||||
" plt.cla()\n",
|
||||
" sap.show_spectrum(*ps, log_scale=not highlight)\n",
|
||||
" plt.xlabel(x)\n",
|
||||
" plt.ylabel(y)\n",
|
||||
" plt.grid(which='minor', color='#BBBBBB', linestyle=':')\n",
|
||||
" plt.grid(True)\n",
|
||||
" \n",
|
||||
" plt.gcf().canvas.draw()\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"def filter_dtm():\n",
|
||||
" x1, x2 = model['x1'], model['x2']\n",
|
||||
" y1, y2 = model['y1'], model['y2']\n",
|
||||
" \n",
|
||||
" # No filter\n",
|
||||
" if x1 == x2 == y1 == y2 == 0:\n",
|
||||
" if not model['overlay']:\n",
|
||||
" model['ax_dtm'].imshow(dtm)\n",
|
||||
" else:\n",
|
||||
" model['ax_dtm'].imshow(dtm_hs, cmap=plt.cm.Greys)\n",
|
||||
" return\n",
|
||||
" \n",
|
||||
" # Filter\n",
|
||||
" t = model['tree_cache']\n",
|
||||
" X = t.get_attribute(model['x'])\n",
|
||||
" Y = t.get_attribute(model['y'])\n",
|
||||
"\n",
|
||||
" img = model['tree_cache'].reconstruct((X < x1) | (X > x2) | (Y < y1) | (Y > y2), filtering='subtractive')\n",
|
||||
" \n",
|
||||
" if not model['overlay']:\n",
|
||||
" model['ax_dtm'].imshow(img)\n",
|
||||
" else:\n",
|
||||
" over = img != model['tree_cache']._alt[-1]\n",
|
||||
" model['ax_dtm'].imshow(dtm_hs, cmap=plt.cm.Greys)\n",
|
||||
" model['ax_dtm'].imshow(over, cmap=cm_hud, alpha=.4)\n",
|
||||
" \n",
|
||||
" plt.gcf().canvas.draw()\n",
|
||||
"\n",
|
||||
"def spectrum_widget(tree=sap.MaxTree, x='area', y='compactness', x_log=True, y_log=False,\n",
|
||||
" hillshade_overlay=False, highlight=True):\n",
|
||||
" # Update only toggle of hillshade_overlay\n",
|
||||
" if hillshade_overlay != model['overlay']:\n",
|
||||
" model['overlay'] = hillshade_overlay\n",
|
||||
" filter_dtm()\n",
|
||||
" return\n",
|
||||
" \n",
|
||||
" # Update model and tree computation if needed\n",
|
||||
" model.update({'x': x, 'y': y, 'x_log': x_log, 'y_log': y_log, 'highlight': highlight})\n",
|
||||
" if not 'tree_cache' in model or not isinstance(model['tree_cache'], tree):\n",
|
||||
" print('Computing the {} of the DTM...'.format(tree))\n",
|
||||
" model['tree_cache'] = tree(dtm)\n",
|
||||
" \n",
|
||||
" show_spectrum()\n",
|
||||
" \n",
|
||||
" # Setup selector spectrum\n",
|
||||
" selector.rs = mpl.widgets.RectangleSelector(model['ax_spectrum'], line_select_callback,\n",
|
||||
" drawtype='box', useblit=True,\n",
|
||||
" button=[1, 3], # don't use middle button\n",
|
||||
" minspanx=5, minspany=5,\n",
|
||||
" spancoords='pixels',\n",
|
||||
" rectprops = dict(facecolor=(0,1,0,.2), edgecolor=(0,1,0), fill=True, ls='--', lw=2),\n",
|
||||
" interactive=True)\n",
|
||||
" \n",
|
||||
" # Setup selector DTM\n",
|
||||
" selector_dtm.rs = mpl.widgets.RectangleSelector(model['ax_dtm'], dtm_select_callback,\n",
|
||||
" drawtype='box', useblit=True,\n",
|
||||
" button=[1, 3], # don't use middle button\n",
|
||||
" minspanx=5, minspany=5,\n",
|
||||
" spancoords='pixels',\n",
|
||||
" rectprops = dict(facecolor=(1,1,0,.2), edgecolor=(1,1,0), fill=True, ls='--', lw=2),\n",
|
||||
" interactive=True)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"plt.close()\n",
|
||||
"fig = plt.figure('Spectrum', figsize=(15, 4), constrained_layout=True)\n",
|
||||
"grid = mpl.gridspec.GridSpec(1, 3, fig)\n",
|
||||
"\n",
|
||||
"model['ax_spectrum'] = fig.add_subplot(grid[:2])\n",
|
||||
"model['ax_dtm'] = fig.add_subplot(grid[-1])\n",
|
||||
"\n",
|
||||
"# Selectors\n",
|
||||
"def selector(event):\n",
|
||||
" pass\n",
|
||||
"\n",
|
||||
"def selector_dtm(event):\n",
|
||||
" pass\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"plt.connect('key_press_event', selector)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"#model['ax_dtm'].imshow(dtm_hs, cmap=plt.cm.Greys)\n",
|
||||
"filter_dtm()\n",
|
||||
"\n",
|
||||
"ipw.interact(spectrum_widget, \n",
|
||||
" tree=inspect.getmembers(sap.trees, lambda t: inspect.isclass(t) and issubclass(t, sap.Tree) and t != sap.Tree),\n",
|
||||
" x=sap.available_attributes().keys(),\n",
|
||||
" y=sap.available_attributes().keys())"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"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",
|
||||
"version": "3.9.2"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 4
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user