Merge remote-tracking branch 'origin/develop' into develop

.gitignore

@@ -28,8 +28,6 @@ __pycache__/
 
 *.dat
 
-*.png
-
 *.DS_Store
 
 alocal.m4
@@ -38,8 +36,6 @@ alocal.m4
 
 *.eps
 
-*.jpg
-
 examples/mare-with-rays-model2/CTEM
 
 examples/global-lunar-bombardment/CTEM

examples/global-lunar-bombardment/console/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/global-lunar-bombardment/frames/.gitignore

@@ -0,0 +1,3 @@
+*.jpg
+*.png
+

examples/global-lunar-bombardment/misc/.gitignore

@@ -0,0 +1 @@
+*.dat

examples/global-lunar-bombardment/rego/.gitignore

@@ -0,0 +1,3 @@
+*.jpg
+*.png
+

examples/global-lunar-bombardment/shaded/.gitignore

@@ -0,0 +1,3 @@
+*.jpg
+*.png
+

examples/global-lunar-bombardment/surf/.gitignore

@@ -0,0 +1,3 @@
+*.jpg
+*.png
+

examples/mare-with-rays-model2/.gitignore

@@ -0,0 +1 @@
+*.idea

examples/mare-with-rays-model2/console/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/mare-with-rays-model2/ctem_io_writers.py

@@ -15,6 +15,8 @@
 import shutil
 import scipy
 from scipy import signal
+from matplotlib.colors import LightSource
+import matplotlib.cm as cm
 
 #Set pixel scaling common for image writing, at 1 pixel/ array element
 dpi = 72.0
@@ -79,112 +81,100 @@ def copy_dists(parameters):
 #Possible references
 #http://nbviewer.jupyter.org/github/ThomasLecocq/geophysique.be/blob/master/2014-02-25%20Shaded%20Relief%20Map%20in%20Python.ipynb   
 
-def image_dem(parameters, surface_dem):
-
-    #Create surface dem map
-    solar_angle = 20.0
-    dem_map = numpy.copy(surface_dem) - numpy.roll(surface_dem, 1, 0)
-    dem_map = (0.5 * numpy.pi) + numpy.arctan2(dem_map, parameters['pix'])
-    dem_map = dem_map - numpy.radians(solar_angle) * (0.5 * numpy.pi)
-    numpy.place(dem_map, dem_map > (0.5 * numpy.pi), 0.5  *numpy.pi)
-    dem_map = numpy.absolute(dem_map)
-    dem_map = 254.0 * numpy.cos(dem_map)
-
-    #Save image to file
+def image_dem(parameters, DEM):
+    dpi = 300.0  # 72.0
+    pix = parameters['pix']
+    gridsize = parameters['gridsize']
+    ve = 1.0
+    azimuth = 300.0 #parameters['azimuth']
+    solar_angle = 20.0 #parameters['solar_angle']
+
+    ls = LightSource(azdeg=azimuth, altdeg=solar_angle)
+    dem_img = ls.hillshade(DEM, vert_exag=ve, dx=pix, dy=pix)
+
+    # Generate image to put into an array
+    height = gridsize / dpi
+    width = gridsize / dpi
+    fig = matplotlib.pyplot.figure(figsize=(width, height), dpi=dpi)
+    ax = matplotlib.pyplot.axes([0, 0, 1, 1])
+    ax.imshow(dem_img, interpolation="nearest", cmap='gray',vmin=0.0,vmax=1.0)
+    matplotlib.pyplot.axis('off')
+    # Save image to file
     filename = parameters['workingdir'] + 'surf' + os.sep + "surf%06d.png" % parameters['ncount']
-    height = parameters['gridsize'] / dpi
-    width = height
-    fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi)
-    fig.figimage(dem_map, cmap = matplotlib.cm.gray, origin = 'lower')
-    matplotlib.pyplot.savefig(filename)    
-
-    return   
-
-def image_regolith(parameters, regolith):
-
-    #Create scaled regolith image
-    minref = parameters['pix'] * 1.0e-4
-    maxreg = numpy.amax(regolith)
-    minreg = numpy.amin(regolith)
-    if (minreg < minref): minreg = minref
-    if (maxreg < minref): maxreg = (minref + 1.0e3)
-    regolith_scaled = numpy.copy(regolith)
-    numpy.place(regolith_scaled, regolith_scaled < minref, minref)
-    regolith_scaled = 254.0 * ((numpy.log(regolith_scaled) - numpy.log(minreg)) / (numpy.log(maxreg) - numpy.log(minreg)))    
-
-    #Save image to file
-    filename = parameters['workingdir'] + 'rego' + os.sep + "rego%06d.png" % parameters['ncount']
-    height = parameters['gridsize'] / dpi
-    width = height
-    fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi)
-    fig.figimage(regolith_scaled, cmap = matplotlib.cm.nipy_spectral, origin = 'lower')
-    matplotlib.pyplot.savefig(filename)   
-
+    matplotlib.pyplot.savefig(filename,dpi=dpi,bbox_inches=0)
+
     return
