Update Notebooks

This commit is contained in:
Florent Guiotte 2018-03-22 07:43:35 +01:00
parent 21f13dc864
commit ae823b3a26
3 changed files with 570 additions and 25 deletions

View File

@ -4,7 +4,7 @@
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": [ "source": [
"# GDAL vs Matplotlib\n", "# GDAL vs Matplotlib (vs OpenCV vs LibTIFF)\n",
"\n", "\n",
"We have in data fusion contest dataset TIFF file with 32 bit per sample. For now TRISKELE only work with less than 17 bps: `BOOST_ASSERT (bits < 17);`. I want to ensure that we can oppen such data with Python.\n" "We have in data fusion contest dataset TIFF file with 32 bit per sample. For now TRISKELE only work with less than 17 bps: `BOOST_ASSERT (bits < 17);`. I want to ensure that we can oppen such data with Python.\n"
] ]
@ -18,7 +18,10 @@
"import matplotlib.pyplot as plt\n", "import matplotlib.pyplot as plt\n",
"import gdal\n", "import gdal\n",
"from pathlib import Path\n", "from pathlib import Path\n",
"import subprocess" "import subprocess\n",
"import numpy as np\n",
"import pandas as pd\n",
"import cv2"
] ]
}, },
{ {
@ -27,7 +30,8 @@
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"file = Path('../Data/phase1_rasters/DSM_C12/UH17c_GEF051_TR.tif')" "file = Path('../Data/phase1_rasters/DSM_C12/UH17c_GEF051_TR.tif')\n",
"ofile = Path('../Res/python_out.tif')"
] ]
}, },
{ {
@ -86,7 +90,6 @@
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"import numpy as np\n",
"(np.power(2, (np.arange(4)[::-1] * 8)) * mat_data).sum(axis=2).shape" "(np.power(2, (np.arange(4)[::-1] * 8)) * mat_data).sum(axis=2).shape"
] ]
}, },
@ -127,13 +130,341 @@
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"plt.figure(figsize=(32,9))\n", "raster.max()"
"plt.imshow(raster * (raster < .25 * 1e38))\n", ]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from ipywidgets import interact, interactive, fixed, interact_manual\n",
"import ipywidgets as widgets\n",
"\n",
"def disp(x):\n",
" display(x)\n",
" plt.figure(figsize=(32,9))\n",
" plt.imshow(raster * (raster < x))\n",
" plt.colorbar()\n",
" plt.show()\n",
" plt.hist(raster.reshape(-1) * (raster.reshape(-1) < x), 100, log=True)\n",
" plt.show()\n",
" \n",
"interact(disp, x=widgets.IntSlider(min=raster.min(),max=80,value=80));"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raster_df = pd.DataFrame(raster.reshape(-1))\n",
"raster_df.describe()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rasterf = raster.copy()\n",
"rasterf[rasterf == rasterf.max()] = 0"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rasterf_df = pd.DataFrame(rasterf.reshape(-1))\n",
"rasterf_df.describe()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.imshow(rasterf)\n",
"plt.colorbar()\n", "plt.colorbar()\n",
"plt.show()\n",
"plt.hist(raster.reshape(-1) * (raster.reshape(-1) < .25 * 1e38), 100, log=True)\n",
"plt.show()" "plt.show()"
] ]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Write TIFF file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### OpenCV"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raster_cv2 = cv2.imread(str(file), -1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(raster_cv2 != raster).sum()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"str(ofile)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ls ../Res"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rasterf_im = np.expand_dims(rasterf, 2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cv2.imwrite(str(ofile), rasterf_im)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### GDAL"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"gdl_data.RasterXSize"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cols = gdl_data.RasterXSize\n",
"rows = gdl_data.RasterYSize\n",
"\n",
"format = 'GTiff'\n",
"driver = gdal.GetDriverByName(format)\n",
"dst_ds = driver.Create(str(ofile), cols, rows, 1, gdal.GDT_Byte)\n",
"dst_ds.GetRasterBand(1).WriteArray(rasterf, 0, 0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dst_ds = None"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Test "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"!tiffinfo ../Res/python_out.tif"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.imshow(rasterf)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"test_ofile = cv2.imread(str(ofile))\n",
"display(test_ofile.shape)\n",
"todf = pd.DataFrame(test_ofile.reshape(-1,3))\n",
"display(todf.describe())\n",
"plt.imshow(test_ofile[:,:,0])\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## LibTiff\n",
"\n",
"### Write"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import libtiff\n",
"\n",
"tiff = libtiff.TIFFimage(rasterf.astype(np.float16), description='test')\n",
"tiff.write_file(ofile)#, compression='lzw')\n",
"del tiff"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"!tiffinfo ../Res/python_out.tif"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ogdl_data = gdal.Open(str(ofile))\n",
"ogdl_data.GetMetadata(), ogdl_data.RasterCount"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"otest = ogdl_data.GetRasterBand(1).ReadAsArray()\n",
"otest.shape"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.imshow(otest)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(otest != rasterf).sum()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Précision au millimètre :"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"np.isclose(otest, rasterf, rtol=1e-3).sum() / rasterf.size"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Read"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tiff_r = libtiff.TIFF.open(file, mode='r')\n",
"tiff_r.read_image().shape"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.imshow(tiff_r.read_image())\n",
"plt.colorbar()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(tiff_r.read_image() != raster).sum()"
]
} }
], ],
"metadata": { "metadata": {

View File

@ -0,0 +1,202 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"from pathlib import Path\n",
"import numpy as np\n",
"import libtiff"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## `float32` to `uint16`"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Load raster"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"infile_path = Path('../Data/phase1_rasters/DSM_C12/UH17c_GEF051_TR.tif')\n",
"outfile_path = Path('../Res/dem_u16.tif')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"infile = libtiff.TIFF.open(infile_path)\n",
"raster = infile.read_image()\n",
"raster.shape"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.imshow(raster)\n",
"plt.colorbar()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Remove extremum"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raster.max(), raster[raster != raster.max()].max()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raster[raster == raster.max()] = raster[raster != raster.max()].max()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.imshow(raster)\n",
"plt.colorbar()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Scale to new representation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rep = np.iinfo(np.uint16)\n",
"\n",
"raster -= raster.min() - rep.min\n",
"raster *= (rep.max - rep.min) / (raster.max() - raster.min())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.imshow(raster)\n",
"plt.colorbar()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"np.iinfo(np.uint16).max, raster.max(), raster.astype(np.uint16).max()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Cast to new type"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raster = raster.astype(np.uint16)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.imshow(raster)\n",
"plt.colorbar()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Save raster"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"outfile = libtiff.TIFFimage(raster, description='TRISKELE raster')\n",
"outfile.write_file(outfile_path)#, compression='lzw')\n",
"outfile = None"
]
}
],
"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"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

