diff --git a/.gitignore b/.gitignore index 74ebc99..00747e1 100644 --- a/.gitignore +++ b/.gitignore @@ -2,5 +2,8 @@ __pycache__* *~ *.pyc /dist/ +/build/ +/.vscode/ /*.egg-info -*.ipynb_checkpoints* \ No newline at end of file +*.ipynb_checkpoints* +.DS_Store \ No newline at end of file diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..c35a428 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,2 @@ +# we include the rsr +include specdal/rsr/*.nc \ No newline at end of file diff --git a/setup.py b/setup.py index 531398d..a0ee5e8 100644 --- a/setup.py +++ b/setup.py @@ -15,6 +15,7 @@ author_email='ylee546@wisc.edu', license='MIT', packages=['specdal', 'specdal.utils'], + include_package_data=True, install_requires=['numpy', 'pandas', 'matplotlib', 'scipy'], zip_safe=False, classifiers=[ diff --git a/specdal/collection.py b/specdal/collection.py index 7dd16d3..c543fd1 100644 --- a/specdal/collection.py +++ b/specdal/collection.py @@ -93,14 +93,14 @@ class Collection(object): """ Represents a dataset consisting of a collection of spectra """ - def __init__(self, name, directory=None, spectra=None, - measure_type='pct_reflect', metadata=None): + def __init__(self, name, directory=None, spectra=None, ext=[".asd", ".sed", ".sig"], + measure_type='pct_reflect', metadata=None, reader=None): self.name = name self.spectra = spectra self.measure_type = measure_type self.metadata = metadata if directory: - self.read(directory, measure_type) + self.read(directory, measure_type, ext=ext, reader=reader) @property def spectra(self): """ @@ -187,7 +187,7 @@ def __contains__(self, item): # reader def read(self, directory, measure_type='pct_reflect', ext=[".asd", ".sed", ".sig"], recursive=False, - verbose=False): + verbose=False, reader=None): """ read all files in a path matching extension """ @@ -204,7 +204,7 @@ def read(self, directory, measure_type='pct_reflect', filepath = os.path.join(dirpath, f) spectrum = Spectrum(name=f_name, filepath=filepath, measure_type=measure_type, - verbose=verbose) + verbose=verbose, reader=reader) self.append(spectrum) ################################################## # wrapper around spectral operations @@ -318,4 +318,37 @@ def std(self, append=False): if append: self.append(spectrum) return spectrum + ################################################## + # method for computing the values for a specific satellite + + def getSatellite(self, satellite="Aqua", sensor="MODIS", rsr_path = __file__.replace("/collection.py","/rsr/")): + c_tmp = Collection(name=self.name, metadata={}) + c_tmp.metadata["satellite"] = satellite + c_tmp.metadata["sensor"] = sensor + # compute reflectance by bande + size_compute = len(self.spectra) + i = 1 + # We iterate over all spectra to compute the reflectance per band + for spectra_tmp in self.spectra: + # we print current spectra + print(f"Spectra {i}/{size_compute}.", end="\r") + c_tmp.append(spectra_tmp.getSatellite(satellite, sensor, rsr_path)) + i+=1 + + return c_tmp + + def savgol_filter(self, window_length, polyorder, deriv=0, + delta=1.0, axis=-1, mode='interp', cval=0.0): + self.metadata["savgol_window_length"] = window_length + self.metadata["savgol_polyorder"] = polyorder + + # We iterate over all spectra + for spectra_tmp in self.spectra: + spectra_tmp.savgol_filter(window_length, + polyorder, deriv, delta, axis, mode, cval) + def walevength_range(self, wlmin=350, wlmax=2500, dtype=None): + # We iterate over all spectra + for spectra_tmp in self.spectra: + spectra_tmp.measurement = spectra_tmp.measurement.loc[wlmin:wlmax] + spectra_tmp.metadata['wavelength_range'] = (wlmin, wlmax) diff --git a/specdal/operator.py b/specdal/operator.py index 9cda34d..f3dea1b 100644 --- a/specdal/operator.py +++ b/specdal/operator.py @@ -2,6 +2,7 @@ # wavelength as index and measurement as values import pandas as pd import numpy as np +from scipy.signal import savgol_filter from .utils.misc import get_monotonic_series ################################################################################ @@ -114,6 +115,20 @@ def vector_normalize(series): ''' return series / np.sqrt(np.sum(series ** 2)) +################################################################################ +# vector_normalize: divide the series such that 2-norm(series) == 1 +def savgol(series, window_length, polyorder, deriv=0, + delta=1.0, axis=-1, mode='interp', cval=0.0): + ''' + Savitzky-Golay filter using scipy implementation + ''' + smooth = savgol_filter(series.values, + window_length, polyorder, deriv, + delta, axis, mode, cval) + smooth = pd.Series(data=smooth, index=series.index) + smooth.index.name = "wavelength" + return smooth + ################################################################################ # DataFrame operations (collection of spectra) ################################################################################ diff --git a/specdal/reader.py b/specdal/reader.py index c0c6d51..579742b 100644 --- a/specdal/reader.py +++ b/specdal/reader.py @@ -1,19 +1,25 @@ # readers.py provides functions to read spectrum files for data and # metadata. +import string import pandas as pd import numpy as np from os.path import abspath, expanduser, splitext from collections import OrderedDict from .utils.reader_utils import * -def read(filepath, read_data=True, read_metadata=True, verbose=False): +def read(filepath, read_data=True, read_metadata=True, verbose=False, reader=None): """ Calls the appropriate reader based on file extension """ SUPPORTED_READERS = {'.asd':read_asd , '.sig':read_sig , '.sed':read_sed } - ext = splitext(filepath)[1] - reader = SUPPORTED_READERS[ext] + # Added option to force reader + if reader == None: + ext = splitext(filepath)[1] + reader = SUPPORTED_READERS[ext] + else: + reader = SUPPORTED_READERS[reader] + return reader(abspath(expanduser(filepath)), read_data, read_metadata, verbose) @@ -193,12 +199,24 @@ def read_asd(filepath, read_data=True, read_metadata=True, verbose=False): gps_flags2 = gps_struct[9:11] # unpack this into bits gps_satellites = gps_struct[11:16] gps_filler = gps_struct[16:18] + # date + time_lock = 160 + second = int.from_bytes(binconts[160:(160 + 2)], byteorder='little') + minute = int.from_bytes(binconts[160+2:160+2+2], byteorder='little') + hour = int.from_bytes(binconts[160+4:160+4+2], byteorder='little') + day = int.from_bytes(binconts[160+6:160+6+2], byteorder='little') + month = int.from_bytes(binconts[160+8:160+8+2], byteorder='little')+1 + year = int.from_bytes(binconts[160+10:160+10+2], byteorder='little')+1900 + # Dark current + dc = int.from_bytes(binconts[181:181+1], byteorder='little') # metadata metadata['integration_time'] = integration_time metadata['measurement_type'] = spectrum_type metadata['gps_time_tgt'] = gps_timestamp metadata['gps_time_ref'] = None metadata['wavelength_range'] = (wavestart, wavestop) + metadata['measurement_date'] = f"{year:04}-{month:02}-{day:02} {hour:02}:{minute:02}:{second:02}" + metadata['dark_current'] = dc # metadata['splice'] = (splice1, splice2) # metadata['resolution'] = wavestep return data, metadata diff --git a/specdal/rsr/adeos_octs_RSR.nc b/specdal/rsr/adeos_octs_RSR.nc new file mode 100644 index 0000000..dacbe44 Binary files /dev/null and b/specdal/rsr/adeos_octs_RSR.nc differ diff --git a/specdal/rsr/aqua_modis_RSR.nc b/specdal/rsr/aqua_modis_RSR.nc new file mode 100644 index 0000000..c474f4b Binary files /dev/null and b/specdal/rsr/aqua_modis_RSR.nc differ diff --git a/specdal/rsr/aviris_RSR.nc b/specdal/rsr/aviris_RSR.nc new file mode 100644 index 0000000..0439838 Binary files /dev/null and b/specdal/rsr/aviris_RSR.nc differ diff --git a/specdal/rsr/coms_goci_RSR.nc b/specdal/rsr/coms_goci_RSR.nc new file mode 100644 index 0000000..ffdd127 Binary files /dev/null and b/specdal/rsr/coms_goci_RSR.nc differ diff --git a/specdal/rsr/envisat_meris_RSR.nc b/specdal/rsr/envisat_meris_RSR.nc new file mode 100644 index 0000000..a02b77e Binary files /dev/null and b/specdal/rsr/envisat_meris_RSR.nc differ diff --git a/specdal/rsr/gcomc_sgli_RSR.nc b/specdal/rsr/gcomc_sgli_RSR.nc new file mode 100644 index 0000000..7ccb2b1 Binary files /dev/null and b/specdal/rsr/gcomc_sgli_RSR.nc differ diff --git a/specdal/rsr/irs-p4_ocm_RSR.nc b/specdal/rsr/irs-p4_ocm_RSR.nc new file mode 100644 index 0000000..ab7606c Binary files /dev/null and b/specdal/rsr/irs-p4_ocm_RSR.nc differ diff --git a/specdal/rsr/iss_hico_RSR.nc b/specdal/rsr/iss_hico_RSR.nc new file mode 100644 index 0000000..c0f69f6 Binary files /dev/null and b/specdal/rsr/iss_hico_RSR.nc differ diff --git a/specdal/rsr/jpss-1_viirs_RSR.nc b/specdal/rsr/jpss-1_viirs_RSR.nc new file mode 100644 index 0000000..2247bba Binary files /dev/null and b/specdal/rsr/jpss-1_viirs_RSR.nc differ diff --git a/specdal/rsr/landsat-8_oli_RSR.nc b/specdal/rsr/landsat-8_oli_RSR.nc new file mode 100644 index 0000000..88a5878 Binary files /dev/null and b/specdal/rsr/landsat-8_oli_RSR.nc differ diff --git a/specdal/rsr/nimbus-7_czcs_RSR.nc b/specdal/rsr/nimbus-7_czcs_RSR.nc new file mode 100644 index 0000000..274650e Binary files /dev/null and b/specdal/rsr/nimbus-7_czcs_RSR.nc differ diff --git a/specdal/rsr/oceansat-2_ocm-2_RSR.nc b/specdal/rsr/oceansat-2_ocm-2_RSR.nc new file mode 100644 index 0000000..2c81a0c Binary files /dev/null and b/specdal/rsr/oceansat-2_ocm-2_RSR.nc differ diff --git a/specdal/rsr/orbview-2_seawifs_RSR.nc b/specdal/rsr/orbview-2_seawifs_RSR.nc new file mode 100644 index 0000000..4734938 Binary files /dev/null and b/specdal/rsr/orbview-2_seawifs_RSR.nc differ diff --git a/specdal/rsr/pace_oci_RSR.nc b/specdal/rsr/pace_oci_RSR.nc new file mode 100644 index 0000000..c5240f4 Binary files /dev/null and b/specdal/rsr/pace_oci_RSR.nc differ diff --git a/specdal/rsr/s3a_olci_RSR.nc b/specdal/rsr/s3a_olci_RSR.nc new file mode 100644 index 0000000..c88f22e Binary files /dev/null and b/specdal/rsr/s3a_olci_RSR.nc differ diff --git a/specdal/rsr/s3b_olci_RSR.nc b/specdal/rsr/s3b_olci_RSR.nc new file mode 100644 index 0000000..8f28624 Binary files /dev/null and b/specdal/rsr/s3b_olci_RSR.nc differ diff --git a/specdal/rsr/seahawk1_hawkeye_RSR.nc b/specdal/rsr/seahawk1_hawkeye_RSR.nc new file mode 100644 index 0000000..cec0dc1 Binary files /dev/null and b/specdal/rsr/seahawk1_hawkeye_RSR.nc differ diff --git a/specdal/rsr/sentinel-2a_msi_RSR.nc b/specdal/rsr/sentinel-2a_msi_RSR.nc new file mode 100644 index 0000000..1d3bbf8 Binary files /dev/null and b/specdal/rsr/sentinel-2a_msi_RSR.nc differ diff --git a/specdal/rsr/sentinel-2b_msi_RSR.nc b/specdal/rsr/sentinel-2b_msi_RSR.nc new file mode 100644 index 0000000..8b4a2e7 Binary files /dev/null and b/specdal/rsr/sentinel-2b_msi_RSR.nc differ diff --git a/specdal/rsr/suomi-npp_viirs_RSR.nc b/specdal/rsr/suomi-npp_viirs_RSR.nc new file mode 100644 index 0000000..678f31c Binary files /dev/null and b/specdal/rsr/suomi-npp_viirs_RSR.nc differ diff --git a/specdal/rsr/terra_misr_RSR.nc b/specdal/rsr/terra_misr_RSR.nc new file mode 100644 index 0000000..69be74a Binary files /dev/null and b/specdal/rsr/terra_misr_RSR.nc differ diff --git a/specdal/rsr/terra_modis_RSR.nc b/specdal/rsr/terra_modis_RSR.nc new file mode 100644 index 0000000..0822047 Binary files /dev/null and b/specdal/rsr/terra_modis_RSR.nc differ diff --git a/specdal/spectrum.py b/specdal/spectrum.py index b6ff5f1..9ec02f0 100644 --- a/specdal/spectrum.py +++ b/specdal/spectrum.py @@ -3,13 +3,17 @@ # pandas.Series. import pandas as pd import numpy as np +import xarray import specdal.operator as op from collections import OrderedDict from .reader import read from .utils.misc import get_pct_reflect import os +from numbers import Number +import numpy.lib.mixins +from scipy.integrate import simps -class Spectrum(object): +class Spectrum(numpy.lib.mixins.NDArrayOperatorsMixin): """Class that represents a single spectrum Parameters @@ -38,7 +42,7 @@ def __init__(self, name=None, filepath=None, measurement=None, measure_type='pct_reflect', metadata=None, interpolated=False, stitched=False, jump_corrected=False, vector_normalized=False, derivative_order=0, - verbose=False): + verbose=False, reader=None): if name is None: assert filepath is not None name = os.path.splitext(os.path.basename(filepath))[0] @@ -52,7 +56,7 @@ def __init__(self, name=None, filepath=None, measurement=None, self.vector_normalized = vector_normalized self.derivative_order = derivative_order if filepath: - self.read(filepath, measure_type, verbose=verbose) + self.read(filepath, measure_type, verbose=verbose, reader=reader) def __str__(self): string = "\nname:\t\t{!s},\n".format(self.name) string += "measure_type:\t{!s}\n".format(self.measure_type) @@ -72,16 +76,16 @@ def __str__(self): return string ################################################## # reader - def read(self, filepath, measure_type, verbose=False): + def read(self, filepath, measure_type, verbose=False, reader=None): ''' Read measurement from a file. ''' - data, meta = read(filepath, verbose=verbose) + data, meta = read(filepath, verbose=verbose, reader=reader) self.metadata = meta if measure_type == 'pct_reflect' and 'pct_reflect' not in data: self.measurement = get_pct_reflect(data) return - assert measure_type in data # TODO: handle this + assert measure_type in data # TODO: handle this self.measurement = data[measure_type] ################################################## # wrappers around spectral operations @@ -110,6 +114,16 @@ def derivative(self): ''' self.measurement = op.derivative(self.measurement) self.derivative_order += 1 + + def savgol_filter(self, window_length, polyorder, deriv=0, + delta=1.0, axis=-1, mode='interp', cval=0.0): + ''' + ''' + self.measurement = op.savgol(self.measurement, + window_length, polyorder, deriv, + delta, axis, mode, cval) + self.savgol_window = window_length + self.savgol_polyorder = polyorder ################################################## # wrapper around plot function @@ -120,35 +134,128 @@ def to_csv(self, *args, **kwargs): '''''' return pd.DataFrame(self.measurement).transpose().to_csv( *args, **kwargs) + ################################################## - # wrapper around pandas series operators - def __add__(self, other): + # wrapper for numpy functions + def __array_ufunc__(self, ufunc, method, *inputs, **kwargs): new_measurement = None new_name = self.name + '+' - if isinstance(other, Spectrum): - assert self.measure_type == other.measure_type - new_measurement = self.measurement.__add__(other.measurement).dropna() - new_name += other.name + if method == '__call__': + if self.metadata is None: + metadata = None + else: + metadata = self.metadata + metadata['file'] = None + metadata["measurement_type"]="TRANS_TYPE" + new_name = ufunc.__name__+"(" + values = [] + for input in inputs: + if isinstance(input, Number): + values.append(input) + new_name = new_name +str(input)+", " + elif isinstance(input, Spectrum): + values.append(input.measurement) + new_name = new_name+input.name+", " + else: + return NotImplemented + + new_name = new_name[:-2]+")" + if metadata is not None: + metadata['name'] = new_name + return Spectrum(name=new_name, measurement=ufunc(*values, **kwargs),metadata=metadata, measure_type = 'TRANS_TYPE') else: - new_measurement = self.measurement.__add__(other) - return Spectrum(name=new_name, measurement=new_measurement, - measure_type=self.measure_type) - def __isub__(self, other): - pass - def __imul__(self, other): - pass - def __itruediv__(self, other): - pass - def __ifloordiv__(self, other): - pass - def __iiadd__(self, other): - pass - def __isub__(self, other): - pass - def __imul__(self, other): - pass - def __itruediv__(self, other): - pass - def __ifloordiv__(self, other): - pass + return NotImplemented + + ################################################## + # wrapper for spectral subset + def walevength_range(self, wlmin=350, wlmax=2500, dtype=None): + self.measurement = self.measurement.loc[wlmin:wlmax] + self.metadata['wavelength_range'] = (wlmin, wlmax) + + + ################################################## + # wrapper for array operations + def __array__(self, dtype=None): + return self.measurement.values + + ################################################## + # method for computing the values for a specific satellite + + def getRSR(self, satellite="aqua", sensor="modis", rsr_path=__file__.replace("/spectrum.py","/rsr/")): + # We build a list of available rsr + available_rsr = [x[:-7] for x in os.listdir(rsr_path) if x[0] != "."] + # Build rsr + if sensor == "aviris": + rsr_path = rsr_path+f"{sensor}_RSR.nc" + else: + rsr_path = rsr_path+f"{satellite}_{sensor}_RSR.nc" + # Read bands from dataframe + try: + df = xarray.open_dataset(rsr_path).to_dataframe() + except (FileNotFoundError, IOError): + print(f"Satellite-sensor combination not available. The options are {available_rsr}") + + # Reshape the dataframe as needed + df.reset_index(inplace=True) + df.drop(["wavelengths"], axis=1, inplace=True) + #df.set_index(["bands","wavelength"], inplace=True) + rsr = df.pivot(index="wavelength", columns="bands") + rsr.columns = rsr.columns.droplevel() + rsr.columns.name=None + rsr.index.name = "Wavelength" + # round index to 1 decimal + rsr.index = rsr.index.values.round(0) + # We sort the dataframe + rsr = rsr[sorted(rsr.columns)] + # Remove duplicated indices and sort by index + rsr = rsr.groupby(level=0).sum().sort_index() + + return rsr + + def getSatellite(self, satellite="aqua", sensor="modis", rsr_path = __file__.replace("/spectrum.py","/rsr/")): + # get relative spectral response + rsr = self.getRSR(satellite, sensor, rsr_path) + # compute reflectance by band + ref = rsr.mul(self.measurement, axis='index').sum(axis="index")/rsr.sum(axis="index") + # save to spectrum + name = self.name + spectrum = Spectrum(name=name, measurement=ref, metadata=self.metadata, measure_type=self.measure_type) + + spectrum.metadata["satellite"] = satellite + spectrum.metadata["sensor"] = sensor + + return spectrum + + + + # wrapper around pandas series operators + # def __add__(self, other): + # new_measurement = None + # new_name = self.name + '+' + # if isinstance(other, Spectrum): + # assert self.measure_type == other.measure_type + # new_measurement = self.measurement.__add__(other.measurement).dropna() + # new_name += other.name + # else: + # new_measurement = self.measurement.__add__(other) + # return Spectrum(name=new_name, measurement=new_measurement, + # measure_type=self.measure_type) + # def __isub__(self, other): + # pass + # def __imul__(self, other): + # pass + # def __itruediv__(self, other): + # pass + # def __ifloordiv__(self, other): + # pass + # def __iiadd__(self, other): + # pass + # def __isub__(self, other): + # pass + # def __imul__(self, other): + # pass + # def __itruediv__(self, other): + # pass + # def __ifloordiv__(self, other): + # pass