-
-def image_shaded_relief(parameters, surface_dem):
-    #The color scale and appearance of this do not quite match the IDL version
 
-    #Create image by convolving DEM with 3x3 illumination matrix
-    light = numpy.array([[1.0, 1.0, 1.0], [0.0, 0.0, 0.0], [-1.0, -1.0, -1.0]])
-    convolved_map = scipy.signal.convolve2d(surface_dem, light, mode = 'same')
-
-    #Adjust output to resemble IDL (north slopes illuminated, south in shadow)
-    convolved_map = convolved_map * -1.0
-    convolved_map[0,:]=0.0
-    convolved_map[-1,:] =0.0
-    convolved_map[:,0]=0.0
-    convolved_map[:,-1]=0.0    
+def image_shaded_relief(parameters, DEM):
+    dpi = 300.0 #72.0
+    pix = parameters['pix']
+    gridsize = parameters['gridsize']
+    ve = 1.0
+    mode = 'overlay'
+    azimuth = 300.0 #parameters['azimuth']
+    solar_angle = 20.0  #parameters['solar_angle']
 
-    #If no shadedmin/max parameters are read in from ctem.dat, determine the values from the data
-    if (parameters['shadedminhdefault'] == 1): shadedminh = numpy.amin(surface_dem)
-    if (parameters['shadedmaxhdefault'] == 1): shadedmaxh = numpy.amax(surface_dem)
+    ls = LightSource(azdeg=azimuth, altdeg=solar_angle)
+    cmap = cm.cividis
 
-    #If min and max appear to be reversed, then fix them
+    # If min and max appear to be reversed, then fix them
     if (parameters['shadedminh'] > parameters['shadedmaxh']):
         temp = parameters['shadedminh']
-        shadedminh = parameters['shadedmaxh']
-        shadedmaxh = temp
+        parameters['shadedminh'] = parameters['shadedmaxh']
+        parameters['shadedmaxh'] = temp
+    else:
+        parameters['shadedminh'] = parameters['shadedminh']
+        parameters['shadedmaxh'] = parameters['shadedmaxh']
+
+    # If no shadedmin/max parameters are read in from ctem.dat, determine the values from the data
+    if (parameters['shadedminhdefault'] == 1):
+        shadedminh = numpy.amin(DEM)
     else:
         shadedminh = parameters['shadedminh']
-        shadedmaxh = parameters['shadedmaxh']
-
-    #If dynamic range is valid, construct a shaded DEM
-    dynamic_range = shadedmaxh - shadedminh
-    if  (dynamic_range != 0):
-        dem_scaled = numpy.copy(surface_dem) - shadedminh
-        numpy.place(dem_scaled, dem_scaled < 0.0, 0.0)
-        numpy.place(dem_scaled, dem_scaled > dynamic_range, dynamic_range)
-        dem_scaled = dem_scaled / dynamic_range
+    if (parameters['shadedmaxhdefault'] == 1):
+        shadedmaxh = numpy.amax(DEM)
     else:
-        dem_scaled = numpy.copy(surface_dem) * 0.0
-
-    #Generate shaded depth map with surface_dem color scaling (RGBA)
-    shaded = scipy.misc.bytescale(convolved_map)
-    if numpy.amax(shaded) == 0:  shaded=255
-    shadedscl = shaded / 255.0
-
-    #shaded_imagearr = matplotlib.cm.jet(dem_scaled)
-    #print dem_scaled[0:4,0:4]
-    #shaded_imagearr[:,:,0] = shaded_imagearr[:,:,0] * shadedscl
-    #shaded_imagearr[:,:,1] = shaded_imagearr[:,:,1] * shadedscl
-    #shaded_imagearr[:,:,2] = shaded_imagearr[:,:,2] * shadedscl
-    #shaded_imagearr[:,:,3] = shaded_imagearr[:,:,3] * shadedscl
-
-    #Delivers nearly proper coloring, but no shaded relief    
-    shaded_imagearr = dem_scaled * shadedscl
-    shaded_imagearr = matplotlib.cm.jet(shaded_imagearr)
-    shaded_imagearr = numpy.around(shaded_imagearr, decimals = 1)
-
-    #Save image to file
+        shadedmaxh = parameters['shadedmaxh']
+
+    dem_img = ls.shade(DEM, cmap=cmap,blend_mode=mode, fraction=1.0,
+                       vert_exag=ve, dx=pix, dy=pix,
+                       vmin=shadedminh, vmax=shadedmaxh)
+
+    #Generate image to put into an array
+    height = gridsize / dpi
+    width = gridsize / dpi
+    fig = matplotlib.pyplot.figure(figsize=(width, height), dpi=dpi)
+    ax = matplotlib.pyplot.axes([0, 0, 1, 1])
+    ax.imshow(dem_img,interpolation="nearest",vmin=0.0,vmax=1.0)
+    matplotlib.pyplot.axis('off')
+    # Save image to file
     filename = parameters['workingdir'] + 'shaded' + os.sep + "shaded%06d.png" % parameters['ncount']
