diff --git a/etc/composites/viirs.cfg b/etc/composites/viirs.cfg index d6dc91d7..cbfd0ff3 100644 --- a/etc/composites/viirs.cfg +++ b/etc/composites/viirs.cfg @@ -9,7 +9,7 @@ sensor=viirs name=histogram_dnb standard_name=equalized_radiance compositor=mpop.composites.viirs.HistogramDNB -prerequisites=DNB,I_SZA +prerequisites=DNB,DNB_SZA sensor=viirs units=1 @@ -20,7 +20,15 @@ compositor=mpop.composites.viirs.AdaptiveDNB adaptive_day=multiple adaptive_mixed=always adaptive_night=never -prerequisites=DNB,I_SZA +prerequisites=DNB,DNB_SZA +sensor=viirs +units=1 + +[composite:dynamic_dnb] +name=dynamic_dnb +standard_name=equalized_radiance +compositor=mpop.composites.viirs.ERFDNB +prerequisites=DNB,DNB_SZA,DNB_LZA sensor=viirs units=1 diff --git a/etc/readers/viirs_sdr.cfg b/etc/readers/viirs_sdr.cfg index 69267f03..dac0d92e 100644 --- a/etc/readers/viirs_sdr.cfg +++ b/etc/readers/viirs_sdr.cfg @@ -256,14 +256,22 @@ navigation=gdnbo file_type=svdnb file_key=radiance_dnb -[dataset:i_sza] -name=I_SZA +[dataset:dnb_sza] +name=DNB_SZA standard_name=solar_zenith_angle resolution=743 -navigation=gitco +navigation=gdnbo file_type=gdnbo file_key=solar_zenith_angle +[dataset:dnb_lza] +name=DNB_LZA +standard_name=lunar_zenith_angle +resolution=743 +navigation=gdnbo +file_type=gdnbo +file_key=lunar_zenith_angle + [metadata:moon_illumination_fraction] name=moon_illumination_fraction file_type=NAVIGATION @@ -451,6 +459,9 @@ variable_name=All_Data/{file_group}_All/BrightnessTemperatureFactors [file_key:solar_zenith_angle] variable_name=All_Data/{file_group}_All/SolarZenithAngle +[file_key:lunar_zenith_angle] +variable_name=All_Data/{file_group}_All/LunarZenithAngle + [file_key:beginning_date] variable_name=Data_Products/{file_group}/{file_group}_Aggr/attr/AggregateBeginningDate diff --git a/mpop/composites/__init__.py b/mpop/composites/__init__.py index 9098e392..eec838cf 100644 --- a/mpop/composites/__init__.py +++ b/mpop/composites/__init__.py @@ -35,23 +35,13 @@ class CompositeBase(InfoObject): - - def __init__(self, name, compositor, prerequisites, default_image_config=None, - **kwargs): + def __init__(self, name, prerequisites=[], optional_prerequisites=[], metadata_requirements=[], **kwargs): # Required info kwargs["name"] = name - kwargs["compositor"] = compositor - kwargs["prerequisites"] = [] - for prerequisite in prerequisites.split(","): - try: - kwargs["prerequisites"].append(float(prerequisite)) - except ValueError: - kwargs["prerequisites"].append(prerequisite) - InfoObject.__init__(self, **kwargs) - if default_image_config is None: - return - for key, value in default_image_config.iteritems(): - self.info.setdefault(key, value) + kwargs["prerequisites"] = prerequisites + kwargs["optional_prerequisites"] = optional_prerequisites + kwargs["metadata_requirements"] = metadata_requirements + super(CompositeBase, self).__init__(**kwargs) @property def prerequisites(self): diff --git a/mpop/composites/viirs.py b/mpop/composites/viirs.py index 431434b6..512a035a 100644 --- a/mpop/composites/viirs.py +++ b/mpop/composites/viirs.py @@ -28,6 +28,7 @@ from mpop.composites import CompositeBase from mpop.projectable import Projectable import numpy as np +from scipy.special import erf LOG = logging.getLogger(__name__) @@ -141,6 +142,15 @@ def __call__(self, datasets, **info): class AdaptiveDNB(HistogramDNB): + """Adaptive histogram equalized DNB composite. + + The logic for this code was taken from Polar2Grid and was originally developed by Eva Schiffer (SSEC). + + This composite separates the DNB data in to 3 main regions: Day, Night, and Mixed. Each region is + equalized separately to bring out the most information from the region due to the high dynamic range + of the DNB data. Optionally, the mixed region can be separated in to multiple smaller regions by + using the `mixed_degree_step` keyword. + """ def __init__(self, *args, **kwargs): """Initialize the compositor with values from the user or from the configuration file. @@ -235,6 +245,90 @@ def __call__(self, datasets, **info): return output_dataset +class ERFDNB(CompositeBase): + """Equalized DNB composite using the error function (erf). + + The logic for this code was taken from Polar2Grid and was originally developed by Curtis Seaman and Steve Miller. + The original code was written in IDL and is included as comments in the code below. + """ + def __init__(self, *args, **kwargs): + self.saturation_correction = kwargs.pop("saturation_correction", False) in [True, "True", "true"] + kwargs.setdefault("metadata_requirements", ["moon_illumination_fraction"]) + super(ERFDNB, self).__init__(*args, **kwargs) + + def __call__(self, datasets, **info): + if len(datasets) != 3: + raise ValueError("Expected 3 datasets, got %d" % (len(datasets),)) + + dnb_data = datasets[0] + sza_data = datasets[1] + lza_data = datasets[2] + good_mask = ~(dnb_data.mask | sza_data.mask) + output_dataset = dnb_data.copy() + output_dataset.mask = ~good_mask + + moon_illum_name = self.info["metadata_requirements"][0] + if moon_illum_name not in dnb_data.info: + raise RuntimeError("Could not find '%s' metadata in DNB dataset" % (moon_illum_name,)) + # convert to decimal instead of % + moon_illum_fraction = np.mean(dnb_data.info[moon_illum_name]) * 0.01 + + ### From Steve Miller and Curtis Seaman + # maxval = 10.^(-1.7 - (((2.65+moon_factor1+moon_factor2))*(1+erf((solar_zenith-95.)/(5.*sqrt(2.0)))))) + # minval = 10.^(-4. - ((2.95+moon_factor2)*(1+erf((solar_zenith-95.)/(5.*sqrt(2.0)))))) + # scaled_radiance = (radiance - minval) / (maxval - minval) + # radiance = sqrt(scaled_radiance) + + ### Version 2: Update from Curtis Seaman + # maxval = 10.^(-1.7 - (((2.65+moon_factor1+moon_factor2))*(1+erf((solar_zenith-95.)/(5.*sqrt(2.0)))))) + # minval = 10.^(-4. - ((2.95+moon_factor2)*(1+erf((solar_zenith-95.)/(5.*sqrt(2.0)))))) + # saturated_pixels = where(radiance gt maxval, nsatpx) + # saturation_pct = float(nsatpx)/float(n_elements(radiance)) + # print, 'Saturation (%) = ', saturation_pct + # + # while saturation_pct gt 0.005 do begin + # maxval = maxval*1.1 + # saturated_pixels = where(radiance gt maxval, nsatpx) + # saturation_pct = float(nsatpx)/float(n_elements(radiance)) + # print, saturation_pct + # endwhile + # + # scaled_radiance = (radiance - minval) / (maxval - minval) + # radiance = sqrt(scaled_radiance) + + moon_factor1 = 0.7 * (1.0 - moon_illum_fraction) + moon_factor2 = 0.0022 * lza_data + erf_portion = 1 + erf((sza_data - 95.0) / (5.0 * np.sqrt(2.0))) + max_val = np.power(10, -1.7 - (2.65 + moon_factor1 + moon_factor2) * erf_portion) + min_val = np.power(10, -4.0 - (2.95 + moon_factor2) * erf_portion) + + # Update from Curtis Seaman, increase max radiance curve until less than 0.5% is saturated + if self.saturation_correction: + saturation_pct = float(np.count_nonzero(dnb_data > max_val)) / dnb_data.size + LOG.debug("Dynamic DNB saturation percentage: %f", saturation_pct) + while saturation_pct > 0.005: + max_val *= 1.1 + saturation_pct = float(np.count_nonzero(dnb_data > max_val)) / dnb_data.size + LOG.debug("Dynamic DNB saturation percentage: %f", saturation_pct) + + inner_sqrt = (dnb_data - min_val) / (max_val - min_val) + # clip negative values to 0 before the sqrt + inner_sqrt[inner_sqrt < 0] = 0 + np.sqrt(inner_sqrt, out=output_dataset) + + info = {} + info.update(self.info) + info["area"] = dnb_data.info["area"] + info["start_time"] = dnb_data.info["start_time"] + info["end_time"] = dnb_data.info["end_time"] + info["name"] = self.info["name"] + info.setdefault("standard_name", "equalized_radiance") + info.setdefault("wavelength_range", self.info.get("wavelength_range", dnb_data.info.get("wavelength_range"))) + info.setdefault("mode", "L") + output_dataset.info = info + return output_dataset + + def make_day_night_masks(solarZenithAngle, good_mask, highAngleCutoff, lowAngleCutoff, stepsDegrees=None): """ given information on the solarZenithAngle for each point, @@ -318,6 +412,7 @@ def histogram_equalization (data, mask_to_equalize, return out + def local_histogram_equalization (data, mask_to_equalize, valid_data_mask=None, number_of_bins=1000, std_mult_cutoff=3.0, do_zerotoone_normalization=True, diff --git a/mpop/scene.py b/mpop/scene.py index b9659c3f..d03b2f78 100644 --- a/mpop/scene.py +++ b/mpop/scene.py @@ -97,7 +97,7 @@ def read_composites_config(self, composite_configs=None, sensor=None, names=None if len(options["sensor"]) == 1: # FIXME: Finalize how multiple sensors and platforms work options["sensor"] = options["sensor"].pop() - comp_cls = options.get("compositor", None) + comp_cls = options.pop("compositor", None) if not comp_cls: raise ValueError("'compositor' missing or empty in config files: %s" % (composite_configs,)) @@ -114,6 +114,32 @@ def read_composites_config(self, composite_configs=None, sensor=None, names=None LOG.warning("Duplicate composite found, previous composite '%s' will be overwritten", options["name"]) + # Pull out prerequisites + if "prerequisites" in options: + prerequisites = options["prerequisites"].split(",") + options["prerequisites"] = [] + for prerequisite in prerequisites: + try: + # prereqs can be wavelengths + options["prerequisites"].append(float(prerequisite)) + except ValueError: + # or names + options["prerequisites"].append(prerequisite) + + if "optional_prerequisites" in options: + prerequisites = options["optional_prerequisites"].split(",") + options["prerequisites"] = [] + for prerequisite in prerequisites: + try: + # prereqs can be wavelengths + kwargs["optional_prerequisites"].append(float(prerequisite)) + except ValueError: + # or names + kwargs["optional_prerequisites"].append(prerequisite) + + if "metadata_requirements" in options: + options["metadata_requirements"] = options["metadata_requirements"].split(",") + try: loader = runtime_import(comp_cls) except ImportError: