Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

i.atcorr: major refactoring of create_iwave.py #3886

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
113 changes: 47 additions & 66 deletions imagery/i.atcorr/create_iwave.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,19 +34,16 @@

def usage():
"""How to use this..."""
print("create_iwave.py <csv file>")
print
print("create_iwave.py <csv file>\n")
print("Generates filter function template for iwave.cpp from csv file. Note:")
print("- csv file must have wl response for each band in each column")
print("- first line must be a header with wl followed by band names")
print("- all following lines will be the data.")
print("If spectral response is null, leave field empty in csv file. Example:")
print
print("If spectral response is null, leave field empty in csv file. Example:\n")
print("WL(nm),band 1,band 2,band 3,band 4")
print("455,0.93,,,")
print("485,0.94,0.00,,")
print("545,0.00,0.87,0.00,")
print
print("545,0.00,0.87,0.00,\n")
print("This script will interpolate the filter functions to 2.5 nm steps")
print("and output a cpp template file in the IWave format to be added to iwave.cpp")

Expand All @@ -62,20 +59,15 @@ def read_input(csvfile):
first column is wavelength
values are those of the discrete band filter functions
"""
infile = open(csvfile, "r")
with open(csvfile, "r") as infile:
# get number of bands and band names
bands = infile.readline().strip().split(",")[1:]

# get number of bands and band names
bands = infile.readline().split(",")
bands.remove(bands[0])
bands[-1] = bands[-1].strip()
print(" > Number of bands found: %d" % len(bands))
infile.close()
print(f" > Number of bands found: {len(bands)}")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A linter should complain that bands could be undefined here.


# create converter dictionary for import
# fix nodata or \n
conv = {}
for b in range(len(bands)):
conv[b + 1] = lambda s: float(s or 0)
conv = {b + 1: lambda s: float(s or 0) for b in range(len(bands))}

values = np.loadtxt(csvfile, delimiter=",", skiprows=1, converters=conv)

Expand All @@ -90,10 +82,8 @@ def interpolate_band(values, step=2.5):
and min, max wl values
values must be numpy array with 2 columns
"""

# removing nodata and invalid values
w = values[:, 1] >= 0
values_clean = values[w]
values_clean = values[values[:, 1] >= 0]

wavelengths = values_clean[:, 0] # 1st column of input array
responses = values_clean[:, 1] # 2nd column
Expand Down Expand Up @@ -176,25 +166,25 @@ def pretty_print(filter_f):
Create pretty string out of filter function
8 values per line, with spaces, commas and all the rest
"""
pstring = ""
pstring = []
for i in range(len(filter_f) + 1):
if i % 8 is 0:
if i is not 0:
if i % 8 == 0:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just a infinitesimal implementation-dependent optimisation possible here: I stumbled upon something on stack overflow this weekend, not exactly related to this, but the example had this form.

This pattern uses at least 2 more opcodes than another possibility, as it has to protect himself from the variable not being an integer when comparing to 0. There's an alternative, but it was open three days ago on my computer, and I'm writing on my phone. It had something to do with either using "in" or nothing at all.

I'm just nerdsniping here, if you are curious to microoptimise here.

if i != 0:
value_wo_leading_zero = ("%.4f" % (filter_f[i - 1])).lstrip("0")
pstring += value_wo_leading_zero
if i > 1 and i < len(filter_f):
pstring += ", "
if i is not 1:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The replacement of this doesn't seem exactly the equivalent. Is the new logic expected behavior? (Old one got triggered for -1, 0, 2, 3, 4..., new gets triggered for 2, 3, 4...)

pstring.append(value_wo_leading_zero)
if i > 1:
if i < len(filter_f):
pstring.append(", ")
# trim the trailing whitespace at the end of line
pstring = pstring.rstrip()
pstring += "\n "
pstring[-1] = pstring[-1].rstrip()
pstring.append("\n ")
else:
value_wo_leading_zero = ("%.4f" % (filter_f[i - 1])).lstrip("0")
pstring += value_wo_leading_zero
pstring.append(value_wo_leading_zero)
if i < len(filter_f):
pstring += ", "
pstring.append(", ")
# trim starting \n and trailing ,
pstring = pstring.lstrip("\n").rstrip(", ")
pstring = "".join(pstring).lstrip("\n").rstrip(", ")
return pstring


Expand All @@ -205,9 +195,26 @@ def write_cpp(bands, values, sensor, folder):
needs other functions: interpolate_bands, pretty_print
"""

def get_min_wavelength(c, rthresh, fi):
echoix marked this conversation as resolved.
Show resolved Hide resolved
"""Get minimum wavelength rounded by threshold.

:param fi: filter function
"""
while c > 0 and fi[c - 1] > rthresh:
c = c - 1
return np.ceil(li[0] * 1000 + (2.5 * c))

def get_max_wavelength(c, rthresh, fi):
"""Get maximum wavelength rounded by threshold.

:param fi: filter function
"""
while c < len(fi) - 1 and fi[c + 1] > rthresh:
c = c + 1
return np.floor(li[0] * 1000 + (2.5 * c))

