#!/usr/bin/env python # # Flare mask from University of Utah L1 sparsity matched filter code (MAG1C) # BSD 3-Clause License # # Copyright (c) 2019, # Scientific Computing and Imaging Institute and # Utah Remote Sensing Applications Lab # All rights reserved. # # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are met: # # 1. Redistributions of source code must retain the above copyright notice, this # list of conditions and the following disclaimer. # # 2. Redistributions in binary form must reproduce the above copyright notice, # this list of conditions and the following disclaimer in the documentation # and/or other materials provided with the distribution. # # 3. Neither the name of the copyright holder nor the names of its # contributors may be used to endorse or promote products derived from # this software without specific prior written permission. # # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE # DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE # FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL # DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR # SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, # OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. # Author: Markus Foote (foote@sci.utah.edu) # Editor: Andrew Thorpe (Andrew.K.Thorpe@jpl.nasa.gov) import argparse import spectral #from spectral import * import numpy as np from skimage import morphology, measure from typing import Tuple, Optional import matplotlib.pyplot as plt # Define the default radiance value that will be used as a saturation threshold to identify flaring SAT_THRESH_DEFAULT = 6.0 # Define the version of script which will be written to ENVI .hdr file SCRIPT_VERSION='0.0.0' # Parser to permit command line options parser = argparse.ArgumentParser(description='Flare mask generated for AVIRIS-NG radiance files based on specified radiance threshold for a specified wavelength range.\n' 'f''v{SCRIPT_VERSION}', epilog='When using this software, please cite: \n' + ' TBD doi:xxxx.xxxx\n', formatter_class=argparse.RawDescriptionHelpFormatter, add_help=False, allow_abbrev=False) parser.add_argument('--txt', type=str, help='Text file and file path containing name of files to batch process.') parser.add_argument('--inpath', type=str, help='File path containing orthocorrected radiance files.') parser.add_argument('--outpath', type=str, help='File path to write outputs to.') parser.add_argument('-T', '--saturationthreshold', type=float, metavar='THRESHOLD', help='specify the threshold used for classifying pixels as saturated ' 'f''(default: {SAT_THRESH_DEFAULT})') parser.add_argument('-W', '--saturationwindow', type=float, nargs=2, metavar=('LOW', 'HIGH'), help='specify the contiguous wavelength window within which to detect saturation, independent ' 'of bands used in the filter (default: 1945, 2485 nanometers)') parser.add_argument('-M', '--maskgrowradius', type=str, metavar='RADIUS', nargs='?', const='150m', default=None, help='radius to use for expanding the saturation mask to cover (and exclude) flare-related ' 'anomalies. This value must include units: meters (abbreviated as m) or pixels ' '(abbreviated as px). If flag is given without a value, %(default)s will be used. This is ' 'a combined flag for enabling mask dilation and setting the distance to dilate.') parser.add_argument('-A', '--mingrowarea', type=int, metavar='PX_AREA', nargs='?', const=5, default=None, help='minimum number of pixels that must constitute a 2-connected saturation region for it to ' 'be grown by the mask-grow-radius value. If flag is provided without a value, ' '%(const)s pixels will be assumed as the value.') parser.add_argument('--saturation-processing-block-length', type=int, metavar='N', default=500, help='control the number of data lines pre-processed at once when using masking options') parser.add_argument('--visible-mask-growing-threshold', type=float, default=9.0, metavar='FLOAT', help='restrict mask dilation to only occur when 500 nm radiance is less than this value') parser.add_argument('-o', '--overwrite', action='store_true', help='Force the output files to overwrite any existing files. (default: %(default)s)') parser.add_argument('-h', '--help', action='help', help='show this help message and exit') args = parser.parse_args() print('Arguments:') print(args) # Text file path txt_path = args.txt # File path containing orthocorrected radiance files in_path = args.inpath # File path to write outputs to out_path = args.outpath # Read in text file of flights with open(txt_path, "r") as fd: files= fd.read().splitlines() # Go through each line of the text file for f in range(0,len(files)): #Go through each row of text file f_txt = str(files[f]) print('Processing flight',f_txt) # Open the specified radiance file as a memory-mapped object rdn_file = spectral.io.envi.open(in_path + f_txt + '.hdr') rdn_file_memmap = rdn_file.open_memmap(interleave='bip', writable=False) # Close memory-mapped object for radiance rdn_file_memmap.flush() def get_saturation_mask(data: np.ndarray, wave: np.ndarray, threshold: Optional[float] = None, waverange: Optional[Tuple[float, float]] = None) -> np.ndarray: """Calculates a mask of pixels that appear saturated (in the SWIR, by default). Pixels containing ANY radiance value above the provided threshold (default 6.0) within the wavelength window provided (default 1945 - 2485 nm). :param data: Radiance image to screen for sensor saturation. :param wave: vector of wavelengths (in nanometers) that correspond to the bands (last dimension) in the data. Caution: No input validation is performed, so this vector MUST be the same length as the data's last dimension. :param threshold: radiance value that defines the edge of saturation. :param waverange: wavelength range, defined as a tuple (low, high), to screen within for saturation. :return: Binary Mask with 1/True where saturation occurs, 0/False for normal pixels """ if threshold is None: threshold = SAT_THRESH_DEFAULT if waverange is None: waverange = (1945, 2485) is_saturated = (data[..., np.logical_and(wave >= waverange[0], wave <= waverange[1])] > threshold).any(axis=-1) return is_saturated def get_radius_in_pixels(value_str, metadata): if value_str.endswith('px'): return np.ceil(float(value_str.split('px')[0])) if value_str.endswith('m'): if 'map info' not in metadata: raise RuntimeError('Image does not have resolution specified. Try giving values in pixels.') if 'meters' not in metadata['map info'][10].lower(): raise RuntimeError('Unknown unit for image resolution.') meters_per_pixel_x = float(metadata['map info'][5]) meters_per_pixel_y = float(metadata['map info'][6]) if meters_per_pixel_x != meters_per_pixel_y: print('Warning: x and y resolutions are not equal, the average resolution will be used.') meters_per_pixel_x = (meters_per_pixel_y + meters_per_pixel_x) / 2.0 pixel_radius = float(value_str.split('m')[0]) / meters_per_pixel_x return np.ceil(pixel_radius) raise RuntimeError('Unknown unit specified.') # Get wavelengths from rdn file wavelengths = np.array(rdn_file.bands.centers) # If thresholding is enabled, calculate the mask and (if enabled) preprocess with dilation sat_mask_full = None print('Detecting saturated pixels...', end='') if args.maskgrowradius is not None: grow_radius_px = get_radius_in_pixels(args.maskgrowradius, rdn_file.metadata) selem = morphology.disk(radius=grow_radius_px, dtype=np.bool) idx_500 = np.argmin(np.absolute(wavelengths - 500)) sat_mask_full = np.zeros((rdn_file.nrows, rdn_file.ncols), dtype=np.uint8) block_overlap = np.ceil((args.mingrowarea if args.mingrowarea is not None else 0) + (grow_radius_px if args.maskgrowradius is not None else 0)).astype(np.int64) block_step = args.saturation_processing_block_length block_length = block_step + block_overlap line_idx_start_values = np.arange(start=0, stop=rdn_file.nrows, step=block_step) for line_block_start in line_idx_start_values: print('.', end='', flush=True) line_block_end = np.minimum(rdn_file.nrows, line_block_start + block_length) block_data = rdn_file.read_subregion((line_block_start, line_block_end), (0, rdn_file.ncols)) sat_mask_block = get_saturation_mask(data=block_data[:, :, :], wave=wavelengths, threshold=args.saturationthreshold, waverange=args.saturationwindow) if args.maskgrowradius is not None: sat_mask_grow_regions = np.zeros_like(sat_mask_block, dtype=np.uint8) for region in measure.regionprops(measure.label(sat_mask_block.astype(np.uint8), connectivity=2)): if args.mingrowarea is None or region.area >= args.mingrowarea: # Mark these large regions in the mask to get dilated for c in region.coords: sat_mask_grow_regions[c[0], c[1]] = 1 if block_data[c[0], c[1], idx_500] < args.visible_mask_growing_threshold else 0 sat_mask_large_grown = morphology.binary_dilation(image=sat_mask_grow_regions.astype(np.bool), selem=selem) sat_mask_out = sat_mask_large_grown.astype(np.uint8) sat_mask_out[sat_mask_block] = np.asarray(2, dtype=np.uint8) sat_mask_block = sat_mask_out sat_mask_full[line_block_start:line_block_end, ...][ np.logical_and(sat_mask_block == 1, sat_mask_full[line_block_start:line_block_end, ...] != 2)] = 1 sat_mask_full[line_block_start:line_block_end, ...][sat_mask_block == 2] = 2 #Specify output type output_dtype = np.uint8 # Create an image file for the output output_metadata = {'description': 'University of Utah Flare Mask Result.', 'band names': ['Saturation Mask (dimensionless'], 'interleave': 'bil', 'lines': rdn_file_memmap.shape[0], 'samples': rdn_file_memmap.shape[1], 'bands': 1, 'data type': spectral.io.envi.dtype_to_envi[np.dtype(output_dtype).char], 'algorithm settings': '{' f'version: {SCRIPT_VERSION}, ' + (f'saturationwindow: {args.saturationwindow if args.saturationwindow is not None else (1945, 2485)}, ' if args.maskgrowradius else '') + (f'saturationthreshold: {args.saturationthreshold if args.saturationthreshold is not None else SAT_THRESH_DEFAULT}, ' if args.maskgrowradius else '') + (f'buffer distance: {args.maskgrowradius if args.maskgrowradius is not None else 0}, ' if args.maskgrowradius else '') + (f'min contiguous px for buffer: {args.mingrowarea if args.mingrowarea is not None else 0}, ' if args.maskgrowradius else '') + (f'500 nm mask buffering threshold: {args.visible_mask_growing_threshold}, ' if args.maskgrowradius else '') + f'parsed cmdline args: {args}' '}'} # Save results output_metadata.update({'map info': rdn_file.metadata['map info']}) output_path = out_path + f_txt[:len('xxxYYYYMMDDtHHMMSS')] + '_flr_' + f_txt[len('xxxYYYYMMDDtHHMMSS_rdn_'):len('xxxYYYYMMDDtHHMMSS_rdn_v2x1')] + '_' + f_txt[len('xxxYYYYMMDDtHHMMSS_rdn_v2x1_'):len('xxxYYYYMMDDtHHMMSS_rdn_v2x1_img')] + '.hdr' spectral.envi.save_image(output_path, sat_mask_full, interleave='bil', ext='', metadata=output_metadata,force=args.overwrite) output_filename = f_txt[:len('xxxYYYYMMDDtHHMMSS')] + '_flr_' + f_txt[len('xxxYYYYMMDDtHHMMSS_rdn_'):len('xxxYYYYMMDDtHHMMSS_rdn_v2x1')] + '_' + f_txt[len('xxxYYYYMMDDtHHMMSS_rdn_v2x1_'):len('xxxYYYYMMDDtHHMMSS_rdn_v2x1_img')] print('Generated ' + output_filename) print('Completed all scenes')