Reconcile save capture 110

micasense/capture.py

@@ -29,6 +29,7 @@
 
 import math
 import os
+import warnings
 
 import cv2
 import imageio
@@ -38,6 +39,8 @@
 import micasense.imageutils as imageutils
 import micasense.plotutils as plotutils
 
+warnings.simplefilter(action="once")
+
 
 class Capture(object):
     """
@@ -80,6 +83,7 @@ def __init__(self, images, panel_corners=None):
         else:
             self.panel_corners = panel_corners
 
+        self.__image_type = None
         self.__aligned_capture = None
 
     def set_panel_corners(self, panel_corners):
@@ -201,6 +205,7 @@ def clear_image_data(self):
         for img in self.images:
             img.clear_image_data()
         self.__aligned_capture = None
+        self.__image_type = None
 
     def center_wavelengths(self):
         """Returns a list of the image center wavelengths in nanometers."""
@@ -468,26 +473,39 @@ def get_warp_matrices(self, ref_index=None):
         warp_matrices = [np.linalg.inv(im.get_homography(ref)) for im in self.images]
         return [w / w[2, 2] for w in warp_matrices]
 
-    def create_aligned_capture(self, irradiance_list=None, warp_matrices=None, normalize=False, img_type=None,
+    def create_aligned_capture(self, img_type, irradiance_list=None, warp_matrices=None, normalize=False,
                                motion_type=cv2.MOTION_HOMOGRAPHY):
         """
         Creates aligned Capture. Computes undistorted radiance or reflectance images if necessary.
+        :param img_type: str 'radiance' or 'reflectance' depending on image metadata.
         :param irradiance_list: List of mean panel region irradiance.
         :param warp_matrices: 2d List of warp matrices derived from Capture.get_warp_matrices()
         :param normalize: FIXME: This parameter isn't used?
-        :param img_type: str 'radiance' or 'reflectance' depending on image metadata.
         :param motion_type: OpenCV import. Also know as warp_mode. MOTION_HOMOGRAPHY or MOTION_AFFINE.
                             For Altum images only use HOMOGRAPHY.
         :return: ndarray with alignment changes
         """
-        if img_type is None and irradiance_list is None and self.dls_irradiance() is None:
+        if img_type == 'radiance':
             self.compute_undistorted_radiance()
-            img_type = 'radiance'
-        elif img_type is None:
-            if irradiance_list is None:
+            self.__image_type = 'radiance'
+        elif img_type == 'reflectance':
+            # TODO: Handle pre-flight Panel cap + post-flight Panel cap + DLS values
+            # Add use_dls or similar user option to properly scale irradiance_list in different configurations.
+            # Alternatively if this should be done elsewhere then dls_irradiance() should probably not be used here.
+
+            # if no irradiance values provided, attempt to use DLS
+            if irradiance_list is None and self.dls_present():
                 irradiance_list = self.dls_irradiance() + [0]
+            elif irradiance_list is None and not self.dls_present():
+                raise RuntimeError('Reflectance output requested, but no irradiance values given and no DLS values '
+                                   'found in image metadata.')
+
             self.compute_undistorted_reflectance(irradiance_list)
-            img_type = 'reflectance'
+            self.__image_type = 'reflectance'
+        else:
+            raise RuntimeError('Unknown img_type output requested: {}\nMust be "radiance" or "reflectance".'
+                               .format(img_type))
+
         if warp_matrices is None:
             warp_matrices = self.get_warp_matrices()
         cropped_dimensions, _ = imageutils.find_crop_bounds(self, warp_matrices, warp_mode=motion_type)
@@ -508,25 +526,46 @@ def aligned_shape(self):
             raise RuntimeError("Call Capture.create_aligned_capture() prior to saving as stack.")
         return self.__aligned_capture.shape
 
-    def save_capture_as_stack(self, out_file_name, sort_by_wavelength=False, photometric='MINISBLACK'):
+    def save_capture_as_stack(self, out_file_name, out_data_type='GDT_UInt16',
+                              sort_by_wavelength=False, photometric='MINISBLACK'):
         """
         Output the Images in the Capture object as GTiff image stack.
         :param out_file_name: str system file path
