This commit is contained in:
Florent Guiotte 2021-03-09 21:16:54 +01:00
parent e75aa11cc1
commit f37c9da47a
2 changed files with 66 additions and 200 deletions

View File

@ -25,7 +25,6 @@
}
],
"source": [
"import numpy as np\n",
"import rasterio as rio\n",
"import sap\n",
"import matplotlib as mpl\n",
@ -33,6 +32,7 @@
"import inspect\n",
"import ipywidgets as ipw\n",
"import higra as hg\n",
"import numpy as np\n",
"\n",
"plt.style.use('dark_background')\n",
"plt.set_cmap('plasma')"
@ -75,8 +75,13 @@
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"execution_count": 3,
"metadata": {
"jupyter": {
"source_hidden": true
},
"tags": []
},
"outputs": [],
"source": [
"def ui8_clip(x, vmin=None, vmax=None):\n",
@ -138,13 +143,18 @@
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"execution_count": 4,
"metadata": {
"jupyter": {
"source_hidden": true
},
"tags": []
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "6a5848b7f64b487fae0ecb76d2ff5f72",
"model_id": "81086011013945a497e4ffb7163998f9",
"version_major": 2,
"version_minor": 0
},
@ -158,10 +168,10 @@
{
"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)>"
"<function __main__.spectrum_widget(tree=<class 'sap.trees.MaxTree'>, x='area', y='compactness', x_log=True, y_log=False, hillshade_overlay=True, highlight=False)>"
]
},
"execution_count": 13,
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
@ -185,7 +195,7 @@
" 'dtm_x2': 0,\n",
" 'dtm_y1': 0,\n",
" 'dtm_y2': 0,\n",
" 'overlay': False,\n",
" 'overlay': True,\n",
" 'pixel_mask': None,\n",
"})\n",
"\n",
@ -207,6 +217,7 @@
" pmask = np.zeros_like(dtm, dtype=np.bool)\n",
" pmask[x1:x2,y1:y2] = True\n",
" model['pixel_mask'] = pmask\n",
" model['highlight'] = True\n",
" show_spectrum()\n",
" \n",
"def show_spectrum():\n",
@ -243,7 +254,6 @@
" 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",
@ -256,13 +266,10 @@
"\n",
" img = model['tree_cache'].reconstruct((X < x1) | (X > x2) | (Y < y1) | (Y > y2), filtering='subtractive')\n",
" \n",
" model['ax_dtm'].clear()\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",
@ -270,7 +277,7 @@
" 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",
" hillshade_overlay=True, highlight=False):\n",
" # Update only toggle of hillshade_overlay\n",
" if hillshade_overlay != model['overlay']:\n",
" model['overlay'] = hillshade_overlay\n",
@ -305,7 +312,7 @@
"\n",
"\n",
"plt.close()\n",
"fig = plt.figure('Spectrum', figsize=(15, 4), constrained_layout=True)\n",
"fig = plt.figure('Spectrum', figsize=(18, 5), constrained_layout=True)\n",
"grid = mpl.gridspec.GridSpec(1, 3, fig)\n",
"\n",
"model['ax_spectrum'] = fig.add_subplot(grid[:2])\n",
@ -319,198 +326,16 @@
" pass\n",
"\n",
"\n",
"plt.connect('key_press_event', selector)\n",
"#plt.connect('key_press_event', selector)\n",
"\n",
"\n",
"#filter_dtm()\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": {

View File

@ -0,0 +1,41 @@
import numpy as np
def ui8_clip(x, vmin=None, vmax=None):
vmin = vmin if vmin else np.nanmin(x)
vmax = vmax if vmax else np.nanmax(x)
ui8 = ((np.clip(x, vmin, vmax) - vmin) / (vmax - vmin) * 255).astype(np.uint8)
return np.ma.array(ui8, mask=np.isnan(x))
def hillshades(x):
hs = np.zeros_like(x)
hs[:-1,:-1] = x[:-1,:-1] - x[1:,1:]
return hs
from sap.spectra import get_bins
def spectrum2d(tree, x_attribute, y_attribute, x_count=100, y_count=100,
x_log=False, y_log=False, weighted=True, normalized=True,
node_mask=None):
x = tree.get_attribute(x_attribute)
y = tree.get_attribute(y_attribute)
bins = (get_bins(x, x_count, 'geo' if x_log else 'lin'),
get_bins(y, y_count, 'geo' if y_log else 'lin'))
weights = tree.get_attribute('area') if weighted else None
weights = weights / tree._image.size if normalized and weighted else weights
s, xedges, yedges = np.histogram2d(x[node_mask].ravel(), y[node_mask].ravel(),
bins=bins, density=None, weights=weights[node_mask].ravel())
return s, xedges, yedges, x_log, y_log
def pixel_to_node(tree, mask):
"""Compute the node mask from the pixel mask."""
node_mask = hg.accumulate_and_min_sequential(tree._tree,
np.ones(tree._tree.num_vertices(), dtype=np.uint8),
mask.ravel(),
hg.Accumulators.min).astype(np.bool)
return node_mask