View File

@ -22,7 +22,8 @@
"import subprocess\n", "import subprocess\n",
"import numpy as np\n", "import numpy as np\n",
"import matplotlib.pyplot as plt\n", "import matplotlib.pyplot as plt\n",
"import pandas as pd" "import pandas as pd\n",
"import libtiff"
] ]
}, },
{ {
@ -55,7 +56,8 @@
"outputs": [], "outputs": [],
"source": [ "source": [
"# Files\n", "# Files\n",
"infile = Path('../Data/phase1_rasters/DEM_C123_TLI/UH17_GEG05_TR.tif')\n", "#infile = Path('../Data/phase1_rasters/DEM_C123_TLI/UH17_GEG05_TR.tif')\n",
"infile = Path('../Res/dem_u16.tif')\n",
"outfile = Path('../Res/out.tiff')\n", "outfile = Path('../Res/out.tiff')\n",
"\n", "\n",
"area_conf = Path('../Res/area.txt')\n", "area_conf = Path('../Res/area.txt')\n",
@ -63,7 +65,7 @@
"inertia_conf = Path('../Res/inertia.txt')\n", "inertia_conf = Path('../Res/inertia.txt')\n",
"\n", "\n",
"# Attributes definition\n", "# Attributes definition\n",
"area = np.array([10,100,1000])\n", "area = np.arange(10) * 100 #* 10000#np.array([10,100,1000])\n",
"deviation = np.array([10,100,1000])\n", "deviation = np.array([10,100,1000])\n",
"inertia = np.array([10,100,1000])\n", "inertia = np.array([10,100,1000])\n",
"np.savetxt(area_conf, area, fmt='%d')\n", "np.savetxt(area_conf, area, fmt='%d')\n",
@ -72,10 +74,10 @@
"\n", "\n",
"process = [triskele_bin.absolute(),\n", "process = [triskele_bin.absolute(),\n",
" '-i', infile.absolute(),\n", " '-i', infile.absolute(),\n",
" '--tos-tree',\n", " '--max-tree',\n",
" '--area', area_conf.absolute(),\n", " '--area', area_conf.absolute(),\n",
" '--standard-deviation', deviation_conf.absolute(),\n", " #'--standard-deviation', deviation_conf.absolute(),\n",
" '--moment-of-inertia', inertia_conf.absolute(),\n", " #'--moment-of-inertia', inertia_conf.absolute(),\n",
" '-o', outfile.absolute()]\n", " '-o', outfile.absolute()]\n",
"\n", "\n",
"display(' '.join(str(v) for v in process))\n", "display(' '.join(str(v) for v in process))\n",
@ -90,8 +92,15 @@
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"A = ' '\n", "process = libtiff.TIFF.open(outfile)\n",
"A.join(str(v) for v in process)" "raster = process.read_image()\n",
"raster.shape\n",
"\n",
"for chan in range(raster.shape[2]):\n",
" plt.figure(figsize=(16,3))\n",
" plt.imshow(raster[:,:,chan])\n",
" plt.colorbar()\n",
" plt.show()"
] ]
}, },
{ {
@ -100,19 +109,22 @@
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"out = plt.imread(outfile.absolute())\n", "rasterf = raster.astype(np.float)\n",
"display(out.shape)\n", "plt.figure(figsize=(16,3))\n",
"\n", "plt.imshow(rasterf[:,:,0] - rasterf[:,:,-1])\n",
"dwidth = 32\n",
"plt.figure(figsize=(dwidth, dwidth * out.shape[0] / out.shape[1] * .80))\n",
"#plt.figure(figsize=(16,9))\n",
"\n",
"display_channel = 2\n",
"plt.imshow(out[:,:,display_channel] * (out[:,:,display_channel] < 50))\n",
"plt.colorbar()\n", "plt.colorbar()\n",
"plt.show()" "plt.show()"
] ]
}, },
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"raster[:,:,0] - raster[:,:,9]"
]
},
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": null,