Interpolate missing data
This commit is contained in:
parent
34aac07fdd
commit
f003715d34
@ -6,11 +6,16 @@
|
|||||||
"metadata": {},
|
"metadata": {},
|
||||||
"outputs": [],
|
"outputs": [],
|
||||||
"source": [
|
"source": [
|
||||||
|
"import sys\n",
|
||||||
"from pathlib import Path\n",
|
"from pathlib import Path\n",
|
||||||
"import numpy as np\n",
|
"import numpy as np\n",
|
||||||
"import matplotlib.pyplot as plt\n",
|
"import matplotlib.pyplot as plt\n",
|
||||||
"import laspy\n",
|
"import laspy\n",
|
||||||
"\n",
|
"\n",
|
||||||
|
"triskele_path = Path('../triskele/python/')\n",
|
||||||
|
"sys.path.append(str(triskele_path.resolve()))\n",
|
||||||
|
"import triskele\n",
|
||||||
|
"\n",
|
||||||
"figsize = (16, 9)"
|
"figsize = (16, 9)"
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
@ -64,12 +69,42 @@
|
|||||||
" print('file: {}, point count: {}'.format(file, infile.header.get_count()))\n",
|
" print('file: {}, point count: {}'.format(file, infile.header.get_count()))\n",
|
||||||
" xdata.append(infile.x)\n",
|
" xdata.append(infile.x)\n",
|
||||||
" ydata.append(infile.y)\n",
|
" ydata.append(infile.y)\n",
|
||||||
" idata.append(infile.intensity)\n",
|
" #idata.append(infile.intensity)\n",
|
||||||
|
" idata.append(infile.z)\n",
|
||||||
|
" \n",
|
||||||
" \n",
|
" \n",
|
||||||
"data = np.array((np.concatenate(xdata), np.concatenate(ydata), np.concatenate(idata))).T\n",
|
"data = np.array((np.concatenate(xdata), np.concatenate(ydata), np.concatenate(idata))).T\n",
|
||||||
"data.shape"
|
"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",
|
"cell_type": "markdown",
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
@ -83,11 +118,21 @@
|
|||||||
"metadata": {},
|
"metadata": {},
|
||||||
"outputs": [],
|
"outputs": [],
|
||||||
"source": [
|
"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",
|
"ground_size = np.array((data[:,0].max() - data[:,0].min(), data[:,1].max() - data[:,1].min()))\n",
|
||||||
"display(ground_size)\n",
|
"display(ground_size)\n",
|
||||||
"\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')"
|
"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",
|
"cell_type": "code",
|
||||||
"execution_count": null,
|
"execution_count": null,
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"outputs": [],
|
"outputs": [],
|
||||||
"source": [
|
"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()"
|
"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": {
|
"metadata": {
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user