adds ERF DNB composite and updates compositor base to allow for metad…

…ata and optional requirements although they are not completely used yet

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
 

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
 

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):

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,

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: