spectralpython / spectral

Python module for hyperspectral image processing
MIT License
563 stars 138 forks source link

Update continuum.py #147

Open Russjas opened 1 year ago

Russjas commented 1 year ago

I am new to github, opensource, python.. a lot of things, so if this is inappropriate, please dont take offense and tell me.

I have written an application of the SQM method for returning true wavelength position of the minimum of an absorbtion feature, the reference is within the docstring. I have been using this with the continuum module of spectral python, and thought it might be something you could use.

I am new to coding, so this may not be well enough written for you to include it, but you have the option.

Again, apologies if any of this is inappropriate or unwelcome

tboggs commented 1 year ago

Hi Russjas. Thanks for the submission. I think this could be included in the package. Could you please provide an example or two of the function being used (e.g. as a comment in this PR, with associated figure).

Also, could you please update the PR to pull into the develop branch (instead of master)?

Thanks, Thomas

Russjas commented 1 year ago

No problem, here is an example. Data is from the Geological Survey Ireland core scanner, and is freely available

Open data as normal using spectral python. I have skipped reflectance correction and normalisation steps and just start with normalised reflectance data

box = envi.open('C:/Users/Russell.Rogers/Documents/Hyperspectral/Datafiles/Reflectance/Corrected_RR/ref_16-2_22/denoisedRun2_ref_16-2_22.hdr')
data = np.array(box.load())
bands = np.array(box.bands.centers)
view =  sp.imshow(data)

Core_run

I use the remove continuum function from spectral-python on a subset of bands. in this instance I know that my data are carbonates and I want to target the 2340 absorbtion feature related to the -CO3 bond

bands = bands[220:]
Cr_tight = sp.remove_continuum(data[:, :, 220:], bands)

mins = np.argmin(Cr_tight, axis=2)
view = sp.imshow(data, (172, 219, 244), classes=mins)

print(bands[mins[20, 30]])

Prints: 2341.76, which the centre of the wavelength band. In this case, nearly all the data has the minima in the same wavelength band argmin

However once sqm is calculated for each pixel, they can be binned into, for example, 10 bins for use in sp.imshow as classes

sqm = get_SQM(Cr_tight, bands)
wav = sqm.wavelengths
dep = sqm.depths

bins = np.arange(np.min(wav), np.max(wav), ((np.max(wav) - np.min(wav)) / 10))
x = np.digitize(wav, bins)
view = sp.imshow(data, classes=x)

print(wav[20, 30])

prints: 2340.149882593446. This is the same pixel that I printed the argmin for above, but the wavelength of the minimum position is different after SQM SQM

And the classes are more useful! In this case the lower classes are indicative of a more dolomitic composition, while higher classes are more calcitic