forked from pmlrsg/arsf_tools
-
Notifications
You must be signed in to change notification settings - Fork 0
/
apply_elc_to_bil.py
executable file
·177 lines (136 loc) · 6.19 KB
/
apply_elc_to_bil.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
#!/usr/bin/env python
"""
Empirical Line Calibration
Takes field and image spectra from white, grey and black target in a CSV file
with the following headings:
wavelength,white,grey,black
From these performs a linar fit for each wavelength in the image (interpolating
field spectra to match).
Coefficients from linear fit are then applied to the image.
Author: Dan Clewley
Creation Date: 2018-02-20
"""
from __future__ import print_function
import argparse
import time
import numpy
from arsf_envi_reader import numpy_bin_reader
from arsf_envi_reader import envi_header
DEFAULT_SCALE_FACTOR = 10000
def get_elc_coefficients(image_spectra_file, field_spectra_file):
"""
Takes images spectra and field spectra from white, grey and black
targets and performs a linear fit for each band.
Returns numpy array containing two columns with ax and bx
"""
image_spectra = numpy.genfromtxt(image_spectra_file,
delimiter=",", names=True)
field_spectra = numpy.genfromtxt(field_spectra_file,
delimiter=",", names=True)
elc_coefficients = []
for i, wavelength in enumerate(image_spectra["wavelength"]):
# Get radiance from image spectra for given wavelength
image_white = image_spectra["white"][i]
image_grey = image_spectra["grey"][i]
image_black = image_spectra["black"][i]
# Get coresponding values from field spectra
# As wavelengths are likely to not exactly match up
# use interpolation to get value for image wavelength
field_white = numpy.interp(wavelength, field_spectra["wavelength"],
field_spectra["white"])
field_grey = numpy.interp(wavelength, field_spectra["wavelength"],
field_spectra["grey"])
field_black = numpy.interp(wavelength, field_spectra["wavelength"],
field_spectra["black"])
# Perform linear fit
ax, bx = numpy.polyfit([image_black, image_grey, image_white],
[field_black, field_grey, field_white],
1)
elc_coefficients.append([ax, bx])
elc_coefficients = numpy.array(elc_coefficients)
return elc_coefficients
def apply_elc(input_image, output_image, image_spectra_file,
field_spectra_file,
scale_factor=DEFAULT_SCALE_FACTOR):
"""
Derive coefficients for an emprical line calibration and apply
to a BIL file.
"""
print("Getting coefficients")
elc_coefficients = get_elc_coefficients(image_spectra_file,
field_spectra_file)
print("Applying to image")
# Read ENVI header
input_header_file = envi_header.find_hdr_file(input_image)
input_header_dict = envi_header.read_hdr_file(input_header_file)
# Open input file
in_data = numpy_bin_reader.BilReader(input_image)
# Open file for output
out_file = open(output_image, "wb")
for i, line in enumerate(in_data):
print(" Line {0:05}/{1:05}".format(i, in_data.lines), end="\r")
# Apply coefficients
elc_line = elc_coefficients[:,0] * line.transpose() \
+ elc_coefficients[:,1]
elc_line = elc_line.transpose()
# Set any negative values to 0
elc_line = numpy.where(elc_line < 0, 0, elc_line)
elc_line = elc_line * scale_factor
# For default scale factor is over 1000 express as integer
if scale_factor >= 1000:
elc_line = elc_line.astype(numpy.uint16)
else:
elc_line = elc_line.astype(numpy.float32)
elc_line.tofile(out_file)
print(" Line {0:05}/{0:05}".format(in_data.lines), end="\n")
output_header_dict = input_header_dict
# See if need to change data type
if scale_factor >= 1000:
output_header_dict["data type"] = 12
else:
output_header_dict["data type"] = 4
# Remove radiance data units
try:
output_header_dict.pop("radiance data units")
except KeyError:
pass
# Add scale factor
output_header_dict["reflectance scale factor"] = scale_factor
procesing_comment = """\nConverted to reflectance using empirical line method. Field spectra: {}, Image spectra: {}\n""".format(field_spectra_file, image_spectra_file)
try:
output_header_dict["_comments"] = output_header_dict["_comments"] + \
procesing_comment
except KeyError:
output_header_dict["_comments"] = procesing_comment
output_header_dict["description"] = "BIL file created by {} on {}".format(__file__,
time.strftime("%Y-%m-%d %H:%M:%S"))
# Copy header
envi_header.write_envi_header(output_image + ".hdr",
output_header_dict)
# Close files
out_file.close()
in_data = None
if __name__ == "__main__":
parser = argparse.ArgumentParser("Use an empirical line calibration to "
"convert image spectra to surface "
"reflectance")
parser.add_argument("inputimage", nargs=1,
type=str, help="Input image")
parser.add_argument("outputimage", nargs=1,
type=str, help="Output image")
parser.add_argument("--image_spectra",
required=True,
type=str, help="CSV file with target spectra from images")
parser.add_argument("--field_spectra",
required=True,
type=str, help="CSV file with target spectra from field")
parser.add_argument("--scale", type=float,
required=False,
default=DEFAULT_SCALE_FACTOR,
help="Scale factor to apply to reflectance values. "
"Used so image can be saved as uint16. "
"(Default {}".format(DEFAULT_SCALE_FACTOR))
args = parser.parse_args()
apply_elc(args.inputimage[0], args.outputimage[0],
args.image_spectra, args.field_spectra)
print("Saved to {}".format(args.outputimage[0]))