{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Blockwise coregistration\n\nOften, biases are spatially variable, and a \"global\" shift may not be enough to coregister a DEM properly.\nIn the `sphx_glr_auto_examples_plot_nuth_kaab.py` example, we saw that the method improved the alignment significantly, but there were still possibly nonlinear artefacts in the result.\nClearly, nonlinear coregistration approaches are needed.\nOne solution is :class:`xdem.coreg.BlockwiseCoreg`, a helper to run any ``Coreg`` class over an arbitrarily small grid, and then \"puppet warp\" the DEM to fit the reference best.\n\nThe ``BlockwiseCoreg`` class runs in five steps:\n\n1. Generate a subdivision grid to divide the DEM in N blocks.\n2. Run the requested coregistration approach in each block.\n3. Extract each result as a source and destination X/Y/Z point.\n4. Interpolate the X/Y/Z point-shifts into three shift-rasters.\n5. Warp the DEM to apply the X/Y/Z shifts.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# sphinx_gallery_thumbnail_number = 2\nimport matplotlib.pyplot as plt\nimport geoutils as gu\nimport numpy as np\nimport xdem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Example files**\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "xdem.examples.download_longyearbyen_examples()\n\nreference_dem = xdem.DEM(xdem.examples.get_path(\"longyearbyen_ref_dem\"), silent=True)\ndem_to_be_aligned = xdem.DEM(xdem.examples.get_path(\"longyearbyen_tba_dem\"), silent=True)\nglacier_outlines = gu.Vector(xdem.examples.get_path(\"longyearbyen_glacier_outlines\"))\n\n# Create a stable ground mask (not glacierized) to mark \"inlier data\"\ninlier_mask = ~glacier_outlines.create_mask(reference_dem)\n\nplt_extent = [\n reference_dem.bounds.left,\n reference_dem.bounds.right,\n reference_dem.bounds.bottom,\n reference_dem.bounds.top,\n]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The DEM to be aligned (a 1990 photogrammetry-derived DEM) has some vertical and horizontal biases that we want to avoid, as well as possible nonlinear distortions.\nThe product is a mosaic of multiple DEMs, so \"seams\" may exist in the data.\nThese can be visualized by plotting a change map:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "diff_before = (reference_dem - dem_to_be_aligned).data\n\nplt.figure(figsize=(8, 5))\nplt.imshow(diff_before.squeeze(), cmap=\"coolwarm_r\", vmin=-10, vmax=10, extent=plt_extent)\nplt.colorbar()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Horizontal and vertical shifts can be estimated using :class:`xdem.coreg.NuthKaab`.\nLet's prepare a coregistration class that calculates 64 offsets, evenly spread over the DEM.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "blockwise = xdem.coreg.BlockwiseCoreg(xdem.coreg.NuthKaab(), subdivision=64)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The grid that will be used can be visualized with a helper function.\nCoregistration will be performed in each block separately.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.title(\"Subdivision grid\")\nplt.imshow(blockwise.subdivide_array(dem_to_be_aligned.shape), cmap=\"gist_ncar\")\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Coregistration is performed with the ``.fit()`` method.\nThis runs in multiple threads by default, so more CPU cores are preferable here.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "blockwise.fit(reference_dem.data, dem_to_be_aligned.data, transform=reference_dem.transform, inlier_mask=inlier_mask)\n\naligned_dem_data = blockwise.apply(dem_to_be_aligned.data, transform=dem_to_be_aligned.transform)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The estimated shifts can be visualized by applying the coregistration to a completely flat surface.\nThis shows the estimated shifts that would be applied in elevation; additional horizontal shifts will also be applied if the method supports it.\nThe :func:`xdem.coreg.BlockwiseCoreg.stats` method can be used to annotate each block with its associated Z shift.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "z_correction = blockwise.apply(np.zeros_like(dem_to_be_aligned.data), transform=dem_to_be_aligned.transform)\nplt.title(\"Vertical correction\")\nplt.imshow(z_correction, cmap=\"coolwarm_r\", vmin=-10, vmax=10, extent=plt_extent)\nfor _, row in blockwise.stats().iterrows():\n plt.annotate(round(row[\"z_off\"], 1), (row[\"center_x\"], row[\"center_y\"]), ha\n=\"center\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, the new difference can be plotted to validate that it improved.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "diff_after = reference_dem.data - aligned_dem_data\n\nplt.figure(figsize=(8, 5))\nplt.imshow(diff_after.squeeze(), cmap=\"coolwarm_r\", vmin=-10, vmax=10, extent=plt_extent)\nplt.colorbar()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can compare the `spatial_stats_nmad` to validate numerically that there was an improvment:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(f\"Error before: {xdem.spatial_tools.nmad(diff_before):.2f} m\")\nprint(f\"Error after: {xdem.spatial_tools.nmad(diff_after):.2f} m\")" ] } ], "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.7.9" } }, "nbformat": 4, "nbformat_minor": 0 }