+        :param out_data_type: str GDT_Float32 or GDT_UInt16
+            Default: GDT_UInt16 will write images as scaled reflectance values. EO (32769=100%)
+            and LWIR in centi-Kelvin (0-65535).
+            GDT_Float32 will write images as floating point reflectance values. EO (1.0=100%)
+            and LWIR in floating point Celsius.
+            https://gdal.org/api/raster_c_api.html#_CPPv412GDALDataType
         :param sort_by_wavelength: boolean
         :param photometric: str GDAL argument for GTiff color matching
         """
-        from osgeo.gdal import GetDriverByName, GDT_UInt16
+        from osgeo.gdal import GetDriverByName, GDT_UInt16, GDT_Float32
         if self.__aligned_capture is None:
             raise RuntimeError("Call Capture.create_aligned_capture() prior to saving as stack.")
 
         rows, cols, bands = self.__aligned_capture.shape
         driver = GetDriverByName('GTiff')
 
-        out_raster = driver.Create(out_file_name, cols, rows, bands, GDT_UInt16,
+        # force correct output datatype
+        if self.__image_type == 'radiance':
+            gdal_type = GDT_Float32  # force floating point values for radiance and degrees C
+        elif self.__image_type == 'reflectance' and out_data_type == 'GDT_UInt16':
+            gdal_type = GDT_UInt16
+        elif self.__image_type == 'reflectance' and out_data_type == 'GDT_Float32':
+            gdal_type = GDT_Float32
+        else:
+            warnings.warn(message='Output data type in Capture.save_capture_as_bands() was called as {}. '
+                                  'Must use "GDT_UInt16" or "GDT_Float32". Defaulting to GDT_UInt16...'
+                          .format(out_data_type),
+                          category=UserWarning)
+            gdal_type = GDT_UInt16
+
+        out_raster = driver.Create(out_file_name, cols, rows, bands, gdal_type,
                                    options=['INTERLEAVE=BAND', 'COMPRESS=DEFLATE', f'PHOTOMETRIC={photometric}'])
         try:
             if out_raster is None:
-                raise IOError("could not load gdal GeoTiff driver")
+                raise IOError("Could not load GDAL GeoTiff driver.")
 
             if sort_by_wavelength:
                 eo_list = list(np.argsort(np.array(self.center_wavelengths())[self.eo_indices()]))
@@ -538,20 +577,87 @@ def save_capture_as_stack(self, out_file_name, sort_by_wavelength=False, photome
                 out_data = self.__aligned_capture[:, :, in_band]
                 out_data[out_data < 0] = 0
                 out_data[out_data > 2] = 2  # limit reflectance data to 200% to allow some specular reflections
-                out_band.WriteArray(out_data * 32768)  # scale reflectance images so 100% = 32768
+                # scale reflectance images so 100% = 32768
+                out_band.WriteArray(out_data * 32768 if gdal_type == 2 else out_data)
                 out_band.FlushCache()
 
             for out_band, in_band in enumerate(self.lw_indices()):
                 out_band = out_raster.GetRasterBand(len(eo_list) + out_band + 1)
                 # scale data from float degC to back to centi-Kelvin to fit into uint16
-                out_data = (self.__aligned_capture[:, :, in_band] + 273.15) * 100
-                out_data[out_data < 0] = 0
-                out_data[out_data > 65535] = 65535
+                out_data = (self.__aligned_capture[:, :, in_band] + 273.15) * 100 if gdal_type == 2 \
+                    else self.__aligned_capture[:, :, in_band]
+                if gdal_type == 2:
+                    out_data[out_data < 0] = 0
+                    out_data[out_data > 65535] = 65535
                 out_band.WriteArray(out_data)
                 out_band.FlushCache()
         finally:
             out_raster = None
 
+    def save_capture_as_bands(self, out_file_name, out_data_type='GDT_UInt16', photometric='MINISBLACK'):
+        """
+        Output the Images in the Capture object as separate GTiffs.
+        :param out_file_name: str system file path without file extension
+        :param out_data_type: str GDT_Float32 or GDT_UInt16
+            Default: GDT_UInt16 will write images as scaled reflectance values. EO (32769=100%)
+            and LWIR in centi-Kelvin (0-65535).
+            GDT_Float32 will write images as floating point reflectance values. EO (1.0=100%)
+            and LWIR in floating point Celsius.
+            https://gdal.org/api/raster_c_api.html#_CPPv412GDALDataType
+        :param photometric: str GDAL argument for GTiff color matching
+        :return: None
+        """
+        from osgeo.gdal import GetDriverByName, GDT_UInt16, GDT_Float32
+        if self.__aligned_capture is None:
+            raise RuntimeError("Call Capture.create_aligned_capture() prior to saving as stack.")
+
+        # force correct output datatype
+        if self.__image_type == 'radiance':
+            gdal_type = GDT_Float32  # force floating point values for radiance and degrees C
+        elif self.__image_type == 'reflectance' and out_data_type == 'GDT_UInt16':
+            gdal_type = GDT_UInt16
+        elif self.__image_type == 'reflectance' and out_data_type == 'GDT_Float32':
+            gdal_type = GDT_Float32
+        else:
+            warnings.warn(message='Output data type in Capture.save_capture_as_bands() was called as {}. '
+                                  'Must use "GDT_UInt16" or "GDT_Float32". Defaulting to GDT_UInt16...'
+                          .format(out_data_type),
+                          category=UserWarning)
+            gdal_type = GDT_UInt16
+
+        # predictably handle accidental .tif.tif values
+        if out_file_name.endswith('.tif'):
+            out_file_path = out_file_name[:-4]
+        else:
+            out_file_path = out_file_name
+
+        rows, cols, bands = self.__aligned_capture.shape
+        driver = GetDriverByName('GTiff')
+
+        for i in self.eo_indices():
+            out_raster = driver.Create(out_file_path + f'_{i + 1}.tif', cols, rows, 1, gdal_type,
+                                       options=['INTERLEAVE=BAND', 'COMPRESS=DEFLATE', f'PHOTOMETRIC={photometric}'])
+            out_band = out_raster.GetRasterBand(1)
+            out_data = self.__aligned_capture[:, :, i]
+            out_data[out_data < 0] = 0
+            out_data[out_data > 2] = 2  # limit reflectance data to 200% to allow some specular reflections
+            # if GDT_UInt16, scale reflectance images so 100% = 32768. GDT_UInt16 resolves to 2.
+            out_band.WriteArray(out_data * 32768 if gdal_type == 2 else out_data)
+            out_band.FlushCache()
+
+        for i in self.lw_indices():
+            out_raster = driver.Create(out_file_path + f'_{i + 1}.tif', cols, rows, 1, gdal_type,
+                                       options=['INTERLEAVE=BAND', 'COMPRESS=DEFLATE', f'PHOTOMETRIC={photometric}'])
+            out_band = out_raster.GetRasterBand(1)
+            # if GDT_UInt16, scale data from float degC to back to centi-Kelvin to fit into UInt16.
+            out_data = (self.__aligned_capture[:, :, i] + 273.15) * 100 if gdal_type == 2 \
+                else self.__aligned_capture[:, :, i]
+            if gdal_type == 2:
+                out_data[out_data < 0] = 0
+                out_data[out_data > 65535] = 65535
+            out_band.WriteArray(out_data)
+            out_band.FlushCache()
+
     def save_capture_as_rgb(self, out_file_name, gamma=1.4, downsample=1, white_balance='norm', hist_min_percent=0.5,
                             hist_max_percent=99.5, sharpen=True, rgb_band_indices=(2, 1, 0)):
         """

