diff --git a/micasense/capture.py b/micasense/capture.py index aa6254c0..86313056 100644 --- a/micasense/capture.py +++ b/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)): """ diff --git a/micasense/imageset.py b/micasense/imageset.py index f9f02401..6f9466fc 100644 --- a/micasense/imageset.py +++ b/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,8 +285,9 @@ 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, @@ -272,7 +295,6 @@ def process_imageset(self, } print('Processing {} Captures ...'.format(len(self.captures))) - # multiprocessing with concurrent futures if multiprocess: parallel_process(function=save_capture, iterable=self.captures, parameters=params,