# keep in sync with IWave::parse()
rthresh = 0.01
print
print(" > Response peaks from interpolation to 2.5 nm steps:")

# getting necessary data
Expand All @@ -219,18 +226,8 @@ def write_cpp(bands, values, sensor, folder):
li = limits
# Get wavelength range for spectral response in band
maxresponse_idx = np.argmax(fi)
# Get minimum wavelength with spectral response
c = maxresponse_idx
while c > 0 and fi[c - 1] > rthresh:
c = c - 1
min_wavelength = np.ceil(li[0] * 1000 + (2.5 * c))
# Get maximum wavelength with spectral response
c = maxresponse_idx
while c < len(fi) - 1 and fi[c + 1] > rthresh:
c = c + 1
max_wavelength = np.floor(li[0] * 1000 + (2.5 * c))
print(" %s (%inm - %inm)" % (bands[b], min_wavelength, max_wavelength))

min_wavelength = get_min_wavelength(maxresponse_idx, rthresh, fi)
max_wavelength = get_max_wavelength(maxresponse_idx, rthresh, fi)
else:
filter_f = []
limits = []
Expand All @@ -241,29 +238,17 @@ def write_cpp(bands, values, sensor, folder):

# Get wavelength range for spectral response in band
maxresponse_idx = np.argmax(fi)
# Get minimum wavelength with spectral response
c = maxresponse_idx
while c > 0 and fi[c - 1] > rthresh:
c = c - 1
min_wavelength = np.ceil(li[0] * 1000 + (2.5 * c))
# Get maximum wavelength with spectral response
c = maxresponse_idx
while c < len(fi) - 1 and fi[c + 1] > rthresh:
c = c + 1
max_wavelength = np.floor(li[0] * 1000 + (2.5 * c))
min_wavelength = get_min_wavelength(maxresponse_idx, rthresh, fi)
max_wavelength = get_max_wavelength(maxresponse_idx, rthresh, fi)
print(" %s (%inm - %inm)" % (bands[b], min_wavelength, max_wavelength))

# writing...
outfile = open(os.path.join(folder, sensor + "_cpp_template.txt"), "w")
outfile.write("/* Following filter function created using create_iwave.py */\n\n")

if len(bands) == 1:
outfile.write("void IWave::%s()\n{\n\n" % (sensor.lower()))
else:
outfile.write("void IWave::%s(int iwa)\n{\n\n" % (sensor.lower()))

# single band case
if len(bands) == 1:
outfile.write("void IWave::%s()\n{\n\n" % (sensor.lower()))
outfile.write(" /* %s of %s */\n" % (bands[0], sensor))
outfile.write(" static const float sr[%i] = {" % (len(filter_f)))
filter_text = pretty_print(filter_f)
Expand All @@ -289,6 +274,7 @@ def write_cpp(bands, values, sensor, folder):
outfile.write("}\n")

else: # more than 1 band
outfile.write("void IWave::%s(int iwa)\n{\n\n" % (sensor.lower()))
# writing bands
for b in range(len(bands)):
outfile.write(" /* %s of %s */\n" % (bands[b], sensor))
Expand All @@ -299,9 +285,8 @@ def write_cpp(bands, values, sensor, folder):
outfile.write(filter_text + "\n };\n\t\n")

# writing band limits
for b in range(len(bands)):
inf = ", ".join(["%.4f" % i[0] for i in limits])
sup = ", ".join(["%.4f" % i[1] for i in limits])
inf = ", ".join(["%.4f" % i[0] for i in limits])
sup = ", ".join(["%.4f" % i[1] for i in limits])

outfile.write(" static const float wli[%i] = {%s};\n" % (len(bands), inf))
outfile.write(" static const float wls[%i] = {%s};\n" % (len(bands), sup))
Expand Down Expand Up @@ -334,7 +319,6 @@ def main():
# getting sensor name from full csv file name
sensor = os.path.splitext(os.path.basename(inputfile))[0]

print
print(" > Getting sensor name from csv file: %s" % (sensor))

# getting data from file
Expand All @@ -344,7 +328,6 @@ def main():
# consider only wavelengths with a reasonably large response
# around the peak response, keep in sync with IWave::parse()
rthresh = 0.01
print
print(" > Response peaks from input file:")
for b in range(1, len(bands) + 1):
lowl = 0
Expand Down Expand Up @@ -372,7 +355,6 @@ def main():
# writing file in same folder of input file
write_cpp(bands, values, sensor, os.path.dirname(inputfile))

print
print(
" > Filter functions exported to %s"
% ("sensors_csv/" + sensor + "_cpp_template.txt")
Expand All @@ -385,7 +367,6 @@ def main():
" > Don't forget to add the necessary data to the files"
" iwave.h, geomcond.h, geomcond.cpp, and to i.atcorr.html"
)
print


if __name__ == "__main__":
Expand Down
Loading