{
    "cells":[
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "# Getting started",
                "\n",
                "In this tutorial we will learn the basics of raster and ",
                "point cloud processing using Pyoints."
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "We begin with loading the required modules."
            ]
        },
        {
            "cell_type":"code",
            "execution_count":null,
            "metadata":{},
            "outputs":[],
            "source":[
                "import numpy as np\n",
                "\n",
                "from pyoints import (\n",
                "\tnptools,\n",
                "\tProj,\n",
                "\tGeoRecords,\n",
                "\tGrid,\n",
                "\tExtent,\n",
                "\ttransformation,\n",
                "\tfilters,\n",
                "\tclustering,\n",
                "\tclassification,\n",
                "\tsmoothing,\n",
                ")"
            ]
        },
        {
            "cell_type":"code",
            "execution_count":null,
            "metadata":{},
            "outputs":[

            ],
            "source":[
                "from mpl_toolkits import mplot3d\n",
                "import matplotlib.pyplot as plt\n",
                "\n",
                "%matplotlib inline"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "## Grid\n",
                "We create a two dimensional raster by providing a ",
                "projection system, a transformation matrix and some data. ",
                "The transformation matrix defines the origin and scale of ",
                "the raster. For this example, we just use the default ",
                "projection system. "
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "print('projection')\n",
                "proj = Proj()\n",
                "print(proj)\n",
                "\n",
                "print('numpy record array:')\n",
                "rec = nptools.recarray(\n",
                    "\t{'indices': np.mgrid[0:50, 0:30].T},\n",
                    "\tdim=2\n",
                ")\n",
                "print(rec.shape)\n",
                "print(rec.dtype.descr)\n",
                "\n",
                "print('transformation matrix')\n",
                "T = transformation.matrix(t=[-15, 10], s=[0.8, -0.8])\n",
                "print(T)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "grid = Grid(proj, rec, T)"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "Let's inspect the properties of the raster."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "print('shape:')\n",
                "print(grid.shape)\n",
                "print('number of cells:')\n",
                "print(grid.count)\n",
                "print('fields:')\n",
                "print(grid.dtype)\n",
                "print('projection:')\n",
                "print(grid.proj)\n",
                "print('transformation matrix:')\n",
                "print(np.round(grid.t, 2))\n",
                "print('origin:')\n",
                "print(np.round(grid.t.origin, 2).tolist())\n",
                "print('extent:')\n",
                "print(np.round(grid.extent(), 2).tolist())"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "Now, we visualize the x 'indices' of the raster."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(10, 10))\n",
                "plt.xlabel('X (pixels)')\n",
                "plt.ylabel('Y (pixels)')\n",
                "\n",
                "plt.imshow(grid.indices[:, :, 0])\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "You might have noticed, that the field 'coords' has been ",
                "implicitly added to the record array. The coordinates ",
                "correspond to the centers of the raster cells."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "print(np.round(grid.coords, 2))"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "Based on these coordinates we create an additional field ",
                "representing a surface."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "x = grid.coords[:, :, 0]\n",
                "y = grid.coords[:, :, 1]\n",
                "dist = np.sqrt(x ** 2 + y ** 2)\n",
                "z = 9 + 10 * (np.sin(0.5 * x) / np.sqrt(dist + 1) + np.cos(0.5 * y) / np.sqrt(dist + 1))\n",
                "grid = grid.add_fields([('z', float)], data=[z])\n",
                "print(grid.dtype.descr)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(10, 10))\n",
                "plt.xlabel('X (pixels)')\n",
                "plt.ylabel('Y (pixels)')\n",
                "\n",
                "plt.imshow(grid.z, cmap='coolwarm')\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "If we like to treat the raster as a point cloud or list of ",
                "points, we call the 'records' function. As a result we ",
                "receive a flattened version of the raster records."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "print('records type:')\n",
                "print(type(grid.records()))\n",
                "print('records shape:')\n",
                "print(grid.records().shape)\n",
                "print('coords:')\n",
                "print(np.round(grid.records().coords, 2))"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "We use these flattened coordinates to visualize the ",
                "centers of the raster cells."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(10, 10))\n",
                "ax = plt.axes(aspect='equal')\n",
                "plt.xlabel('X (m)')\n",
                "plt.ylabel('Y (m)')\n",
                "\n",
                "plt.scatter(*grid.records().coords.T, c=grid.records().z, cmap='coolwarm')\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "## GeoRecords\n",
                "The 'Grid' class presented before extends the 'GeoRecords' ",
                "class. Grid objects, like rasters or voxels are well ",
                "structured, while GeoRecords in general are just a ",
                "collection of points."
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "To understand the usage of the GeoRecords, we crate a ",
                "three dimensional point cloud using the coordinates ",
                "derived before. We also use the same coordinate reference ",
                "system. The creation of a GeoRecord array requires for a ",
                "record array with at least a field 'coords' specifying the ",
                "point coordinates."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "rec = nptools.recarray({\n",
                    "\t'coords': np.vstack([\n",
                        "\t\tgrid.records().coords[:, 0],\n",
                        "\t\tgrid.records().coords[:, 1],\n",
                        "\t\tgrid.records().z\n",
                    "\t]).T,\n",
                    "\t'z': grid.records().z\n",
                "})\n",
                "geoRecords = GeoRecords(grid.proj, rec)"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "We inspect the properties of the point cloud first."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "print('shape:')\n",
                "print(geoRecords.shape)\n",
                "print('number of points:')\n",
                "print(geoRecords.count)\n",
                "print('fields:')\n",
                "print(geoRecords.dtype)\n",
                "print('projection:')\n",
                "print(geoRecords.proj)\n",
                "print('transformation matrix:')\n",
                "print(np.round(geoRecords.t, 2))\n",
                "print('origin:')\n",
                "print(np.round(geoRecords.t.origin, 2).tolist())\n",
                "print('extent:')\n",
                "print(np.round(geoRecords.extent(), 2).tolist())"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "Before we visualize the point cloud, we define the axis ",
                "limits."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "axes_lims = Extent([\n",
                "\tgeoRecords.extent().center - 0.5 * geoRecords.extent().ranges.max(),\n",
                "\tgeoRecords.extent().center + 0.5 * geoRecords.extent().ranges.max()\n",
                "])"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(15, 15))\n",
                "ax = plt.axes(projection='3d')\n",
                "ax.set_xlim(axes_lims[0], axes_lims[3])\n",
                "ax.set_ylim(axes_lims[1], axes_lims[4])\n",
                "ax.set_zlim(axes_lims[2], axes_lims[5])\n",
                "ax.set_xlabel('X (m)')\n",
                "ax.set_ylabel('Y (m)')\n",
                "ax.set_zlabel('Z (m)')\n",
                "\n",
                "ax.scatter(*geoRecords.coords.T, c=geoRecords.z, cmap='coolwarm')\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "## Transformation\n",
                "For some applications we might like to transform the raster ",
                "coordinates a bit."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "T = transformation.matrix(t=[15, -10], s=[1.5, 2], r=10*np.pi/180, order='trs')\n",
                "tcoords = transformation.transform(grid.records().coords, T)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(10, 10))\n",
                "ax = plt.axes(aspect='equal')\n",
                "plt.xlabel('X (m)')\n",
                "plt.ylabel('Y (m)')\n",
                "\n",
                "plt.scatter(*tcoords.T, c=grid.records().z, cmap='coolwarm')\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "Or we roto-translate the raster."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "T = transformation.matrix(t=[1, 2], r=20*np.pi/180)\n",
                "grid.transform(T)\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(10, 10))\n",
                "ax = plt.axes(aspect='equal')\n",
                "plt.xlabel('X (m)')\n",
                "plt.ylabel('Y (m)')\n",
                "\n",
                "plt.scatter(*grid.records().coords.T, c=grid.records().z, cmap='coolwarm')\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "## IndexKD\n",
                "The 'GeoRecords' class provides a 'IndexKD' instance to ",
                "perform spatial neighborhood queries efficiently."
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "### Radial filtering\n",
                "We begin with filtering the points within a sphere around ",
                "some points. As a result, we receive a list of point ",
                "indices which can be used for sub-sampling."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "coords = [[-5, 0, 8], [10, -5, 5]]\n",
                "r = 6.0\n"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "Once in 3D ..."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fids_list = geoRecords.indexKD().ball(coords, r)\n",
                "print(len(fids_list))\n",
                "print(fids_list[0])\n",
                "print(fids_list[1])"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(15, 15))\n",
                "ax = plt.axes(projection='3d')\n",
                "ax.set_xlim(axes_lims[0], axes_lims[3])\n",
                "ax.set_ylim(axes_lims[1], axes_lims[4])\n",
                "ax.set_zlim(axes_lims[2], axes_lims[5])\n",
                "ax.set_xlabel('X (m)')\n",
                "ax.set_ylabel('Y (m)')\n",
                "ax.set_zlabel('Z (m)')\n",
                "\n",
                "ax.scatter(*geoRecords.coords.T, c=geoRecords.z, cmap='coolwarm', marker='.')\n",
                "for fids in fids_list:\n",
                "\tax.scatter(*geoRecords[fids].coords.T, s=100)\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "... and once in 2D."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fids_list = geoRecords.indexKD(2).ball(coords, r)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(15, 15))\n",
                "ax = plt.axes(projection='3d')\n",
                "ax.set_xlim(axes_lims[0], axes_lims[3])\n",
                "ax.set_ylim(axes_lims[1], axes_lims[4])\n",
                "ax.set_zlim(axes_lims[2], axes_lims[5])\n",
                "ax.set_xlabel('X (m)')\n",
                "ax.set_ylabel('Y (m)')\n",
                "ax.set_zlabel('Z (m)')\n",
                "\n",
                "ax.scatter(*geoRecords.coords.T, c=geoRecords.z, cmap='coolwarm', marker='.')\n",
                "for fids in fids_list:\n",
                "\tax.scatter(*geoRecords[fids].coords.T, s=100)\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "Of course, we can do the same with the raster."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fids_list = grid.indexKD().ball(coords, r)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(10, 10))\n",
                "ax = plt.axes(aspect='equal')\n",
                "plt.xlabel('X (m)')\n",
                "plt.ylabel('Y (m)')\n",
                "\n",
                "plt.scatter(*grid.records().coords.T, c=grid.records().z, cmap='coolwarm', marker='.')\n",
                "for fids in fids_list:\n",
                "\tax.scatter(*grid.records()[fids].coords.T)\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "### Nearest neighbor filtering\n",
                "We can also filter the nearest neighbors of the points ",
                "given before. Next to a list of point indices, we receive ",
                "a list of point distances of the same shape."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "k=50"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "Once in 3D ..."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "dists_list, fids_list = geoRecords.indexKD().knn(coords, k)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(15, 15))\n",
                "ax = plt.axes(projection='3d')\n",
                "ax.set_xlim(axes_lims[0], axes_lims[3])\n",
                "ax.set_ylim(axes_lims[1], axes_lims[4])\n",
                "ax.set_zlim(axes_lims[2], axes_lims[5])\n",
                "ax.set_xlabel('X (m)')\n",
                "ax.set_ylabel('Y (m)')\n",
                "ax.set_zlabel('Z (m)')\n",
                "\n",
                "ax.scatter(*geoRecords.coords.T, c=geoRecords.z, cmap='coolwarm', marker='.')\n",
                "for fids in fids_list:\n",
                "\tax.scatter(*geoRecords[fids].coords.T, s=100)\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "... and once in 2D."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "dists_list, fids_list = geoRecords.indexKD(2).knn(coords, k)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(15, 15))\n",
                "ax = plt.axes(projection='3d')\n",
                "ax.set_xlim(axes_lims[0], axes_lims[3])\n",
                "ax.set_ylim(axes_lims[1], axes_lims[4])\n",
                "ax.set_zlim(axes_lims[2], axes_lims[5])\n",
                "ax.set_xlabel('X (m)')\n",
                "ax.set_ylabel('Y (m)')\n",
                "ax.set_zlabel('Z (m)')\n",
                "\n",
                "ax.scatter(*geoRecords.coords.T, c=geoRecords.z, cmap='coolwarm', marker='.')\n",
                "for fids in fids_list:\n",
                "\tax.scatter(*geoRecords[fids].coords.T, s=100)\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "And again, once with the raster."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "dists_list, fids_list = grid.indexKD(2).knn(coords, k)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(10, 10))\n",
                "ax = plt.axes(aspect='equal')\n",
                "plt.xlabel('X (m)')\n",
                "plt.ylabel('Y (m)')\n",
                "\n",
                "ax.scatter(*grid.records().coords.T, c=geoRecords.z, cmap='coolwarm', marker='.')\n",
                "for fids in fids_list:\n",
                "\tax.scatter(*grid.records()[fids].coords.T)\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "### Point counting\n",
                "We have the option to count the number of points within a ",
                "given radius. For this purpose we select a subset of the ",
                "raster first. Then, using the point cloud, we count the ",
                "number of raster cells within the given radius."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "grid_subset = grid[15:25, 30:40]\n",
                "count = grid_subset.indexKD(2).ball_count(r, geoRecords.coords)\n",
                "print(count)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(15, 15))\n",
                "ax = plt.axes(projection='3d')\n",
                "ax.set_xlim(axes_lims[0], axes_lims[3])\n",
                "ax.set_ylim(axes_lims[1], axes_lims[4])\n",
                "ax.set_zlim(axes_lims[2], axes_lims[5])\n",
                "ax.set_xlabel('X (m)')\n",
                "ax.set_ylabel('Y (m)')\n",
                "ax.set_zlabel('Z (m)')\n",
                "\n",
                "ax.scatter(*geoRecords.coords.T, c=count, cmap='YlOrRd')\n",
                "ax.scatter(*grid_subset.records().coords.T, color='black')\n",
                "ax.scatter(*grid.records().coords.T, color='gray', marker='.')\n",

                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "You might have noticed, that this was a fist step of data ",
                "fusion, since we related the raster cells to the point ",
                "cloud. This can also be done for nearest neighbor or ",
                "similar spatial queries regardless of dimension."
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "## Point filters\n",
                "To create a subset of points, we typically use some kind ",
                "of point filters. We begin with a duplicate point filter. ",
                "To ease the integration of such filters into your own ",
                "algorithms, an iterator is returned instead of a list of ",
                "point indices."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fids = list(filters.ball(geoRecords.indexKD(), 2.5))\n",
                "print(fids)\n",
                "print(len(fids))"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(15, 15))\n",
                "ax = plt.axes(projection='3d')\n",
                "ax.set_xlim(axes_lims[0], axes_lims[3])\n",
                "ax.set_ylim(axes_lims[1], axes_lims[4])\n",
                "ax.set_zlim(axes_lims[2], axes_lims[5])\n",
                "ax.set_xlabel('X (m)')\n",
                "ax.set_ylabel('Y (m)')\n",
                "ax.set_zlabel('Z (m)')\n",
                "\n",
                "ax.scatter(*geoRecords.coords.T, c=geoRecords.z, cmap='coolwarm')\n",
                "ax.scatter(*geoRecords[fids].coords.T, color='red', s=100)\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "Sometimes we like to filter local maxima of an attribute ",
                "using a given radius ..."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fids = list(filters.extrema(geoRecords.indexKD(2), geoRecords.z, 1.5))\n",
                "print(fids)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(15, 15))\n",
                "ax = plt.axes(projection='3d')\n",
                "ax.set_xlim(axes_lims[0], axes_lims[3])\n",
                "ax.set_ylim(axes_lims[1], axes_lims[4])\n",
                "ax.set_zlim(axes_lims[2], axes_lims[5])\n",
                "ax.set_xlabel('X (m)')\n",
                "ax.set_ylabel('Y (m)')\n",
                "ax.set_zlabel('Z (m)')\n",
                "\n",
                "ax.scatter(*geoRecords.coords.T, c=geoRecords.z, cmap='coolwarm')\n",
                "ax.scatter(*geoRecords[fids].coords.T, color='red', s=100)\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "... or find local minima."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fids = list(filters.extrema(geoRecords.indexKD(2), geoRecords.z, 1.5, inverse=True))\n",
                "print(fids)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(15, 15))\n",
                "ax = plt.axes(projection='3d')\n",
                "ax.set_xlim(axes_lims[0], axes_lims[3])\n",
                "ax.set_ylim(axes_lims[1], axes_lims[4])\n",
                "ax.set_zlim(axes_lims[2], axes_lims[5])\n",
                "ax.set_xlabel('X (m)')\n",
                "ax.set_ylabel('Y (m)')\n",
                "ax.set_zlabel('Z (m)')\n",
                "\n",
                "ax.scatter(*geoRecords.coords.T, c=geoRecords.z, cmap='coolwarm')\n",
                "ax.scatter(*geoRecords[fids].coords.T, color='blue', s=100)\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "## Smoothing\n",
                "To compensate for noise, or receive just a smoother result ",
                "we can use smoothing algorithms. The algorithm presented ",
                "here averages the coordinates of the nearest neighbors."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "scoords = smoothing.mean_knn(geoRecords.coords, 20, num_iter=3)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(15, 15))\n",
                "ax = plt.axes(projection='3d')\n",
                "ax.set_xlim(axes_lims[0], axes_lims[3])\n",
                "ax.set_ylim(axes_lims[1], axes_lims[4])\n",
                "ax.set_zlim(axes_lims[2], axes_lims[5])\n",
                "ax.set_xlabel('X (m)')\n",
                "ax.set_ylabel('Y (m)')\n",
                "ax.set_zlabel('Z (m)')\n",
                "\n",
                "ax.scatter(*geoRecords.coords.T, c=geoRecords.z, cmap='coolwarm', marker='.')\n",
                "ax.plot_trisurf(*scoords.T, cmap='gist_earth')\n",
                "plt.show()"
            ]
        },