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