Wrote an override of the matplotlib LightSource class to remove the intensity scaling step which was causing the movies of the surface evolution to flicker.

python/ctem/ctem/util.py

@@ -5,6 +5,62 @@
 import matplotlib.cm as cm
 import matplotlib.pyplot as plt
 
+
+class CTEMLightSource(LightSource):
+    """
+    Override one function in LightSource to prevent the contrast from being rescaled.
+    """
+    def shade_normals(self, normals, fraction=1.):
+        """
+        Calculate the illumination intensity for the normal vectors of a
+        surface using the defined azimuth and elevation for the light source.
+
+        Imagine an artificial sun placed at infinity in some azimuth and
+        elevation position illuminating our surface. The parts of the surface
+        that slope toward the sun should brighten while those sides facing away
+        should become darker.
+
+        Changes by David Minton: The matplotlib version of this rescales the intensity
+        in a way that causes the brightness level of the images to change between
+        simulation outputs. This makes movies of the surface evolution appear to flicker.
+
+        Parameters
+        ----------
+        fraction : number, optional
+            Increases or decreases the contrast of the hillshade.  Values
+            greater than one will cause intermediate values to move closer to
+            full illumination or shadow (and clipping any values that move
+            beyond 0 or 1). Note that this is not visually or mathematically
+            the same as vertical exaggeration.
+
+        Returns
+        -------
+        ndarray
+            A 2D array of illumination values between 0-1, where 0 is
+            completely in shadow and 1 is completely illuminated.
+        """
+
+        intensity = normals.dot(self.direction)
+
+        # Apply contrast stretch
+        imin, imax = intensity.min(), intensity.max()
+        intensity *= fraction
+
+        # Rescale to 0-1, keeping range before contrast stretch
+        # If constant slope, keep relative scaling (i.e. flat should be 0.5,
+        # fully occluded 0, etc.)
+        # if (imax - imin) > 1e-6:
+        #     # Strictly speaking, this is incorrect. Negative values should be
+        #     # clipped to 0 because they're fully occluded. However, rescaling
+        #     # in this manner is consistent with the previous implementation and
+        #     # visually appears better than a "hard" clip.
+        #     intensity -= imin
+        #     intensity /= (imax - imin)
+        intensity = np.clip(intensity, 0, 1)
+
+        return intensity
+
+
 # Set pixel scaling common for image writing, at 1 pixel/ array element
 dpi = 72.0
 
@@ -125,16 +181,17 @@ def create_rplot(user, odist, pdist, tdist, ph1):
 
     return
 
+
 def image_dem(user, DEM):
     dpi = 300.0  # 72.0
     pix = user['pix']
     gridsize = user['gridsize']
     ve = 1.0
     azimuth = 300.0  # user['azimuth']
-    solar_angle = 15.0  # user['solar_angle']
+    solar_angle = 20.0  # user['solar_angle']
 
-    ls = LightSource(azdeg=azimuth, altdeg=solar_angle)
-    dem_img = ls.hillshade(DEM, vert_exag=ve, dx=pix, dy=pix)
+    ls = CTEMLightSource(azdeg=azimuth, altdeg=solar_angle)
+    dem_img = ls.hillshade(DEM, vert_exag=ve, dx=pix, dy=pix, fraction=1)
 
     # Generate image to put into an array
     height = gridsize / dpi
@@ -180,11 +237,11 @@ def image_shaded_relief(user, DEM):
     pix = user['pix']
     gridsize = user['gridsize']
     ve = 1.0
-    mode = 'hsv'
+    mode = 'overlay'
     azimuth = 300.0  # user['azimuth']
-    solar_angle = 15.0  # user['solar_angle']
+    solar_angle = 20.0  # user['solar_angle']
 
-    ls = LightSource(azdeg=azimuth, altdeg=solar_angle)
+    ls = CTEMLightSource(azdeg=azimuth, altdeg=solar_angle)
     cmap = cm.cividis
 
     # If min and max appear to be reversed, then fix them