-
Notifications
You must be signed in to change notification settings - Fork 4
/
analise_histograma_imagens.py
175 lines (118 loc) · 5.69 KB
/
analise_histograma_imagens.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
# -*- coding: utf-8 -*-
"""
Created on Tue Dec 27 11:59:07 2016
Programa para analisar o histograma de diversas imagens
Plotar o hist acumulado de várias imagens, a média, mediana etc
@author: m330625
"""
import glob
from osgeo import gdal
from osgeo import ogr
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
#pasta = 'I:\parana\8bit*.tif'
pasta = 'I:\\SUL\\epsg4326_21s_udm\\*.tif'
articulacao = 'C:\\geodb\\articulacao\\RapidEye_Tiles\\RapidEye_Tiles_Brasil.shp'
imagens = glob.glob(pasta)
def cospeHist(imagem):
"""Cospe o histograma da imagem, para a banda específica
banda assume valores de 1 a n... não começa no 0 (zero)"""
arq = gdal.Open(imagem)
hist_r = arq.GetRasterBand(1).GetHistogram()
hist_g = arq.GetRasterBand(2).GetHistogram()
hist_b = arq.GetRasterBand(3).GetHistogram()
arq = None
return hist_r, hist_g, hist_b
hist_rgb = [cospeHist(i) for i in imagens]
hist_r = [i[0] for i in hist_rgb if sum(i[0]) > 10] # if é para tirar histogramas com valor 0 para tudo
hist_g = [i[1] for i in hist_rgb if sum(i[1]) > 10]
hist_b = [i[2] for i in hist_rgb if sum(i[2]) > 10]
cdf_r = np.array([np.cumsum(i)/float(sum(i)) for i in hist_r])
cdf_g = np.array([np.cumsum(i)/float(sum(i)) for i in hist_g])
cdf_b = np.array([np.cumsum(i)/float(sum(i)) for i in hist_b])
def geraGraficos():
""" funçao para gerar os gráficos e analisar a bagaça"""
plt.subplot(221)
plt.plot(np.transpose(cdf_r), color='gray', ls='dashed')
plt.plot(np.mean(cdf_r, axis=0), color='red')
plt.plot(np.median(cdf_r, axis=0), color='blue')
plt.title('CDF for red band with mean and median values')
plt.subplot(222)
plt.plot(np.transpose(cdf_g), color='gray', ls='dashed')
plt.plot(np.mean(cdf_r, axis=0), color='red')
plt.plot(np.median(cdf_r, axis=0), color='blue')
plt.title('CDF for green band with mean and median values')
plt.subplot(223)
plt.plot(np.transpose(cdf_b), color='gray', ls='dashed')
plt.plot(np.mean(cdf_r, axis=0), color='red')
plt.plot(np.median(cdf_r, axis=0), color='blue')
plt.title('CDF for blue band with mean and median values')
return
# gerando função de transformação para cada banda
# usando valores medianos
# inv_cdf_cor relaciona percentual de ocorrência ao valor de 0 a 255
# depois aplica a função inversa na cdf de cada imagem para obter uma
# função que transforma os valores da imagem aos correspondentes na mediana
median_cdf_r = np.median(cdf_r, axis=0)
median_cdf_g = np.median(cdf_g, axis=0)
median_cdf_b = np.median(cdf_b, axis=0)
inv_cdf_r = interp1d(median_cdf_r, range(len(median_cdf_r)), bounds_error=0)
inv_cdf_g = interp1d(median_cdf_g, range(len(median_cdf_g)), bounds_error=0)
inv_cdf_b = interp1d(median_cdf_b, range(len(median_cdf_b)), bounds_error=0)
def executaBalanco():
for i in range(len(imagens)):
arq = gdal.Open(imagens[i])
saida = "I:\\SUL\\epsg_4326_21s_udm_balanco\\" + str(i) + '.tif'
transf_r = inv_cdf_r(cdf_r[i])
transf_g = inv_cdf_g(cdf_r[i])
transf_b = inv_cdf_b(cdf_r[i])
np.insert(transf_r, 0, 0)
np.insert(transf_g, 0, 0)
np.insert(transf_b, 0, 0)
transf_func_r = lambda x:transf_r[x]
transf_func_g = lambda x:transf_g[x]
transf_func_b = lambda x:transf_b[x]
file_type = 'GTiff'
driver = gdal.GetDriverByName(file_type)
outfile = driver.CreateCopy(saida, arq, 0, ['COMPRESS=LZW'])
corr_r = transf_func_r(arq.GetRasterBand(1).ReadAsArray())
corr_g = transf_func_g(arq.GetRasterBand(2).ReadAsArray())
corr_b = transf_func_b(arq.GetRasterBand(3).ReadAsArray())
arq = None
outfile.GetRasterBand(1).WriteArray(corr_r)
outfile.GetRasterBand(2).WriteArray(corr_g)
outfile.GetRasterBand(3).WriteArray(corr_b)
outfile = None
return
# esse procedimento não deu muito certo. Não melhorou em relação ao anterior
# imagens apresentam muitas diferenças entre datas.
# olhar de perto 4 cenas que apresentam distinção, para ver histograma melhor
cenas_str = ['2227621', '2227622', '2227521', '2227522']
cenas = []
for i in cenas_str:
cenas.append([j for j in enumerate(imagens) if i in j[1]])
def dataCenas(imagens, articulacao, saida):
"""Extrai a data das cenas e adiciona no shape de articulacao"""
entrada = ogr.Open(articulacao)
driver = ogr.GetDriverByName("ESRI Shapefile")
outfile = driver.CopyDataSource(entrada, saida)
entrada = None
layer = outfile.GetLayer()
layer.CreateField(ogr.FieldDefn('Data', ogr.OFTDate))
for i in imagens:
data_i = i[41:51].split('-')
cena = i[33:40]
layer.SetAttributeFilter('TILE_ID = %i' % int(cena))
for feature in layer:
feature.SetField('Data', int(data_i[0]), int(data_i[1]), int(data_i[2]), 0, 0, 0, 0)
layer.SetFeature(feature)
outfile= None
return
dataCenas(imagens, articulacao, 'I:\\articulacao_rapideye_2014.shp')
# olhando a data das cenas e a imagem com o balanço feito pelo 2SD
# aparentemente as diferenças não são só causadas pela data
# o porcentual de nuvens também pode atrapalhar a estimativa dos valores de SD
# o que bagunça o ajuste de cores
# devo implementar uma função para ignorar pixels com valores muito elevados (nuves) ou baixos (bordas)
# no cálculo do SD