From f1d6a25e5ae203e3382d640c21aa5b9af0071f02 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Hollstein?= Date: Wed, 6 Sep 2017 10:45:24 +0200 Subject: [PATCH] Implement save function for EnMAPL1Product_SensorGeo --- enpt/io/reader.py | 8 ++--- enpt/model/images.py | 65 ++++++++++++++++++++++++++++++------ enpt/model/metadata.py | 7 +++- enpt/utils/path_generator.py | 2 +- tests/test_l1b_reader.py | 34 +++++++++++++------ 5 files changed, 90 insertions(+), 26 deletions(-) diff --git a/enpt/io/reader.py b/enpt/io/reader.py index e3a6ce2..fe9e69c 100644 --- a/enpt/io/reader.py +++ b/enpt/io/reader.py @@ -36,24 +36,24 @@ class L1B_Reader(object): # read VNIR data # call L1_obj.vnir.arr.setter which sets L1_obj.swir.arr to an instance of GeoArray class - L1_obj.vnir.data = L1_obj.paths.vnir.imagedata + L1_obj.vnir.data = L1_obj.paths.vnir.data L1_obj.vnir.mask_clouds = L1_obj.paths.vnir.mask_clouds L1_obj.vnir.deadpixelmap = L1_obj.paths.vnir.deadpixelmap L1_obj.vnir.detector_meta = L1_obj.meta.vnir # read SWIR data # call L1_obj.swir.arr.setter which sets L1_obj.swir.arr to an instance of GeoArray class - L1_obj.swir.data = L1_obj.paths.swir.imagedata + L1_obj.swir.data = L1_obj.paths.swir.data L1_obj.swir.mask_clouds = L1_obj.paths.swir.mask_clouds L1_obj.swir.deadpixelmap = L1_obj.paths.swir.deadpixelmap L1_obj.swir.detector_meta = L1_obj.meta.swir # compute radiance and calculate snr L1_obj.DN2TOARadiance() - if snr_vnir is not None: + 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: + 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) diff --git a/enpt/model/images.py b/enpt/model/images.py index f287348..cf818c1 100644 --- a/enpt/model/images.py +++ b/enpt/model/images.py @@ -4,6 +4,9 @@ import logging from types import SimpleNamespace import numpy as np +from os import path, sep, makedirs +from shutil import copyfile +from xml.etree import ElementTree from geoarray import GeoArray, NoDataMask, CloudMask @@ -378,7 +381,7 @@ class EnMAP_Detector_SensorGeo(_EnMAP_Image): paths = SimpleNamespace() paths.root_dir = self._root_dir paths.metaxml = pathGen.get_path_metaxml() - paths.imagedata = pathGen.get_path_imagedata() + paths.data = pathGen.get_path_data() paths.mask_clouds = pathGen.get_path_cloudmask() paths.deadpixelmap = pathGen.get_path_deadpixelmap() paths.quicklook = pathGen.get_path_quicklook() @@ -394,15 +397,17 @@ class EnMAP_Detector_SensorGeo(_EnMAP_Image): :return: None """ # TODO move to processors.radiometric_transform? - if self.detector_meta.unitcode != 'DN': - raise RuntimeError("'DN2TOARadiance' is intended to convert digital numbers into TOA radiance. The current " - "unitcode must be 'DN'! Found %s." % self.detector_meta.unitcode) - - self.logger.info('Converting DN values to radiance for %s detector...' % self.detector_name) - self.data = (self.detector_meta.l_min + (self.detector_meta.l_max - self.detector_meta.l_min) / - (2 ** 16 - 1) * self.data[:]) - self.detector_meta.unit = "mW m^-2 sr^-1 nm^-1" - self.detector_meta.unitcode = "TOARad" + if self.detector_meta.unitcode == 'DN': + self.logger.info('Converting DN values to radiance for %s detector...' % self.detector_name) + self.data = (self.detector_meta.l_min + (self.detector_meta.l_max - self.detector_meta.l_min) / + (2 ** 16 - 1) * self.data[:]) + self.detector_meta.unit = "mW m^-2 sr^-1 nm^-1" + self.detector_meta.unitcode = "TOARad" + else: + self.logger.info( + "No is DN to Radiance conversion is performed since unitcode is not DN (found: {code}).".format( + code=self.detector_meta.unitcode) + ) class EnMAPL1Product_SensorGeo(object): @@ -458,6 +463,46 @@ class EnMAPL1Product_SensorGeo(object): self.vnir.DN2TOARadiance() self.swir.DN2TOARadiance() + def save(self, outdir: str, suffix="") -> str: + """Save this product to disk using almost the same format as for reading. + + :param outdir: Path to output directory + :return: Root path of written product + """ + product_dir = path.join( + path.abspath(outdir), + "{name}{suffix}".format( + name=[ff for ff in self.paths.root_dir.split(sep) if ff != ''][-1], + suffix=suffix) + ) + self.logger.info("Write product to: %s" % product_dir) + makedirs(product_dir, exist_ok=True) + + for detector_name in self.detector_attrNames: + detector = getattr(self, detector_name) + detector_paths = getattr(self.paths, detector_name) + + for atts, fmt in ((("deadpixelmap", "mask_clouds"), "GTIff"), + (("data",), "ENVI")): + for att in atts: + getattr(detector, att).save( + path.join(product_dir, path.basename(getattr(detector_paths, att))), fmt=fmt) + + copyfile( + src=detector_paths.quicklook, + dst=path.join(product_dir, path.basename(detector_paths.quicklook)) + ) + + xml = ElementTree.parse(self.paths.metaxml) + for xml_name, real_name in (("detector1", "vnir"), ("detector2", "swir")): + ele = xml.getroot().find(xml_name) + new_ele = ElementTree.Element("unitcode") + new_ele.text = getattr(self.meta, real_name).unitcode + ele.append(new_ele) + xml.write(path.join(product_dir, path.basename(self.paths.metaxml))) + + return product_dir + #################################### # EnPT EnMAP objects in map geometry diff --git a/enpt/model/metadata.py b/enpt/model/metadata.py index 60ae405..6290689 100644 --- a/enpt/model/metadata.py +++ b/enpt/model/metadata.py @@ -105,7 +105,12 @@ class EnMAP_Metadata_L1B_Detector_SensorGeo(object): 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()) - self.unitcode = 'DN' + + try: + self.unitcode = xml.findall("%s/unitcode" % lbl)[0].text + except IndexError: + self.unitcode = 'DN' + self.snr = None def calc_smile(self): diff --git a/enpt/utils/path_generator.py b/enpt/utils/path_generator.py index 045cb82..56e02c6 100644 --- a/enpt/utils/path_generator.py +++ b/enpt/utils/path_generator.py @@ -27,7 +27,7 @@ class PathGenL1BProduct(object): """Return the path of the metadata XML file.""" return glob(os.path.join(self.root_dir, "*_header.xml"))[0] - def get_path_imagedata(self): + def get_path_data(self): """Return the path of the image data file.""" return os.path.join(self.root_dir, self._find_in_metaxml("%s/filename" % self.detector_label)) diff --git a/tests/test_l1b_reader.py b/tests/test_l1b_reader.py index 4a45e63..991072e 100644 --- a/tests/test_l1b_reader.py +++ b/tests/test_l1b_reader.py @@ -11,6 +11,7 @@ Tests for `l1b_reader` module. import unittest from glob import glob import os +from os import path import tempfile import zipfile from datetime import datetime @@ -28,11 +29,14 @@ class Test_L1B_Reader(unittest.TestCase): def tearDown(self): pass - def test_read_inputdata(self): + def test_read_and_save_inputdata(self): from enpt.io.reader import L1B_Reader from enpt.model.images import EnMAPL1Product_SensorGeo print("Test reading EnMAP Level-1B products") + + rd = L1B_Reader(**self.user_config) + for l1b_file in self.pathList_testimages: with tempfile.TemporaryDirectory() as tmpdir: print("Tmp dir: %s" % tmpdir) @@ -41,18 +45,28 @@ class Test_L1B_Reader(unittest.TestCase): with zipfile.ZipFile(l1b_file, "r") as zf: zf.extractall(tmpdir) - + # without snr root_dir = os.path.join(tmpdir, os.path.basename(l1b_file).split(".zip")[0]) - L1_obj_no_snr = L1B_Reader(**self.user_config)\ - .read_inputdata(root_dir, observation_time=datetime(2015, 12, 7, 10)) + 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(tmpdir, "_no_snr")) + L1_obj = rd.read_inputdata(root_dir, observation_time=datetime(2015, 12, 7, 10)) + self.assertIsInstance(L1_obj, EnMAPL1Product_SensorGeo) - L1_obj_with_snr = L1B_Reader(**self.user_config) \ - .read_inputdata(root_dir, observation_time=datetime(2015, 12, 7, 10), - snr_vnir=os.path.join(tmpdir, "SNR_D1.hdr"), - snr_swir=os.path.join(tmpdir, "SNR_D2.hdr")) + # with snr + L1_obj = rd.read_inputdata( + root_dir, observation_time=datetime(2015, 12, 7, 10), + snr_vnir=os.path.join(tmpdir, "SNR_D1.hdr"), + snr_swir=os.path.join(tmpdir, "SNR_D2.hdr")) + self.assertIsInstance(L1_obj, EnMAPL1Product_SensorGeo) - self.assertIsInstance(L1_obj_no_snr, EnMAPL1Product_SensorGeo) - self.assertIsInstance(L1_obj_with_snr, EnMAPL1Product_SensorGeo) + root_dir = L1_obj.save(path.join(tmpdir), "with_snr") + L1_obj = rd.read_inputdata( + root_dir, observation_time=datetime(2015, 12, 7, 10), + snr_vnir=os.path.join(tmpdir, "SNR_D1.hdr"), + snr_swir=os.path.join(tmpdir, "SNR_D2.hdr")) + self.assertIsInstance(L1_obj, EnMAPL1Product_SensorGeo) if __name__ == "__main__": -- GitLab