+    matplotlib.pyplot.savefig(filename,dpi=dpi,bbox_inches=0)
+    return
+
+
+def image_regolith(parameters, regolith):
+    # Create scaled regolith image
+    minref = parameters['pix'] * 1.0e-4
+    maxreg = numpy.amax(regolith)
+    minreg = numpy.amin(regolith)
+    if (minreg < minref): minreg = minref
+    if (maxreg < minref): maxreg = (minref + 1.0e3)
+    regolith_scaled = numpy.copy(regolith)
+    numpy.place(regolith_scaled, regolith_scaled < minref, minref)
+    regolith_scaled = 254.0 * (
+                (numpy.log(regolith_scaled) - numpy.log(minreg)) / (numpy.log(maxreg) - numpy.log(minreg)))
+
+    # Save image to file
+    filename = parameters['workingdir'] + 'rego' + os.sep + "rego%06d.png" % parameters['ncount']
     height = parameters['gridsize'] / dpi
     width = height
-    fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi)
-    fig.figimage(shaded_imagearr, cmap = matplotlib.cm.jet, origin = 'lower')
-    matplotlib.pyplot.savefig(filename)    
+    fig = matplotlib.pyplot.figure(figsize=(width, height), dpi=dpi)
+    fig.figimage(regolith_scaled, cmap=matplotlib.cm.nipy_spectral, origin='lower')
+    matplotlib.pyplot.savefig(filename)
+
     return
-    
+
 def create_rplot(parameters,odist,pdist,tdist,ph1):
     #Parameters: empirical saturation limit and dfrac
     satlimit = 3.12636

examples/mare-with-rays-model2/dist/.gitignore

@@ -0,0 +1 @@
+*.dat

examples/mare-with-rays-model2/frames/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/mare-with-rays-model2/misc/.gitignore

@@ -0,0 +1 @@
+*.dat

examples/mare-with-rays-model2/rego/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/mare-with-rays-model2/shaded/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/mare-with-rays-model2/surf/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/morphology_test_cases/complex/console/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/morphology_test_cases/complex/dist/.gitignore

@@ -0,0 +1 @@
+*.dat

examples/morphology_test_cases/complex/frames/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/morphology_test_cases/complex/misc/.gitignore

@@ -0,0 +1 @@
+*.dat

examples/morphology_test_cases/complex/rego/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/morphology_test_cases/complex/shaded/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/morphology_test_cases/complex/surf/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/morphology_test_cases/simple/console/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/morphology_test_cases/simple/dist/.gitignore

@@ -0,0 +1 @@
+*.dat

examples/morphology_test_cases/simple/frames/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/morphology_test_cases/simple/misc/.gitignore

@@ -0,0 +1 @@
+*.dat

examples/morphology_test_cases/simple/rego/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/morphology_test_cases/simple/shaded/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/morphology_test_cases/simple/surf/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/ray-generation-test-environment/dist/.gitignore

@@ -0,0 +1 @@
+*.dat

examples/ray-generation-test-environment/frames/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/ray-generation-test-environment/misc/.gitignore

@@ -0,0 +1 @@
+*.dat

examples/ray-generation-test-environment/rego/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/ray-generation-test-environment/shaded/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/ray-generation-test-environment/surf/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/regolith-transport-run-example/dist/.gitignore

@@ -0,0 +1 @@
+*.dat

examples/regolith-transport-run-example/frames/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/regolith-transport-run-example/misc/.gitignore

@@ -0,0 +1 @@
+*.dat

examples/regolith-transport-run-example/rego/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/regolith-transport-run-example/shaded/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png

examples/regolith-transport-run-example/surf/.gitignore

@@ -0,0 +1,2 @@
+*.jpg
+*.png