-
Notifications
You must be signed in to change notification settings - Fork 25
/
ApplyMNFcoefficients.py
executable file
·313 lines (290 loc) · 12.4 KB
/
ApplyMNFcoefficients.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
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
#! /usr/bin/env python
########################################################################################################################
#
# ApplyMNFcoefficients.py
# A python script to perform MNF transformation to one image and then apply the coefficients to other images.
#
# WARNING!: this assume that all images are in the same format
#
# Info: The script perform MNF transformation to all raster images stored in a folder.
#
# Author: Javier Lopatin
# Email: javierlopatin@gmail.com
# Date: 09/08/2016
# Version: 1.0
#
# Usage:
#
# python MNF.py -i <Input raster from which copy the MNF coefficients>
# -c <Number of components>
# -m <Method option [default = 1]>
# -p <Preprocessing: Brightness Normalization of Hyperspectral data [Optional]>
# -s <Apply Savitzky Golay filtering [Optional]>
# -v <Accumulated explained variance [Optional]>
#
# --inputImage [-i]: input raster from which copy the MNF coefficients
#
# -- method [-m]: Method options: 1 (default) regular MNF transformation.
# 2 MNF inverse transformation.
#
# --preprop [-p]: Brightness Normalization presented in Feilhauer et al., 2010
#
# --SavitzkyGolay [-s]: Apply Savitzky Golay filtering
#
# --variance [-v]: Get the accumulative explained variance of MNF components
#
# examples:
# # Get the accumulated explained variance
# python ApplyMNFcoefficients.py -i image.tif -c 1 -v
#
# # with Brightness Normalization
# python ApplyMNFcoefficients.py -i image.tif-c 1 -v -p
#
# # Get the regular MNF transformation
# python ApplyMNFcoefficients.py -i image.tif -c 10
# python ApplyMNFcoefficients.py -i image.tif -c 10 -s # with Savitzky Golay
#
# # with Brightness Normalization
# python ApplyMNFcoefficients.py -i image.tif -c 10 -p
#
# # Get the reduced nose MNF with inverse transformation
# python ApplyMNFcoefficients.py -i image.tif-c 10 -m 2
# python ApplyMNFcoefficients.py -i image.tif-c 10 -m 2 -s # with Savitzky Golay
#
# # with Brightness Normalization
# python ApplyMNFcoefficients.py -i image.tif-c 10 -m 2 -p
#
#
# Bibliography:
#
# Feilhauer, H., Asner, G. P., Martin, R. E., Schmidtlein, S. (2010): Brightness-normalized Partial Least Squares
# Regression for hyperspectral data. Journal of Quantitative Spectroscopy and Radiative Transfer 111(12-13),
# pp. 1947–1957. 10.1016/j.jqsrt.2010.03.007
#
# C-I Change and Q Du. 1999. Interference and Noise-Adjusted Principal Components Analysis.
# IEEE TGRS, Vol 36, No 5.
#
########################################################################################################################
import os, glob, argparse
import numpy as np
from sklearn.decomposition import PCA
try:
import rasterio
except ImportError:
print("ERROR: Could not import Rasterio Python library.")
print("Check if Rasterio is installed.")
try:
import pysptools.noise as ns
except ImportError:
print("ERROR: Could not import Pysptools Python library.")
print("Check if Pysptools is installed.")
## Functions
def BrigthnessNormalization(img):
r = img / np.sqrt( np.sum((img**2), 0) )
return r
def explained_variance(img):
from sklearn.decomposition import PCA
w = ns.Whiten()
wdata = w.apply(img)
numBands = r.count
h, w, numBands = wdata.shape
X = np.reshape(wdata, (w*h, numBands))
pca = PCA()
mnf = pca.fit_transform(X)
return print(np.cumsum(np.round(pca.explained_variance_ratio_, decimals=4)*100))
def reshape_as_image(arr):
"""Returns the source array reshaped into the order
expected by image processing and visualization software
(matplotlib, scikit-image, etc)
by swapping the axes order from (bands, rows, columns)
to (rows, columns, bands)
Parameters
----------
source : array-like in a of format (bands, rows, columns)
"""
# swap the axes order from (bands, rows, columns) to (rows, columns, bands)
im = np.ma.transpose(arr, [1,2,0])
return im
def reshape_as_raster(arr):
"""Returns the array in a raster order
by swapping the axes order from (rows, columns, bands)
to (bands, rows, columns)
Parameters
----------
arr : array-like in the image form of (rows, columns, bands)
"""
# swap the axes order from (rows, columns, bands) to (bands, rows, columns)
im = np.transpose(arr, [2,0,1])
return im
def saveMNF(img, inputRaster):
# Save TIF image to a nre directory of name MNF
img2 = reshape_as_raster(img)
if args["preprop"]==True:
output = "BN_MNF/" + name[:-4] + "_BN_MNF.tif"
else:
output = "MNF/" + name[:-4] + "_MNF.tif"
new_dataset = rasterio.open(output, 'w', driver='GTiff',
height=inputRaster.shape[0], width=inputRaster.shape[1],
count=int(n_components), dtype=str(img.dtype),
crs=inputRaster.crs, transform=inputRaster.transform)
new_dataset.write(img2)
new_dataset.close()
###############################
if __name__ == "__main__":
# create the arguments for the algorithm
parser = argparse.ArgumentParser()
parser.add_argument('-i','--inputImage',
help='Input raster from which copy the MNF coefficients', type=str)
parser.add_argument('-c','--components',
help='Number of components', type=int, required=True)
parser.add_argument('-m','--method',
help='MNF method to apply: 1 (default) = regular MNF transformation; 2 = MNF invers transformation',
type=int, default=1)
parser.add_argument('-p','--preprop',
help='Preprocessing: Brightness Normalization of Hyperspectral data [Optional]',
action="store_true", default=False)
parser.add_argument('-s','--SavitzkyGolay',
help='Apply Savitzky Golay filtering [Optional]', action="store_true", default=False)
parser.add_argument('-v','--variance',
help='Accumulated explained variance', action="store_true", default=False)
parser.add_argument('--version', action='version', version='%(prog)s 1.0')
args = vars(parser.parse_args())
# Define number of components for the MNF
n_components = args['components']
# Input image
inImage = args['inputImage']
# list of .tif files in the Input File Path
imageList = glob.glob('*.'+inImage[-3:])
imageList.remove(inImage) # remove the Input image
# Create folders to store results if thay do no exist
if args["preprop"]==True:
if not os.path.exists("BN_MNF"):
os.makedirs("BN_MNF")
else:
if not os.path.exists("MNF"):
os.makedirs("MNF")
if args['variance']==True:
# Show the accumulated explained variance
name = os.path.basename(inImage)
r = rasterio.open(inImage)
r2 = r.read()
# Apply Brightness Normalization if the option -p is added
if args["preprop"]==True:
bn = np.apply_along_axis(BrigthnessNormalization, 0, r2)
r2 = reshape_as_image(bn)
img = reshape_as_image(r2)
print("Accumulated explained variances of " + name + "are:")
explained_variance(img)
else:
if args['method']==1:
# Load raster/convert to ndarray format
name = os.path.basename(inImage)
r = rasterio.open(inImage)
r2 = r.read()
# Apply Brightness Normalization if the option -p is added
if args["preprop"]==True:
r2 = np.apply_along_axis(BrigthnessNormalization, 0, r2)
img = reshape_as_image(r2)
# Apply MNF -m 1
print("Creating MNF components of " + name)
# Apply MNF namualy acording to pysptools
w = ns.Whiten()
wdata = w.apply(img)
numBands = r.count
h, w, numBands = wdata.shape
X = np.reshape(wdata, (w*h, numBands))
pca = PCA()
mnf = pca.fit_transform(X)
mnf = np.reshape(mnf, (h, w, numBands))
if args["SavitzkyGolay"]==True:
dn = ns.SavitzkyGolay()
mnf[:,:,1:2] = dn.denoise_bands(mnf[:,:,1:2], 15, 2)
mnf = mnf[:,:,:n_components]
saveMNF(mnf, r)
# Apply MNF coefficients to the other images
for i in range(len(imageList)):
name = os.path.basename(imageList[i])
r = rasterio.open(imageList[i])
r2 = r.read()
# Apply Brightness Normalization if the option -p is added
if args["preprop"]==True:
r2 = np.apply_along_axis(BrigthnessNormalization, 0, r2)
img = reshape_as_image(r2)
# Apply MNF -m 1
print("Creating MNF components of " + name)
# Apply MNF namualy acording to pysptools
w = ns.Whiten()
wdata = w.apply(img)
numBands = r.count
h, w, numBands = wdata.shape
Y = np.reshape(wdata, (w*h, numBands))
mnf = pca.fit_transform(Y)
mnf = np.reshape(mnf, (h, w, numBands))
if args["SavitzkyGolay"]==True:
dn = ns.SavitzkyGolay()
mnf[:,:,1:2] = dn.denoise_bands(mnf[:,:,1:2], 15, 2)
mnf = mnf[:,:,:n_components]
saveMNF(mnf, r)
elif args['method']==2:
# Load raster/convert to ndarray format
name = os.path.basename(inImage)
r = rasterio.open(inImage)
r2 = r.read()
# Apply Brightness Normalization if the option -p is added
if args["preprop"]==True:
r2 = np.apply_along_axis(BrigthnessNormalization, 0, r2)
img = reshape_as_image(r2)
# Apply MNF -m 1
print("Creating MNF components of " + name)
# Apply MNF namualy acording to pysptools
w = ns.Whiten()
wdata = w.apply(img)
numBands = r.count
h, w, numBands = wdata.shape
X = np.reshape(wdata, (w*h, numBands))
pca = PCA()
mnf = pca.fit_transform(X)
mnf = np.reshape(mnf, (h, w, numBands))
if args["SavitzkyGolay"]==True:
dn = ns.SavitzkyGolay()
mnf[:,:,1:2] = dn.denoise_bands(mnf[:,:,1:2], 15, 2)
a = pca.inverse_transform(mnf)
mnf = a[:,:,1:n_components+1]
saveMNF(mnf, r)
# Apply MNF coefficients to the other images
for i in range(len(imageList)):
name = os.path.basename(imageList[i])
r = rasterio.open(imageList[i])
r2 = r.read()
# Apply Brightness Normalization if the option -p is added
if args["preprop"]==True:
r2 = np.apply_along_axis(BrigthnessNormalization, 0, r2)
img = reshape_as_image(r2)
# Apply MNF -m 1
print("Creating MNF components of " + name)
# Apply MNF namualy acording to pysptools
w = ns.Whiten()
wdata = w.apply(img)
numBands = r.count
h, w, numBands = wdata.shape
Y = np.reshape(wdata, (w*h, numBands))
mnf = pca.fit_transform(Y)
mnf = np.reshape(mnf, (h, w, numBands))
if args["SavitzkyGolay"]==True:
dn = ns.SavitzkyGolay()
mnf[:,:,1:2] = dn.denoise_bands(mnf[:,:,1:2], 15, 2)
a = pca.inverse_transform(mnf)
mnf = a[:,:,1:n_components+1]
saveMNF(mnf, r)
else:
print('ERROR!. Command should have the form:')
print('python MNF.py -f <Imput raster formar> -c <Number of components> -m <Method option>[optional] -v <Accumulated explained variance>[optional]')
print("")
print("Method options: 1 (default) regular MNF transformation")
print(" 2 inverse transform")
print("")
print("-p or --preprop: Apply Broghtness Normalization of hyperspectral data")
print("")
print("-s or --Savitzky Golay: Use Savitzky Golay methods")
print("")
print("example: python MNF_cmd.py -f tif -c 10")