premiere version analyse Sigfreid
authorStephane Burgos <stephane.burgos@bfh.ch>
Tue, 28 Nov 2017 13:53:33 +0000 (14:53 +0100)
committerStephane Burgos <stephane.burgos@bfh.ch>
Tue, 28 Nov 2017 13:53:33 +0000 (14:53 +0100)
siegfried_extract.py [new file with mode: 0644]

diff --git a/siegfried_extract.py b/siegfried_extract.py
new file mode 100644 (file)
index 0000000..0fe5b57
--- /dev/null
@@ -0,0 +1,149 @@
+# -*- coding: utf-8 -*-\r
+"""\r
+Created on Tue Oct 04 13:33:48 2016\r
+\r
+@author: bgs2\r
+"""\r
+\r
+import matplotlib.pyplot as plt\r
+import numpy as np\r
+import cv2\r
+from scipy import misc\r
+import os\r
+import gdal\r
+import scipy.ndimage.filters as spfilters\r
+import scipy.ndimage as ndimage\r
+import skimage.morphology as morf\r
+from skimage.morphology import square\r
+\r
+#sudo mount -t cifs //smb.ceph.bfh.ch/bfh/Geodata/ ~/mnt/geodaten/ -o user=bgs2,uid=1000\r
+\r
+#sudo mount -t cifs //smb.ceph.bfh.ch/hafl/AGR-Bodenkunde/ ~/mnt/AGR-Bodenkunde/ -o user=bgs2,uid=1000\r
+\r
+\r
+#%%\r
+setdw="/home/stephane/mnt/geodaten/Switzerland/national/Siegfriedkarte/GeoTIFF/TA25/"\r
+\r
+#datasource = "P:/LFE/HAFL/AGR-Bodenkunde/2016_carto-seeland/2016_cartographie des sols seeland/cartesiegfried_treiten.tif"\r
+datasource= setdw+"TA25_LV03-136-1915-1.tif"\r
+#ds = cv2.imread("U:/RaD/carto_moos/colorchart.JPG")\r
+gdal.AllRegister()\r
+ds=gdal.Open(datasource, 0)\r
+geotr=ds.GetGeoTransform()\r
+bandR=ds.GetRasterBand(1)\r
+bandG=ds.GetRasterBand(2)\r
+bandB=ds.GetRasterBand(3)\r
+dsR=bandR.ReadAsArray()\r
+dsG=bandG.ReadAsArray()\r
+dsB=bandB.ReadAsArray()\r
+cols=ds.RasterXSize\r
+rows=ds.RasterYSize\r
+\r
+#cols=ds.shape[0]\r
+#rows=ds.shape[1]\r
+    \r
+#%%\r
+imrgbr=np.dstack((dsR,dsG,dsB))\r
+plt.imshow(imrgbr)\r
+plt.colorbar()\r
+plt.show()\r
+#imrgbr=lumequal(imrgbr)\r
+linedeb=4100\r
+linefin=4300\r
+deb=2400\r
+fin=2600\r
+#plt.imshow(imrgbr[4000:6000,2000:4000])\r
+plt.imshow(imrgbr[linedeb:linefin,deb:fin])\r
+plt.show()\r
+#%% conversion into HSV after gaussian filter\r
+dsGau=cv2.medianBlur(imrgbr,3)\r
+dshsv=cv2.cvtColor(dsGau,cv2.COLOR_RGB2HSV)\r
+plt.imshow(dshsv[linedeb:linefin,deb:fin,1] )\r
+plt.colorbar()\r
+plt.show()\r
+\r
+#%%\r
+def fextract(filename, h1,h2, s1,s2, v1 ,v2, value):\r
+    low_brown = np.array([h1,s1,v1])\r
+    high_brown = np.array([h2,s2,v2])\r
+    mask=cv2.inRange(dshsv, low_brown, high_brown)/255*value\r
+    return mask\r
+#%%\r
+\r
+#%%#%%\r
+line=4200\r
+deb=2400\r
+fin=2600\r
+plt.plot(imrgbr[line,deb:fin,0], 'red')\r
+plt.plot(imrgbr[line,deb:fin,1], 'green')\r
+plt.plot(imrgbr[line,deb:fin,2], 'blue')\r
+plt.show()\r
+#%%\r
+line=4200\r
+\r
+plt.plot(dshsv[line,deb:fin,0], 'red')\r
+plt.plot(dshsv[line,deb:fin,1], 'green')\r
+plt.plot(dshsv[line,deb:fin,2], 'blue')\r
+plt.colorbar()\r
+plt.show()\r
+#%%\r
+line=2250\r
+\r
+plt.plot(dshsv[linedeb:linefin,line,0], 'red')\r
+plt.plot(dshsv[linedeb:linefin,line,1], 'green')\r
+plt.plot(dshsv[linedeb:linefin,line,2], 'blue')\r
+plt.show()\r
+#%%\r
+#%%filration and removing for noise\r
+masktot=np.zeros((np.shape(dsR)[0],np.shape(dsR)[1]), dtype='float32')\r
+maskblue = fextract(dshsv[:,:,:], 70,150,60,230,120,230,1)\r
+maskblack = fextract(dshsv[:,:,:], 0,220,0,220,0,150,1)\r
+mask=np.ones((2,2), dtype='int')\r
+maskblackClose=morf.erosion(maskblack, mask)\r
+maskblueClose=morf.erosion(maskblue, mask)\r
+mask=np.ones((5,5), dtype='int')\r
+#maskblueClose2=morf.closing(maskblueClose, mask)\r
+maskblackClose2=morf.dilation(maskblackClose, mask)\r
+maskblueClose2=morf.dilation(maskblueClose, mask)\r
+plt.imshow(maskblueClose2[linedeb:linefin,deb:fin], 'Blues')\r
+plt.show()\r
+plt.colorbar()\r
+#%%\r
+#maskmini=np.array([[0,  1, 0], [ 1,  1,  1],  [0,  1, 1]], dtype='int')\r
+maskmini=np.ones((1,1), dtype='int')\r
+maskbluesquel=morf.skeletonize(maskblueClose2)\r
+maskbluesqueldil=morf.dilation(maskbluesquel, maskmini)\r
+plt.imshow(maskbluesquel[linedeb:linefin,deb:fin], 'Blues')\r
+plt.show()\r
+#%%\r
+def fcanny(nomfile,n,m):\r
+    candata = cv2.Canny(nomfile,n,m)\r
+    return candata\r
+\r
+#%%\r
+datacanny=fcanny(np.uint8(masktotmed),3,3)\r
+datacannyClose=morf.closing(datacanny, square(2))\r
+plt.imshow(datacannyClose[3500:5500,3500:7000],'gray')#pl.colorbar()\r
+#%%\r
+driver=gdal.GetDriverByName('GTiff')\r
+bavoisdiffbinlab=driver.Create("/home/stephane/mnt/AGR-Bodenkunde/2016_carto-seeland/analyse/TA25_LV03-136-1915-1-blue_squel2.tif", cols,rows,1,gdal.GDT_Float32)\r
+bavoisdiffbinlab.SetGeoTransform(geotr)\r
+bavoisdiffbinlab.GetRasterBand(1).WriteArray(maskbluesqueldil)\r
+bavoisdiffbinlab=None\r
+#%%\r
+driver=gdal.GetDriverByName('GTiff')\r
+bavoisdiffbinlab=driver.Create("/home/stephane/mnt/AGR-Bodenkunde/2016_carto-seeland/analyse/TA25_LV03-136-1915-1-black.tif", cols,rows,1,gdal.GDT_Float32)\r
+bavoisdiffbinlab.SetGeoTransform(geotr)\r
+bavoisdiffbinlab.GetRasterBand(1).WriteArray(maskblackClose2)\r
+bavoisdiffbinlab=None\r
+\r
+#%% labelisation des objets\r
+mask_lab = np.array([[0, 1, 0],\r
+                     [1, 1, 1],\r
+                     [0, 1, 0]])\r
+                     \r
+labeled_array, num_features = ndimage.measurements.label(maskClose, structure=mask_lab)\r
+plt.imshow(labeled_array)\r
+plt.colorbar()  \r
+\r
+#%%\r