micasense/imageset.py

@@ -82,27 +82,26 @@ def parallel_process(function, iterable, parameters, progress_callback=None, use
 
 def save_capture(params, cap):
     """
-    Process an ImageSet according to program parameters. Saves rgb
-    :param params: dict of program parameters from ImageSet.process_imageset()
+    Convenience function to save Captures in an ImageSet according to input parameters.
+    :param params: dict of input parameters from ImageSet.process_imageset()
     :param cap: micasense.capture.Capture object
     """
     try:
         # align capture
         if len(cap.images) == params['capture_len']:
             cap.create_aligned_capture(
-                irradiance_list=params['irradiance'],
-                warp_matrices=params['warp_matrices'],
-                img_type=params['img_type']
+                img_type=params['img_type'],
+                irradiance_list=params['irradiance_list'],
+                warp_matrices=params['warp_matrices']
             )
         else:
             print(f"\tCapture {cap.uuid} only has {len(cap.images)} Images. Should have {params['capture_len']}. "
                   f"Skipping...")
             return
-
         if params['output_stack_dir']:
             output_stack_file_path = os.path.join(params['output_stack_dir'], cap.uuid + '.tif')
             if params['overwrite'] or not os.path.exists(output_stack_file_path):
-                cap.save_capture_as_stack(output_stack_file_path)
+                cap.save_capture_as_stack(output_stack_file_path, out_data_type=params['out_data_type'])
         if params['output_rgb_dir']:
             output_rgb_file_path = os.path.join(params['output_rgb_dir'], cap.uuid + '.jpg')
             if params['overwrite'] or not os.path.exists(output_rgb_file_path):
@@ -224,28 +223,40 @@ def dls_irradiance(self):
         return series
 
     def process_imageset(self,
+                         img_type,
                          output_stack_directory=None,
                          output_rgb_directory=None,
                          warp_matrices=None,
-                         irradiance=None,
-                         img_type=None,
+                         irradiance_list=None,
+                         out_data_type='GDT_UInt16',
                          multiprocess=True,
                          overwrite=False,
                          progress_callback=None,
                          use_tqdm=False):
         """
         Write band stacks and rgb thumbnails to disk.
+        :param img_type: str 'radiance' or 'reflectance'. Desired image output type.
+            - A 'radiance' value will output EO Images as radiance floating point values,
+             and LWIR Images in floating point Celsius.
+            - A 'reflectance' value will output EO Images as reflectance  UInt16
+             and LWIR Images as UInt16 centi-Kelvin using the provided irradiance_list.
         :param warp_matrices: 2d List of warp matrices derived from Capture.get_warp_matrices()
         :param output_stack_directory: str system file path to output stack directory
         :param output_rgb_directory: str system file path to output thumbnail directory
-        :param irradiance: List returned from Capture.dls_irradiance() or Capture.panel_irradiance()    <-- TODO: Write a better docstring for this
-        :param img_type: str 'radiance' or 'reflectance'. Desired image output type.
+        :param irradiance_list: List returned from Capture.dls_irradiance() or Capture.panel_irradiance()    <-- TODO: Write a better docstring for this
+        :param out_data_type: str GDT_Float32 or GDT_UInt16
+            Default: GDT_UInt16 will write images as scaled reflectance values. EO (32769=100%)
+            and LWIR in centi-Kelvin (0-65535).
+            GDT_Float32 will write images as floating point reflectance values. EO (1.0=100%)
+            and LWIR in floating point Celsius.
+            https://gdal.org/api/raster_c_api.html#_CPPv412GDALDataType
         :param multiprocess: boolean True to use multiprocessing module
         :param overwrite: boolean True to overwrite existing files
         :param progress_callback: function to report progress to
         :param use_tqdm: boolean True to use tqdm progress bar
         """
 
+        # progress_callback deprecation warning
         if progress_callback is not None:
             warnings.warn(message='The progress_callback parameter will be deprecated in favor of use_tqdm',
                           category=PendingDeprecationWarning)
@@ -254,6 +265,17 @@ def process_imageset(self,
         if output_stack_directory is None and output_rgb_directory is None:
             raise RuntimeError('No output requested for the ImageSet.')
 
+        # ensure output type makes sense
+        if img_type not in ['radiance', 'reflectance']:
+            raise RuntimeError('Unknown img_type output requested: {}\nMust be "radiance" or "reflectance".'
+                               .format(img_type))
+
+        # warn user if only DLS metadata used in reflectance processing.
+        if img_type == 'reflectance' and irradiance_list is None:
+            warnings.warn(message='Reflectance output requested, but no irradiance_list value given. '
+                                  'Attempting to process reflectance from DLS metadata...',
+                          category=UserWarning)
+
         # make output dirs if not exist
         if output_stack_directory is not None and not os.path.exists(output_stack_directory):
             os.mkdir(output_stack_directory)
@@ -263,16 +285,16 @@ def process_imageset(self,
         # processing parameters
         params = {
             'warp_matrices': warp_matrices,
-            'irradiance': irradiance,
+            'irradiance_list': irradiance_list,
             'img_type': img_type,
+            'out_data_type': out_data_type,
             'capture_len': len(self.captures[0].images),
             'output_stack_dir': output_stack_directory,
             'output_rgb_dir': output_rgb_directory,
             'overwrite': overwrite,
         }
 
         print('Processing {} Captures ...'.format(len(self.captures)))
-
         # multiprocessing with concurrent futures
         if multiprocess:
             parallel_process(function=save_capture, iterable=self.captures, parameters=params,