efface plot rgb-image
[staff/bgs2/sendbox.git] / siegfried_extract.py
CommitLineData
29878d29
SB
1# -*- coding: utf-8 -*-\r
2"""\r
3Created on Tue Oct 04 13:33:48 2016\r
4\r
5@author: bgs2\r
6"""\r
7\r
8import matplotlib.pyplot as plt\r
9import numpy as np\r
10import cv2\r
11from scipy import misc\r
12import os\r
13import gdal\r
14import scipy.ndimage.filters as spfilters\r
15import scipy.ndimage as ndimage\r
16import skimage.morphology as morf\r
17from skimage.morphology import square\r
18\r
19#sudo mount -t cifs //smb.ceph.bfh.ch/bfh/Geodata/ ~/mnt/geodaten/ -o user=bgs2,uid=1000\r
20\r
21#sudo mount -t cifs //smb.ceph.bfh.ch/hafl/AGR-Bodenkunde/ ~/mnt/AGR-Bodenkunde/ -o user=bgs2,uid=1000\r
22\r
23\r
24#%%\r
25setdw="/home/stephane/mnt/geodaten/Switzerland/national/Siegfriedkarte/GeoTIFF/TA25/"\r
26\r
27#datasource = "P:/LFE/HAFL/AGR-Bodenkunde/2016_carto-seeland/2016_cartographie des sols seeland/cartesiegfried_treiten.tif"\r
28datasource= setdw+"TA25_LV03-136-1915-1.tif"\r
29#ds = cv2.imread("U:/RaD/carto_moos/colorchart.JPG")\r
30gdal.AllRegister()\r
31ds=gdal.Open(datasource, 0)\r
32geotr=ds.GetGeoTransform()\r
33bandR=ds.GetRasterBand(1)\r
34bandG=ds.GetRasterBand(2)\r
35bandB=ds.GetRasterBand(3)\r
36dsR=bandR.ReadAsArray()\r
37dsG=bandG.ReadAsArray()\r
38dsB=bandB.ReadAsArray()\r
39cols=ds.RasterXSize\r
40rows=ds.RasterYSize\r
41\r
42#cols=ds.shape[0]\r
43#rows=ds.shape[1]\r
44 \r
45#%%\r
46imrgbr=np.dstack((dsR,dsG,dsB))\r
47plt.imshow(imrgbr)\r
48plt.colorbar()\r
49plt.show()\r
50#imrgbr=lumequal(imrgbr)\r
51linedeb=4100\r
52linefin=4300\r
53deb=2400\r
54fin=2600\r
55#plt.imshow(imrgbr[4000:6000,2000:4000])\r
56plt.imshow(imrgbr[linedeb:linefin,deb:fin])\r
57plt.show()\r
58#%% conversion into HSV after gaussian filter\r
59dsGau=cv2.medianBlur(imrgbr,3)\r
60dshsv=cv2.cvtColor(dsGau,cv2.COLOR_RGB2HSV)\r
61plt.imshow(dshsv[linedeb:linefin,deb:fin,1] )\r
62plt.colorbar()\r
63plt.show()\r
64\r
65#%%\r
66def fextract(filename, h1,h2, s1,s2, v1 ,v2, value):\r
67 low_brown = np.array([h1,s1,v1])\r
68 high_brown = np.array([h2,s2,v2])\r
69 mask=cv2.inRange(dshsv, low_brown, high_brown)/255*value\r
70 return mask\r
71#%%\r
72\r
73#%%#%%\r
74line=4200\r
75deb=2400\r
76fin=2600\r
77plt.plot(imrgbr[line,deb:fin,0], 'red')\r
78plt.plot(imrgbr[line,deb:fin,1], 'green')\r
79plt.plot(imrgbr[line,deb:fin,2], 'blue')\r
80plt.show()\r
81#%%\r
29878d29 82\r
29878d29
SB
83line=2250\r
84\r
85plt.plot(dshsv[linedeb:linefin,line,0], 'red')\r
86plt.plot(dshsv[linedeb:linefin,line,1], 'green')\r
87plt.plot(dshsv[linedeb:linefin,line,2], 'blue')\r
88plt.show()\r
89#%%\r
90#%%filration and removing for noise\r
91masktot=np.zeros((np.shape(dsR)[0],np.shape(dsR)[1]), dtype='float32')\r
92maskblue = fextract(dshsv[:,:,:], 70,150,60,230,120,230,1)\r
93maskblack = fextract(dshsv[:,:,:], 0,220,0,220,0,150,1)\r
94mask=np.ones((2,2), dtype='int')\r
95maskblackClose=morf.erosion(maskblack, mask)\r
96maskblueClose=morf.erosion(maskblue, mask)\r
97mask=np.ones((5,5), dtype='int')\r
98#maskblueClose2=morf.closing(maskblueClose, mask)\r
99maskblackClose2=morf.dilation(maskblackClose, mask)\r
100maskblueClose2=morf.dilation(maskblueClose, mask)\r
101plt.imshow(maskblueClose2[linedeb:linefin,deb:fin], 'Blues')\r
102plt.show()\r
103plt.colorbar()\r
104#%%\r
105#maskmini=np.array([[0, 1, 0], [ 1, 1, 1], [0, 1, 1]], dtype='int')\r
106maskmini=np.ones((1,1), dtype='int')\r
107maskbluesquel=morf.skeletonize(maskblueClose2)\r
108maskbluesqueldil=morf.dilation(maskbluesquel, maskmini)\r
109plt.imshow(maskbluesquel[linedeb:linefin,deb:fin], 'Blues')\r
110plt.show()\r
111#%%\r
112def fcanny(nomfile,n,m):\r
113 candata = cv2.Canny(nomfile,n,m)\r
114 return candata\r
115\r
116#%%\r
117datacanny=fcanny(np.uint8(masktotmed),3,3)\r
118datacannyClose=morf.closing(datacanny, square(2))\r
119plt.imshow(datacannyClose[3500:5500,3500:7000],'gray')#pl.colorbar()\r
120#%%\r
121driver=gdal.GetDriverByName('GTiff')\r
122bavoisdiffbinlab=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
123bavoisdiffbinlab.SetGeoTransform(geotr)\r
124bavoisdiffbinlab.GetRasterBand(1).WriteArray(maskbluesqueldil)\r
125bavoisdiffbinlab=None\r
126#%%\r
127driver=gdal.GetDriverByName('GTiff')\r
128bavoisdiffbinlab=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
129bavoisdiffbinlab.SetGeoTransform(geotr)\r
130bavoisdiffbinlab.GetRasterBand(1).WriteArray(maskblackClose2)\r
131bavoisdiffbinlab=None\r
132\r
133#%% labelisation des objets\r
134mask_lab = np.array([[0, 1, 0],\r
135 [1, 1, 1],\r
136 [0, 1, 0]])\r
137 \r
138labeled_array, num_features = ndimage.measurements.label(maskClose, structure=mask_lab)\r
139plt.imshow(labeled_array)\r
140plt.colorbar() \r
141\r
142#%%\r