From f003715d348f77f9d1adc462eb6e7f80d609621f Mon Sep 17 00:00:00 2001 From: Karamaz0V1 Date: Wed, 13 Jun 2018 07:10:40 +0200 Subject: [PATCH] Interpolate missing data --- .../Custom rasters from LiDAR data.ipynb | 314 +++++++++++++++++- 1 file changed, 310 insertions(+), 4 deletions(-) diff --git a/Notebooks/Custom rasters from LiDAR data.ipynb b/Notebooks/Custom rasters from LiDAR data.ipynb index 80104e7..41d8248 100644 --- a/Notebooks/Custom rasters from LiDAR data.ipynb +++ b/Notebooks/Custom rasters from LiDAR data.ipynb @@ -6,11 +6,16 @@ "metadata": {}, "outputs": [], "source": [ + "import sys\n", "from pathlib import Path\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import laspy\n", "\n", + "triskele_path = Path('../triskele/python/')\n", + "sys.path.append(str(triskele_path.resolve()))\n", + "import triskele\n", + "\n", "figsize = (16, 9)" ] }, @@ -64,12 +69,42 @@ " print('file: {}, point count: {}'.format(file, infile.header.get_count()))\n", " xdata.append(infile.x)\n", " ydata.append(infile.y)\n", - " idata.append(infile.intensity)\n", + " #idata.append(infile.intensity)\n", + " idata.append(infile.z)\n", + " \n", " \n", "data = np.array((np.concatenate(xdata), np.concatenate(ydata), np.concatenate(idata))).T\n", "data.shape" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.hist(data[:,2], 1000)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "extremum = .001\n", + "tresholds = np.percentile(data[:,2], [extremum, 1 - extremum])\n", + "display(tresholds)\n", + "\n", + "filtered_data = data.copy()\n", + "filtered_data[:,2][data[:,2] < tresholds[0]] = tresholds[0]\n", + "filtered_data[:,2][data[:,2] > tresholds[1]] = tresholds[1]\n", + "\n", + "plt.hist(filtered_data[:,2], 1000)\n", + "plt.show()" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -83,11 +118,21 @@ "metadata": {}, "outputs": [], "source": [ - "resolution = 2\n", + "resolution = 0.5\n", "ground_size = np.array((data[:,0].max() - data[:,0].min(), data[:,1].max() - data[:,1].min()))\n", "display(ground_size)\n", "\n", - "bins = ground_size / resolution" + "bins = np.rint(ground_size / resolution).astype(int)\n", + "bins" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "data[:,0].max(), data[:,0].min(), data[:,1].max(), data[:,1].min()" ] }, { @@ -178,15 +223,276 @@ "plt.imsave('../Res/hist_int.png', idisp, origin='lower')" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Create rasters DFC comptabible\n", + "\n", + "[Official description](http://www.grss-ieee.org/community/technical-committees/data-fusion/data-fusion-contest/) announce .5 m resolution on the DSMs.\n", + "\n", + "## Compare with existing DSMs\n", + "\n", + "### Shape" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "plt.scatter(infile.x, infile.y)\n", + "dfc_dsm = triskele.read('../Data/phase1_rasters/Intensity_C1/UH17_GI1F051_TR.tif')\n", + "dfc_dsm.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "idisp.shape" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Offset\n", + "\n", + "#### Filter DFC raster" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nodata_filter = dfc_dsm > 1e4\n", + "dfc_dsm[nodata_filter] = dfc_dsm[nodata_filter == False].max()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.hist(dfc_dsm.reshape(-1), 1000, label='DFC')\n", + "plt.hist(idisp[f], 1000, label='our', alpha=.8)\n", + "\n", + "plt.legend()\n", "plt.show()" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Visual test" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "offset_disp = np.moveaxis(np.array((np.flip(idisp, 0), dfc_dsm, np.zeros_like(idisp))), 0, 2)\n", + "display(offset_disp.shape)\n", + "\n", + "plt.figure(figsize=figsize)\n", + "plt.imshow(offset_disp / np.nanmax(offset_disp), origin='upper', extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Seems good." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.imsave('../Res/offset.png',offset_disp / np.nanmax(offset_disp), origin='upper')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Better raster with interpolation of missing data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.interpolate import griddata" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Test Griddata" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "gridx, gridy = np.mgrid[0:10, 0:10] + 10\n", + "gridx, gridy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sampled_data = np.array([[0, 0, 5], [9, 9, 5], [0, 9, 10], [9, 0, 10]])\n", + "sampled_data[:, :2] += 10\n", + "display(sampled_data)\n", + "\n", + "inter_data = griddata(sampled_data[:,:2], sampled_data[:,2], (gridx, gridy), method='cubic') # linear, nearest, cubic\n", + "\n", + "plt.imshow(inter_data, origin='lower')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Automatic grid computation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "test_resolution = .01\n", + "\n", + "sampled_data = np.array([[0, 0, 5], [9, 9, 5], [0, 9, 10], [20, 0, 10]])\n", + "\n", + "def rasterize(coords, values, resolution, method='cubic'):\n", + " xmin, xmax = coords[:,0].min(), coords[:,0].max()\n", + " ymin, ymax = coords[:,1].min(), coords[:,1].max()\n", + "\n", + " gridx, gridy = np.mgrid[xmin:xmax:resolution, ymin:ymax:resolution]\n", + " \n", + " return griddata(coords, values, (gridx, gridy), method=method) # linear, nearest, cubic\n", + "\n", + "plt.imshow(rasterize(sampled_data[:,:2], sampled_data[:,2], test_resolution, 'nearest').T, origin='lower')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.hist(data[:,2], 1000)\n", + "plt.show()\n", + "\n", + "should_have_keep_idata = data[:,2].copy()\n", + "should_have_keep_idata[data[:,2] > 1e4] = 1e4\n", + "\n", + "plt.hist(should_have_keep_idata, 700)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "treshold = -1\n", + "raster = rasterize(data[:treshold,:2], should_have_keep_idata[:treshold], 0.5, 'nearest')\n", + "\n", + "plt.figure(figsize=figsize)\n", + "plt.imshow(raster.T, origin='lower')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.imsave('../Res/my_raster.png', raster.T, origin='lower')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "idisp_complete = griddata(data[:,:2], data[:,2], test)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "X, Y = np.meshgrid(*([np.linspace(-1,1,200)] * 2))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "test = np.mgrid[0:bins[1], 0:bins[0]]\n", + "test" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "data[:,:2]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Save TIFF" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "triskele.write('../Res/my_intensity.tiff', np.flip(idisp, 0))" + ] } ], "metadata": {