538 lines
15 KiB
Plaintext
538 lines
15 KiB
Plaintext
{
|
|
"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": 5,
|
|
"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",
|
|
"alpha = ((dtm_hs.astype(float) - 127) / 127) ** 2 * .4\n",
|
|
"cX, cY = np.meshgrid(np.arange(dtm.shape[1]), np.arange(dtm.shape[0]))\n",
|
|
"kwargs = {}\n",
|
|
"kwargs['vmin'], kwargs['vmax'] = np.quantile(dtm, (.01, .99))\n",
|
|
"\n",
|
|
"model = {}"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Interactive GUI"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 13,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"application/vnd.jupyter.widget-view+json": {
|
|
"model_id": "6a5848b7f64b487fae0ecb76d2ff5f72",
|
|
"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=True, highlight=True)>"
|
|
]
|
|
},
|
|
"execution_count": 13,
|
|
"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",
|
|
" # 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, **kwargs)\n",
|
|
" else:\n",
|
|
" #model['ax_dtm'].imshow(dtm_hs, cmap=plt.cm.Greys)\n",
|
|
" model['ax_dtm'].imshow(dtm, **kwargs)\n",
|
|
" model['ax_dtm'].imshow(dtm_hs, cmap=plt.cm.Greys, alpha=alpha)\n",
|
|
"\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",
|
|
" model['ax_dtm'].clear()\n",
|
|
" model['ax_dtm'].imshow(dtm, **kwargs)\n",
|
|
" model['ax_dtm'].imshow(dtm_hs, cmap=plt.cm.Greys, alpha=alpha)\n",
|
|
" model['ax_dtm'].contour(cX, cY, img, [model['tree_cache']._alt[-1]], colors='lime')\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=True, 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",
|
|
"#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": 70,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"(2001, 2001)"
|
|
]
|
|
},
|
|
"execution_count": 70,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"model['tree_cache'].reconstruct().shape"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 73,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"2001"
|
|
]
|
|
},
|
|
"execution_count": 73,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"len(X)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 66,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"16.93000030517578"
|
|
]
|
|
},
|
|
"execution_count": 66,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"model['tree_cache']._alt[-1]"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 62,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"ename": "NameError",
|
|
"evalue": "name 'img' is not defined",
|
|
"output_type": "error",
|
|
"traceback": [
|
|
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
|
|
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
|
|
"\u001b[0;32m<ipython-input-62-77a7949ba908>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mimg\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
|
|
"\u001b[0;31mNameError\u001b[0m: name 'img' is not defined"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"img"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 82,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"<matplotlib.contour.QuadContourSet at 0x7fe61585e430>"
|
|
]
|
|
},
|
|
"execution_count": 82,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"plt.close()\n",
|
|
"plt.figure()\n",
|
|
"ax = plt.gca()\n",
|
|
"\n",
|
|
"plt.imshow(dtm, **kwargs)\n",
|
|
"plt.colorbar()\n",
|
|
"ax.imshow(dtm_hs, cmap=plt.cm.Greys, alpha=alpha)\n",
|
|
"ax.contour(X, Y, dtm, [30], colors='lime')"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 21,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"(2001, 2001)"
|
|
]
|
|
},
|
|
"execution_count": 21,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"dtm_hs.shape"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 22,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"(2001, 2001)"
|
|
]
|
|
},
|
|
"execution_count": 22,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"dtm.shape"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 37,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 42,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"0.12586920037702484"
|
|
]
|
|
},
|
|
"execution_count": 42,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"alpha.mean()"
|
|
]
|
|
},
|
|
{
|
|
"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
|
|
}
|