
        {
            "cell_type":"markdown",
            "metadata":{},
            "source":[
                "## Clustering\n",
                "A common problem is to cluster point clouds. Here we use ",
                "a clustering algorithm, which assigns points iteratively ",
                "to the most dominant class within a given radius. By ",
                "iterating from top to bottom, the points are assigned to ",
                "the hills of the surface."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "order = np.argsort(geoRecords.z)[::-1]\n",
                "cluster_indices = clustering.majority_clusters(geoRecords.indexKD(), 5.0, order=order)\n",
                "print(cluster_indices)\n",
                "cluster_dict = classification.classes_to_dict(cluster_indices)\n",
                "print(cluster_dict)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig = plt.figure(figsize=(15, 15))\n",
                "ax = plt.axes(projection='3d')\n",
                "ax.set_xlim(axes_lims[0], axes_lims[3])\n",
                "ax.set_ylim(axes_lims[1], axes_lims[4])\n",
                "ax.set_zlim(axes_lims[2], axes_lims[5])\n",
                "ax.set_xlabel('X (m)')\n",
                "ax.set_ylabel('Y (m)')\n",
                "ax.set_zlabel('Z (m)')\n",
                "\n",
                "for fids in cluster_dict.values():\n",
                    "\tax.scatter(*geoRecords[fids].coords.T, s=100)\n",
                "plt.show()"
            ]
        }
    ],
    "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.6.4"
        }
    },
    "nbformat":4,
    "nbformat_minor":2
}