From e703f3810740cab500fb2d67615351b70f91804f Mon Sep 17 00:00:00 2001 From: Daniel Scheffler Date: Fri, 6 Apr 2018 16:56:37 +0200 Subject: [PATCH 1/3] Revised SNR computation (not yet working). --- enpt/io/reader.py | 22 ++++++++----- enpt/model/images.py | 8 +++++ enpt/model/metadata.py | 70 ++++++++++++++++++++++------------------ tests/test_l1b_reader.py | 58 +++++++++++++++++++-------------- 4 files changed, 93 insertions(+), 65 deletions(-) diff --git a/enpt/io/reader.py b/enpt/io/reader.py index 9b7633e..f537d14 100644 --- a/enpt/io/reader.py +++ b/enpt/io/reader.py @@ -3,6 +3,8 @@ from datetime import datetime import logging +import tempfile +import zipfile import numpy as np from scipy.interpolate import interp1d @@ -20,7 +22,7 @@ class L1B_Reader(object): self.logger = logger or logging.getLogger(__name__) def read_inputdata(self, root_dir, observation_time: datetime, lon_lat_smpl: tuple=(15, 15), nsmile_coef: int=5, - snr_vnir: str=None, snr_swir: str=None): + compute_snr: bool=True): # TODO move to init? """Read L1B, DEM and spatial reference data. @@ -51,14 +53,18 @@ class L1B_Reader(object): L1_obj.swir.deadpixelmap = L1_obj.paths.swir.deadpixelmap L1_obj.swir.detector_meta = L1_obj.meta.swir - # compute radiance and calculate snr + # compute radiance L1_obj.DN2TOARadiance() - if snr_vnir is not None and L1_obj.meta.vnir.unit == 'mW m^-2 sr^-1 nm^-1': - self.logger.info("Compute SNR for vnir: %s" % snr_vnir) - L1_obj.vnir.detector_meta.calc_snr_vnir(detector=L1_obj.vnir, snr_data_fn=snr_vnir) - if snr_swir is not None and L1_obj.meta.vnir.unit == 'mW m^-2 sr^-1 nm^-1': - self.logger.info("Compute SNR for swir: %s" % snr_swir) - L1_obj.swir.detector_meta.calc_snr_swir(detector=L1_obj.swir, snr_data_fn=snr_swir) + + # compute SNR + if compute_snr: + with tempfile.TemporaryDirectory() as tmpDir, zipfile.ZipFile(self.cfg.path_l1b_snr_model, "r") as zf: + zf.extractall(tmpDir) + + if L1_obj.meta.vnir.unitcode == 'TOARad': + L1_obj.vnir.detector_meta.calc_snr_from_radiance(rad_data=L1_obj.vnir.data, dir_snr_models=tmpDir) + if L1_obj.meta.swir.unitcode == 'TOARad': + L1_obj.swir.detector_meta.calc_snr_from_radiance(rad_data=L1_obj.swir.data, dir_snr_models=tmpDir) return L1_obj diff --git a/enpt/model/images.py b/enpt/model/images.py index 009f5b4..332e6b2 100644 --- a/enpt/model/images.py +++ b/enpt/model/images.py @@ -500,9 +500,17 @@ class EnMAPL1Product_SensorGeo(object): xml = ElementTree.parse(self.paths.metaxml) for xml_name, real_name in (("detector1", "vnir"), ("detector2", "swir")): ele = xml.getroot().find(xml_name) + + # add unitcode new_ele = ElementTree.Element("unitcode") new_ele.text = getattr(self.meta, real_name).unitcode ele.append(new_ele) + + # add unit + new_ele = ElementTree.Element("unit") + new_ele.text = getattr(self.meta, real_name).unit + ele.append(new_ele) + xml.write(path.join(product_dir, path.basename(self.paths.metaxml))) return product_dir diff --git a/enpt/model/metadata.py b/enpt/model/metadata.py index f4cf41c..f16cc96 100644 --- a/enpt/model/metadata.py +++ b/enpt/model/metadata.py @@ -5,9 +5,11 @@ from datetime import datetime from xml.etree import ElementTree import logging import os +from typing import Union import numpy as np from scipy.interpolate import interp2d import spectral as sp +from geoarray import GeoArray from ..options.config import EnPTConfig from .srf import SRF @@ -112,12 +114,16 @@ class EnMAP_Metadata_L1B_Detector_SensorGeo(object): for corner in ("UL", "UR", "LL", "LR")] self.lats = self.interpolate_corners(*self.lat_UL_UR_LL_LR, *lon_lat_smpl) self.lons = self.interpolate_corners(*self.lon_UL_UR_LL_LR, *lon_lat_smpl) - self.unit = 'none' # '" ".join(xml.findall("%s/radiance_unit" % lbl)[0].text.split()) try: self.unitcode = xml.findall("%s/unitcode" % lbl)[0].text + # '" ".join(xml.findall("%s/radiance_unit" % lbl)[0].text.split()) + self.unit = xml.findall("%s/unit" % lbl)[0].text except IndexError: self.unitcode = 'DN' + self.unit = 'none' + except Exception: + raise self.snr = None @@ -135,39 +141,39 @@ class EnMAP_Metadata_L1B_Detector_SensorGeo(object): self.smile_coef # shape = (nwvl, nsmile_coef) ) # shape = (ncols, nwvl) - def calc_snr_vnir(self, detector, snr_data_fn: str) -> None: - """Compute EnMAP-VNIR SNR from radiance data. + def calc_snr_from_radiance(self, rad_data: Union[GeoArray, np.ndarray], dir_snr_models: str): + """Compute EnMAP SNR from radiance data for the given detector. - :param snr_data_fn: filename to SNR file - :param detector: vnir instance of EnMAP_Detector_SensorGeo + :param rad_data: image radiance data of EnMAP_Detector_SensorGeo + :param dir_snr_models: root directory where SNR model data is stored (must contain SNR_D1.bsq/SNR_D2.bsq) """ - assert self.unit == 'mW m^-2 sr^-1 nm^-1' - self.logger.info("Compute snr for: %s using: %s" % (self.detector_name, snr_data_fn)) - ds = sp.open_image(snr_data_fn) - p_hg = ds[0:3, :, :] - p_lg = ds[0:3, :, :] - l_th = np.squeeze(ds[6, :, :]) - - self.snr = np.zeros(detector.data.shape) - for irow in range(self.nrows): - hg = detector.data[irow, :, :] < l_th - l_hg = detector.data[irow, hg] - self.snr[irow, :, :][hg] = p_hg[0, hg] + p_hg[1, hg] * l_hg + p_hg[2, hg] * l_hg**2 - - lg = detector.data[irow, :, :] >= l_th - l_lg = detector.data[irow, lg] - self.snr[irow, :, :][lg] = p_lg[0, lg] + p_lg[1, lg] * l_lg + p_lg[2, lg] * l_lg ** 2 - - def calc_snr_swir(self, detector, snr_data_fn: str) -> None: - """Compute EnMAP SNR from radiance data. - - :param snr_data_fn: filename to SNR file - :param detector: swir instance of EnMAP_Detector_SensorGeo - """ - assert self.unit == 'mW m^-2 sr^-1 nm^-1' - self.logger.info("Compute snr for: %s using: %s" % (self.detector_name, snr_data_fn)) - pp = sp.open_image(snr_data_fn)[:, :, :] - self.snr = pp[0, :, :] + pp[1, :, :] * detector.data[:, :, :] + pp[2, :, :] * detector.data[:, :, :]**2 + path_snr_model = os.path.join(dir_snr_models, "SNR_D1.hdr" if self.detector_name == 'VNIR' else "SNR_D2.hdr") + + assert self.unitcode == 'TOARad' + self.logger.info("Computing SNR for: %s using: %s" % (self.detector_name, path_snr_model)) + + if self.detector_name == 'VNIR': + ds = sp.open_image(path_snr_model) + p_hg = ds[0:3, :, :] + p_lg = ds[0:3, :, :] + l_th = np.squeeze(ds[6, :, :]) + + self.snr = np.zeros(rad_data.shape) + for irow in range(self.nrows): + hg = rad_data[irow, :, :] < l_th + l_hg = rad_data[irow, hg] + try: + self.snr[irow, :, :][hg] = p_hg[0, hg] + p_hg[1, hg] * l_hg + p_hg[2, hg] * l_hg ** 2 + except: + a=1 + + lg = rad_data[irow, :, :] >= l_th + l_lg = rad_data[irow, lg] + self.snr[irow, :, :][lg] = p_lg[0, lg] + p_lg[1, lg] * l_lg + p_lg[2, lg] * l_lg ** 2 + + else: + pp = sp.open_image(path_snr_model)[:, :, :] + self.snr = pp[0, :, :] + pp[1, :, :] * rad_data[:, :, :] + pp[2, :, :] * rad_data[:, :, :] ** 2 @staticmethod def interpolate_corners(ul: float, ur: float, ll: float, lr: float, nx: int, ny: int): diff --git a/tests/test_l1b_reader.py b/tests/test_l1b_reader.py index c790d68..921aafe 100644 --- a/tests/test_l1b_reader.py +++ b/tests/test_l1b_reader.py @@ -42,33 +42,41 @@ class Test_L1B_Reader(unittest.TestCase): for l1b_file in self.pathList_testimages: print("Tmp dir: %s" % self.tmpdir) - with zipfile.ZipFile(self.config.path_l1b_snr_model, "r") as zf: + with zipfile.ZipFile(l1b_file, "r") as zf: zf.extractall(self.tmpdir) - with zipfile.ZipFile(l1b_file, "r") as zf: - zf.extractall(self.tmpdir) - # without snr - root_dir = os.path.join(self.tmpdir, os.path.basename(l1b_file).split(".zip")[0]) - L1_obj = rd.read_inputdata(root_dir, observation_time=datetime(2015, 12, 7, 10)) - self.assertIsInstance(L1_obj, EnMAPL1Product_SensorGeo) - - root_dir = L1_obj.save(path.join(self.tmpdir, "_no_snr")) - L1_obj = rd.read_inputdata(root_dir, observation_time=datetime(2015, 12, 7, 10)) - self.assertIsInstance(L1_obj, EnMAPL1Product_SensorGeo) - - # with snr - L1_obj = rd.read_inputdata( - root_dir, observation_time=datetime(2015, 12, 7, 10), - snr_vnir=os.path.join(self.tmpdir, "SNR_D1.hdr"), - snr_swir=os.path.join(self.tmpdir, "SNR_D2.hdr")) - self.assertIsInstance(L1_obj, EnMAPL1Product_SensorGeo) - - root_dir = L1_obj.save(path.join(self.tmpdir), "with_snr") - L1_obj = rd.read_inputdata( - root_dir, observation_time=datetime(2015, 12, 7, 10), - snr_vnir=os.path.join(self.tmpdir, "SNR_D1.hdr"), - snr_swir=os.path.join(self.tmpdir, "SNR_D2.hdr")) - self.assertIsInstance(L1_obj, EnMAPL1Product_SensorGeo) + root_dir = os.path.join(self.tmpdir, os.path.basename(l1b_file).split(".zip")[0]) + + ############### + # without snr # + ############### + + # read and write L1 data + L1_obj = rd.read_inputdata(root_dir, observation_time=datetime(2015, 12, 7, 10), compute_snr=False) + self.assertIsInstance(L1_obj, EnMAPL1Product_SensorGeo) + self.assertIsNone(L1_obj.vnir.detector_meta.snr) + self.assertIsNone(L1_obj.swir.detector_meta.snr) + root_dir_written_L1_data = L1_obj.save(path.join(self.tmpdir, "no_snr")) + + # read self written L1 data + L1_obj = rd.read_inputdata(root_dir_written_L1_data, observation_time=datetime(2015, 12, 7, 10), + compute_snr=False) + self.assertIsInstance(L1_obj, EnMAPL1Product_SensorGeo) + + ############ + # with snr # + ############ + + # read and write L1 data + L1_obj = rd.read_inputdata(root_dir, observation_time=datetime(2015, 12, 7, 10)) + self.assertIsInstance(L1_obj, EnMAPL1Product_SensorGeo) + self.assertIsNotNone(L1_obj.vnir.detector_meta.snr) + self.assertIsNotNone(L1_obj.swir.detector_meta.snr) + root_dir_written_L1_data = L1_obj.save(path.join(self.tmpdir, "with_snr")) + + # read self written L1 data + L1_obj = rd.read_inputdata(root_dir_written_L1_data, observation_time=datetime(2015, 12, 7, 10)) + self.assertIsInstance(L1_obj, EnMAPL1Product_SensorGeo) if __name__ == "__main__": -- GitLab From 4c7b6b696d92c3f83ae0a9f6fc7ab86ef7a8303b Mon Sep 17 00:00:00 2001 From: Daniel Scheffler Date: Mon, 9 Apr 2018 14:03:09 +0200 Subject: [PATCH 2/3] Removed deprecated equation. --- enpt/model/metadata.py | 1 - 1 file changed, 1 deletion(-) diff --git a/enpt/model/metadata.py b/enpt/model/metadata.py index f16cc96..f3b67c9 100644 --- a/enpt/model/metadata.py +++ b/enpt/model/metadata.py @@ -244,7 +244,6 @@ class EnMAP_Metadata_L1B_SensorGeo(object): """ # FIXME observation time is currently missing in the XML self.observation_datetime = observation_time - # self.earthSunDist = np.cos(np.deg2rad(self.geom_sun_zenith)) # this seems to be wrong -> Andre? self.earthSunDist = self.get_earth_sun_distance(self.observation_datetime) def get_earth_sun_distance(self, acqDate: datetime): -- GitLab From da448a8539e5f914f14b94dd970a98d569ea4508 Mon Sep 17 00:00:00 2001 From: Daniel Scheffler Date: Mon, 9 Apr 2018 23:52:05 +0200 Subject: [PATCH 3/3] Fixed calc_snr_from_radiance(). Updated minimal version of geoarray. --- enpt/model/metadata.py | 49 +++++++++++++------ requirements.txt | 2 +- setup.py | 2 +- .../context/environment_enpt.yml | 2 +- 4 files changed, 36 insertions(+), 19 deletions(-) diff --git a/enpt/model/metadata.py b/enpt/model/metadata.py index f3b67c9..c78a24a 100644 --- a/enpt/model/metadata.py +++ b/enpt/model/metadata.py @@ -144,36 +144,53 @@ class EnMAP_Metadata_L1B_Detector_SensorGeo(object): def calc_snr_from_radiance(self, rad_data: Union[GeoArray, np.ndarray], dir_snr_models: str): """Compute EnMAP SNR from radiance data for the given detector. + SNR equation: SNR = p0 + p1 * LTOA + p2 * LTOA ^ 2 [W/(m^2 sr nm)]. + NOTE: The SNR model files (SNR_D1.bsq/SNR_D2.bsq) contain polynomial coefficients needed to compute SNR. + SNR_D1.bsq: SNR model for EnMAP SWIR detector (contains high gain and low gain model coefficients) + - 1000 columns (for 1000 EnMAP columns) + - 88 bands for 88 EnMAP VNIR bands + - 7 lines: - 3 lines: high gain coefficients + - 3 lines: low gain coefficients + - 1 line: threshold needed to decide about high gain or low gain + SNR_D2.bsq: SNR model for EnMAP SWIR detector + - 1000 columns (for 1000 EnMAP columns) + - x bands for x EnMAP SWIR bands + - 3 lines for 3 coefficients + :param rad_data: image radiance data of EnMAP_Detector_SensorGeo :param dir_snr_models: root directory where SNR model data is stored (must contain SNR_D1.bsq/SNR_D2.bsq) """ path_snr_model = os.path.join(dir_snr_models, "SNR_D1.hdr" if self.detector_name == 'VNIR' else "SNR_D2.hdr") + rad_data = np.array(rad_data) assert self.unitcode == 'TOARad' self.logger.info("Computing SNR for: %s using: %s" % (self.detector_name, path_snr_model)) if self.detector_name == 'VNIR': - ds = sp.open_image(path_snr_model) - p_hg = ds[0:3, :, :] - p_lg = ds[0:3, :, :] - l_th = np.squeeze(ds[6, :, :]) + gA = sp.open_image(path_snr_model) + coeffs_highgain = gA[0:3, :, :] + coeffs_lowgain = gA[3:6, :, :] + gain_threshold = np.squeeze(gA[6, :, :]) self.snr = np.zeros(rad_data.shape) for irow in range(self.nrows): - hg = rad_data[irow, :, :] < l_th - l_hg = rad_data[irow, hg] - try: - self.snr[irow, :, :][hg] = p_hg[0, hg] + p_hg[1, hg] * l_hg + p_hg[2, hg] * l_hg ** 2 - except: - a=1 - - lg = rad_data[irow, :, :] >= l_th - l_lg = rad_data[irow, lg] - self.snr[irow, :, :][lg] = p_lg[0, lg] + p_lg[1, lg] * l_lg + p_lg[2, lg] * l_lg ** 2 + highgain_mask = rad_data[irow, :, :] < gain_threshold # a single row + rad_highgain = rad_data[irow, highgain_mask] + self.snr[irow, :, :][highgain_mask] = \ + coeffs_highgain[0, highgain_mask] + \ + coeffs_highgain[1, highgain_mask] * rad_highgain + \ + coeffs_highgain[2, highgain_mask] * rad_highgain ** 2 + + lowgain_mask = rad_data[irow, :, :] >= gain_threshold + rad_lowgain = rad_data[irow, lowgain_mask] + self.snr[irow, :, :][lowgain_mask] = \ + coeffs_lowgain[0, lowgain_mask] + \ + coeffs_lowgain[1, lowgain_mask] * rad_lowgain + \ + coeffs_lowgain[2, lowgain_mask] * rad_lowgain ** 2 else: - pp = sp.open_image(path_snr_model)[:, :, :] - self.snr = pp[0, :, :] + pp[1, :, :] * rad_data[:, :, :] + pp[2, :, :] * rad_data[:, :, :] ** 2 + coeffs = sp.open_image(path_snr_model)[:, :, :] + self.snr = coeffs[0, :, :] + coeffs[1, :, :] * rad_data[:, :, :] + coeffs[2, :, :] * rad_data[:, :, :] ** 2 @staticmethod def interpolate_corners(ul: float, ur: float, ll: float, lr: float, nx: int, ny: int): diff --git a/requirements.txt b/requirements.txt index 0b5f062..74df5e2 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,6 +1,6 @@ numpy scipy -geoarray>=0.6.12 +geoarray>=0.7.15 spectral>=0.16 cerberus jsmin diff --git a/setup.py b/setup.py index 12b89a8..685e018 100644 --- a/setup.py +++ b/setup.py @@ -14,7 +14,7 @@ with open("enpt/version.py") as version_file: exec(version_file.read(), version) requirements = [ # put package requirements here - 'numpy', 'scipy', 'geoarray>=0.6.12', 'spectral>=0.16', 'cerberus', 'jsmin', 'matplotlib' + 'numpy', 'scipy', 'geoarray>=0.7.15', 'spectral>=0.16', 'cerberus', 'jsmin', 'matplotlib' ] test_requirements = ['coverage', 'nose', 'nose-htmloutput', 'rednose'] diff --git a/tests/gitlab_CI_docker/context/environment_enpt.yml b/tests/gitlab_CI_docker/context/environment_enpt.yml index f82fd68..f1b7bbf 100644 --- a/tests/gitlab_CI_docker/context/environment_enpt.yml +++ b/tests/gitlab_CI_docker/context/environment_enpt.yml @@ -24,7 +24,7 @@ dependencies: - bokeh - pip: - scipy - - geoarray>=0.6.12 + - geoarray>=0.7.15 - spectral>=0.16 - cerberus - jsmin -- GitLab