From 9512de6df6a4315a46fa9ab91fa8fe234df200c9 Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 18 Feb 2022 10:40:37 -0500 Subject: [PATCH 01/19] Tweaked order of vertices when generating faces --- python/ctem/ctem/ctem_viewer_3d.py | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/python/ctem/ctem/ctem_viewer_3d.py b/python/ctem/ctem/ctem_viewer_3d.py index be967d10..4d545b69 100644 --- a/python/ctem/ctem/ctem_viewer_3d.py +++ b/python/ctem/ctem/ctem_viewer_3d.py @@ -68,25 +68,39 @@ def ctem2mesh(self, surface_file): for i in range(s - 1): idx = (s - 1) * j + i # faces[idx,:] = [ j*s+i, j*s+i+1, (j+1)*s+i+1, (j+1)*s+i ] - faces[idx, :] = [j * s + i, j * s + i + 1, (j + 1) * s + i] + faces[idx, :] = [j * s + i, j * s + i + 1, (j + 1) * s + i + 1] idx += (s - 1) ** 2 - faces[idx, :] = [j * s + i + 1, (j + 1) * s + i + 1, (j + 1) * s + i] + faces[idx, :] = [(j + 1) * s + i + 1, (j + 1) * s + i, j * s + i ] self.mesh.vertices = open3d.utility.Vector3dVector(verts) self.mesh.triangles = open3d.utility.Vector3iVector(faces) self.mesh.compute_vertex_normals() + return + def visualize(self): + vis = open3d.visualization.Visualizer() + vis.create_window() + vis.add_geometry(self.mesh) + + + # zmax = np.amax(surf.dem) + # zmin = np.amin(surf.dem) + # cnorm = (surf.dem.flatten() - zmin) / (zmax - zmin) + # cval = plt.cm.terrain(cnorm)[:, :3] + + self.mesh.paint_uniform_color([0.5, 0.5, 0.5]) + vis.run() + vis.destroy_window() + + if __name__ == '__main__': import matplotlib.pyplot as plt surf = polysurface() surf.ctem2mesh(surface_file="surface_dem.dat") + surf.visualize() + + - zmax = np.amax(surf.dem) - zmin = np.amin(surf.dem) - cnorm = (surf.dem.flatten() - zmin) / (zmax - zmin) - cval = plt.cm.terrain(cnorm)[:,:3] - surf.mesh.vertex_colors = open3d.utility.Vector3dVector(cval) - open3d.visualization.draw_geometries([surf.mesh]) From 2fa61db9c477ad43c47f9ae1a11ea36143191aa0 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 18 Feb 2022 10:41:26 -0500 Subject: [PATCH 02/19] Added setup script so that the ctem Python packaage may be installed --- python/ctem/setup.py | 6 ++++++ 1 file changed, 6 insertions(+) create mode 100644 python/ctem/setup.py diff --git a/python/ctem/setup.py b/python/ctem/setup.py new file mode 100644 index 00000000..f7003601 --- /dev/null +++ b/python/ctem/setup.py @@ -0,0 +1,6 @@ +from setuptools import setup, find_packages + +setup(name='ctem', + version='0.1', + author='David A. Minton', + packages=find_packages()) From 528a1182098c69594da6478990ba190dd438cfa6 Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 18 Feb 2022 12:47:33 -0500 Subject: [PATCH 03/19] Added real2float function to convert Fortran double precision real ASCII that uses d or D for exponent to float with E --- .gitignore | 2 ++ python/ctem/ctem/ctem_io_readers.py | 41 +++++++++++++++++++++-------- 2 files changed, 32 insertions(+), 11 deletions(-) diff --git a/.gitignore b/.gitignore index 1650087f..f7ab64c6 100644 --- a/.gitignore +++ b/.gitignore @@ -69,3 +69,5 @@ examples/global-lunar-bombardment/scale.ipynb *.code-workspace *.png + +python/ctem/ctem.egg-info/ diff --git a/python/ctem/ctem/ctem_io_readers.py b/python/ctem/ctem/ctem_io_readers.py index 84cf7e06..0935fac9 100644 --- a/python/ctem/ctem/ctem_io_readers.py +++ b/python/ctem/ctem/ctem_io_readers.py @@ -12,6 +12,25 @@ import numpy + +def real2float(realstr): + """ + Converts a Fortran-generated ASCII string of a real value into a numpy float type. Handles cases where double precision + numbers in exponential notation use 'd' or 'D' instead of 'e' or 'E' + + Parameters + ---------- + realstr : string + Fortran-generated ASCII string of a real value. + + Returns + ------- + : float + The converted floating point value of the input string + """ + return float(realstr.replace('d', 'E').replace('D', 'E')) + + def read_ctemin(parameters,notset): #Read and parse ctem.in file inputfile = parameters['workingdir'] + parameters['ctemfile'] @@ -26,14 +45,14 @@ def read_ctemin(parameters,notset): for line in lines: fields = line.split() if len(fields) > 0: - if ('pix' == fields[0].lower()): parameters['pix']=float(fields[1]) + if ('pix' == fields[0].lower()): parameters['pix']=real2float(fields[1]) if ('gridsize' == fields[0].lower()): parameters['gridsize']=int(fields[1]) if ('seed' == fields[0].lower()): parameters['seed']=int(fields[1]) if ('sfdfile' == fields[0].lower()): parameters['sfdfile']=fields[1] if ('impfile' == fields[0].lower()): parameters['impfile']=fields[1] - if ('maxcrat' == fields[0].lower()): parameters['maxcrat']=float(fields[1]) + if ('maxcrat' == fields[0].lower()): parameters['maxcrat']=real2float(fields[1]) if ('sfdcompare' == fields[0].lower()): parameters['sfdcompare']=fields[1] - if ('interval' == fields[0].lower()): parameters['interval']=float(fields[1]) + if ('interval' == fields[0].lower()): parameters['interval']=real2float(fields[1]) if ('numintervals' == fields[0].lower()): parameters['numintervals']=int(fields[1]) if ('popupconsole' == fields[0].lower()): parameters['popupconsole']=fields[1] if ('saveshaded' == fields[0].lower()): parameters['saveshaded']=fields[1] @@ -43,10 +62,10 @@ def read_ctemin(parameters,notset): if ('runtype' == fields[0].lower()): parameters['runtype']=fields[1] if ('restart' == fields[0].lower()): parameters['restart']=fields[1] if ('shadedminh' == fields[0].lower()): - parameters['shadedminh'] = float(fields[1]) + parameters['shadedminh'] = real2float(fields[1]) parameters['shadedminhdefault'] = 0 if ('shadedmaxh' == fields[0].lower()): - parameters['shadedmaxh'] = float(fields[1]) + parameters['shadedmaxh'] = real2float(fields[1]) parameters['shadedmaxhdefault'] = 0 #Test values for further processing @@ -109,19 +128,19 @@ def read_ctemdat(parameters, seedarr): #Parse file lines and update parameter fields fields = lines[0].split() if len(fields) > 0: - parameters['totalimpacts'] = float(fields[0]) + parameters['totalimpacts'] = real2float(fields[0]) parameters['ncount'] = int(fields[1]) - parameters['curyear'] = float(fields[2]) + parameters['curyear'] = real2float(fields[2]) parameters['restart'] = fields[3] - parameters['fracdone'] = float(fields[4]) - parameters['masstot'] = float(fields[5]) + parameters['fracdone'] = real2float(fields[4]) + parameters['masstot'] = real2float(fields[5]) #Parse remainder of file to build seed array nlines = len(lines) index = 1 while (index < nlines): fields = lines[index].split() - seedarr[index - 1] = float(fields[0]) + seedarr[index - 1] = real2float(fields[0]) index += 1 parameters['seedn'] = index - 1 @@ -137,7 +156,7 @@ def read_impact_mass(filename): fields = line[0].split() if (len(fields) > 0): - mass = float(fields[0]) + mass = real2float(fields[0]) else: mass = 0 From b273acba928270961c125284674f2280eb4a2527 Mon Sep 17 00:00:00 2001 From: David Minton Date: Sat, 19 Feb 2022 14:23:06 -0500 Subject: [PATCH 04/19] Updated the viewer to plot elevation points as square pixels as it is characterized internally --- python/ctem/ctem/ctem_viewer_3d.py | 78 ++++++++++++++++++++------ python/ctem/tests/viewer3d/polytest.py | 13 +++-- 2 files changed, 69 insertions(+), 22 deletions(-) diff --git a/python/ctem/ctem/ctem_viewer_3d.py b/python/ctem/ctem/ctem_viewer_3d.py index 4d545b69..093ae1e3 100644 --- a/python/ctem/ctem/ctem_viewer_3d.py +++ b/python/ctem/ctem/ctem_viewer_3d.py @@ -1,7 +1,7 @@ import numpy as np import os import open3d -import ctem_io_readers +from ctem import ctem_io_readers import os @@ -57,31 +57,76 @@ def ctem2mesh(self, surface_file): # Build mesh grid s = parameters['gridsize'] pix = parameters['pix'] - ygrid, xgrid = np.mgrid[-s / 2:s / 2, -s / 2:s / 2] * pix + ygrid, xgrid = np.mgrid[-s/2:s/2, -s/2:s/2] * pix + y0, x0 = ygrid - pix/2, xgrid - pix/2 + y1, x1 = ygrid - pix/2, xgrid + pix/2 + y2, x2 = ygrid + pix/2, xgrid + pix/2 + y3, x3 = ygrid + pix/2, xgrid - pix/2 + - xvals = xgrid.flatten() - yvals = ygrid.flatten() - zvals = self.dem.flatten() + xvals = np.append( + np.append( + np.append(x0.flatten(), + x1.flatten()), + x2.flatten()), + x3.flatten()) + yvals = np.append( + np.append( + np.append(y0.flatten(), + y1.flatten()), + y2.flatten()), + y3.flatten()) + zvals = np.append( + np.append( + np.append(self.dem.flatten(), + self.dem.flatten()), + self.dem.flatten()), + self.dem.flatten()) verts = np.array((xvals, yvals, zvals)).T - faces = np.empty([2 * (s - 1) ** 2, 3], dtype=np.int64) - for j in range(s - 1): - for i in range(s - 1): - idx = (s - 1) * j + i - # faces[idx,:] = [ j*s+i, j*s+i+1, (j+1)*s+i+1, (j+1)*s+i ] - faces[idx, :] = [j * s + i, j * s + i + 1, (j + 1) * s + i + 1] - idx += (s - 1) ** 2 - faces[idx, :] = [(j + 1) * s + i + 1, (j + 1) * s + i, j * s + i ] + nface_triangles = 10 + faces = np.full([nface_triangles*s**2, 3], -1, dtype=np.int64) + for j in range(s): + for i in range(s): + i0 = s*j + i + i1 = i0 + s**2 + i2 = i1 + s**2 + i3 = i2 + s**2 + + fidx = np.arange(nface_triangles*i0,nface_triangles*(i0+1)) + # Make the two top triangles + faces[fidx[0],:] = [i0, i1, i2] + faces[fidx[1],:] = [i0, i2, i3] + # Make the two west side triangles + if i > 0: + faces[fidx[0],:] = [i0, i3, i3-1] + faces[fidx[1],:] = [i0, i3-1, i0-1] + # Make the two south side triangles + if j > 0: + faces[fidx[4],:] = [i1, i0, i3-s ] + faces[fidx[5],:] = [i1, i3-s, i2-s] + # Make the two east side triangles + if i < (s - 1): + faces[fidx[6],:] = [i2, i1, i0+1] + faces[fidx[7],:] = [i2, i0+1, i3+1] + #Make the two north side triangles + if j < (s -1): + faces[fidx[8],:] = [i3, i2, i1+s] + faces[fidx[9],:] = [i3, i1+s, i0+s] + nz = faces[:,0] != -1 + f2 = faces[nz,:] self.mesh.vertices = open3d.utility.Vector3dVector(verts) - self.mesh.triangles = open3d.utility.Vector3iVector(faces) + self.mesh.triangles = open3d.utility.Vector3iVector(f2) self.mesh.compute_vertex_normals() - + self.mesh.compute_triangle_normals() + return def visualize(self): vis = open3d.visualization.Visualizer() vis.create_window() vis.add_geometry(self.mesh) - + opt = vis.get_render_option() + opt.background_color = np.asarray([0, 0, 0]) # zmax = np.amax(surf.dem) # zmin = np.amin(surf.dem) @@ -99,6 +144,7 @@ def visualize(self): surf = polysurface() surf.ctem2mesh(surface_file="surface_dem.dat") surf.visualize() + #open3d.visualization.draw_geometries_with_editing([surf.mesh]) diff --git a/python/ctem/tests/viewer3d/polytest.py b/python/ctem/tests/viewer3d/polytest.py index 13cee521..c528374c 100644 --- a/python/ctem/tests/viewer3d/polytest.py +++ b/python/ctem/tests/viewer3d/polytest.py @@ -1,8 +1,9 @@ import numpy as np from ctem import ctem_viewer_3d - -surf = ctem_viewer_3d.surf() -surf.io_read_ply("216kleopatra-mesh.ply") -v = np.asarray(surf.mesh) -t = np.asarray(surf.mesh) -surf.render() \ No newline at end of file +import open3d +surf = ctem_viewer_3d.polysurface() +#surf.io_read_ply("216kleopatra-mesh.ply") +surf.mesh = open3d.geometry.TriangleMesh.create_box() +v = np.asarray(surf.mesh.vertices) +t = np.asarray(surf.mesh.triangles) +surf.visualize() \ No newline at end of file From 1b9e858a830f81ec6486d4de14393af583b21b14 Mon Sep 17 00:00:00 2001 From: David Minton Date: Sat, 19 Feb 2022 14:34:36 -0500 Subject: [PATCH 05/19] Fixed bad indices --- python/ctem/ctem/ctem_viewer_3d.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/ctem/ctem/ctem_viewer_3d.py b/python/ctem/ctem/ctem_viewer_3d.py index 093ae1e3..45612354 100644 --- a/python/ctem/ctem/ctem_viewer_3d.py +++ b/python/ctem/ctem/ctem_viewer_3d.py @@ -98,8 +98,8 @@ def ctem2mesh(self, surface_file): faces[fidx[1],:] = [i0, i2, i3] # Make the two west side triangles if i > 0: - faces[fidx[0],:] = [i0, i3, i3-1] - faces[fidx[1],:] = [i0, i3-1, i0-1] + faces[fidx[2],:] = [i0, i3, i2-1] + faces[fidx[3],:] = [i0, i2-1, i1-1] # Make the two south side triangles if j > 0: faces[fidx[4],:] = [i1, i0, i3-s ] From 93f168dfc23b01c9c02a1b97cf808a535697f466 Mon Sep 17 00:00:00 2001 From: David Minton Date: Sat, 19 Feb 2022 15:09:17 -0500 Subject: [PATCH 06/19] Added two different mesh styles (block vs triangle) and separated input file reading out into its own method --- python/ctem/ctem/ctem_viewer_3d.py | 50 +++++++++++++++++++++++------- 1 file changed, 38 insertions(+), 12 deletions(-) diff --git a/python/ctem/ctem/ctem_viewer_3d.py b/python/ctem/ctem/ctem_viewer_3d.py index 45612354..b4517054 100644 --- a/python/ctem/ctem/ctem_viewer_3d.py +++ b/python/ctem/ctem/ctem_viewer_3d.py @@ -13,13 +13,12 @@ def __init__(self): self.mesh = open3d.geometry.TriangleMesh() return - def ctem2mesh(self, surface_file): - + def readctem(self, surface_file): # Create and initialize data dictionaries for parameters and options from CTEM.in notset = '-NOTSET-' currentdir = os.getcwd() + os.sep - parameters = {'restart': notset, + self.parameters = {'restart': notset, 'runtype': notset, 'popupconsole': notset, 'saveshaded': notset, @@ -48,22 +47,24 @@ def ctem2mesh(self, surface_file): 'impfile': notset, 'sfdcompare': notset, 'sfdfile': notset} - + # Read ctem.in to initialize parameter values based on user input - ctem_io_readers.read_ctemin(parameters, notset) + ctem_io_readers.read_ctemin(self.parameters, notset) # Read surface dem(shaded relief) and ejecta data files - dem_file = parameters['workingdir'] + surface_file - self.dem = ctem_io_readers.read_unformatted_binary(dem_file, parameters['gridsize']) + dem_file = self.parameters['workingdir'] + surface_file + self.dem = ctem_io_readers.read_unformatted_binary(dem_file, self.parameters['gridsize']) + return + + def ctem2blockmesh(self): # Build mesh grid - s = parameters['gridsize'] - pix = parameters['pix'] + s = self.parameters['gridsize'] + pix = self.parameters['pix'] ygrid, xgrid = np.mgrid[-s/2:s/2, -s/2:s/2] * pix y0, x0 = ygrid - pix/2, xgrid - pix/2 y1, x1 = ygrid - pix/2, xgrid + pix/2 y2, x2 = ygrid + pix/2, xgrid + pix/2 y3, x3 = ygrid + pix/2, xgrid - pix/2 - xvals = np.append( np.append( np.append(x0.flatten(), @@ -121,6 +122,31 @@ def ctem2mesh(self, surface_file): return + def ctem2trimesh(self): + # Build mesh grid + s = self.parameters['gridsize'] + pix = self.parameters['pix'] + ygrid, xgrid = np.mgrid[-s/2:s/2, -s/2:s/2] * pix + + xvals = xgrid.flatten() + yvals = ygrid.flatten() + zvals = self.dem.flatten() + verts = np.array((xvals, yvals, zvals)).T + faces = np.full([2 * (s-1)**2, 3], -1, dtype=np.int64) + for j in range(s - 1): + for i in range(s - 1): + idx = (s - 1) * j + i + # faces[idx,:] = [ j*s+i, j*s+i+1, (j+1)*s+i+1, (j+1)*s+i ] + faces[idx, :] = [j * s + i, j * s + i + 1, (j + 1) * s + i + 1] + idx += (s - 1) ** 2 + faces[idx, :] = [(j + 1) * s + i + 1, (j + 1) * s + i, j * s + i ] + self.mesh.vertices = open3d.utility.Vector3dVector(verts) + self.mesh.triangles = open3d.utility.Vector3iVector(faces) + self.mesh.compute_vertex_normals() + self.mesh.compute_triangle_normals() + + return + def visualize(self): vis = open3d.visualization.Visualizer() vis.create_window() @@ -142,9 +168,9 @@ def visualize(self): import matplotlib.pyplot as plt surf = polysurface() - surf.ctem2mesh(surface_file="surface_dem.dat") + surf.readctem(surface_file="surface_dem.dat") + surf.ctem2trimesh() surf.visualize() - #open3d.visualization.draw_geometries_with_editing([surf.mesh]) From ead9d0cad7c6f1be351275580073aa6fbd49e509 Mon Sep 17 00:00:00 2001 From: David Minton Date: Mon, 21 Feb 2022 11:03:11 -0500 Subject: [PATCH 07/19] Began process of restructuring the Python frontend with a basic class definition for a CTEM simulation --- python/ctem/ctem/__init__.py | 1 + python/ctem/ctem/frontend.py | 46 ++++++++++++++++ python/ctem/ctem/io.py | 87 ++++++++++++++++++++++++++++++ python/ctem/tests/testio/ctem.in | 70 ++++++++++++++++++++++++ python/ctem/tests/testio/testio.py | 7 +++ 5 files changed, 211 insertions(+) create mode 100644 python/ctem/ctem/frontend.py create mode 100644 python/ctem/ctem/io.py create mode 100755 python/ctem/tests/testio/ctem.in create mode 100644 python/ctem/tests/testio/testio.py diff --git a/python/ctem/ctem/__init__.py b/python/ctem/ctem/__init__.py index 4774c6af..b9aa95ac 100644 --- a/python/ctem/ctem/__init__.py +++ b/python/ctem/ctem/__init__.py @@ -1,3 +1,4 @@ # from ctem.ctem_io_readers import * # from ctem.ctem_io_writers import * # from ctem.ctem_driver import * +from ctem.frontend import * \ No newline at end of file diff --git a/python/ctem/ctem/frontend.py b/python/ctem/ctem/frontend.py new file mode 100644 index 00000000..6340d5ec --- /dev/null +++ b/python/ctem/ctem/frontend.py @@ -0,0 +1,46 @@ +import numpy +import os +import subprocess +import shutil +from ctem import io + +class Simulation: + """ + This is a class that defines the basic CTEM simulation object + """ + def __init__(self, param_file="ctem.in"): + currentdir = os.getcwd() + self.parameters = { + 'restart': None, + 'runtype': None, + 'popupconsole': None, + 'saveshaded': None, + 'saverego': None, + 'savepres': None, + 'savetruelist': None, + 'seedn': 1, + 'totalimpacts': 0, + 'ncount': 0, + 'curyear': 0.0, + 'fracdone': 1.0, + 'masstot': 0.0, + 'interval': 0.0, + 'numintervals': 0, + 'pix': -1.0, + 'gridsize': -1, + 'seed': 0, + 'maxcrat': 1.0, + 'shadedminhdefault': 1, + 'shadedmaxhdefault': 1, + 'shadedminh': 0.0, + 'shadedmaxh': 0.0, + 'workingdir': currentdir, + 'ctemfile': param_file, + 'datfile': 'ctem.dat', + 'impfile': None, + 'sfdcompare': None, + 'sfdfile': None + } + self.parameters = io.read_param(self.parameters) + + return \ No newline at end of file diff --git a/python/ctem/ctem/io.py b/python/ctem/ctem/io.py new file mode 100644 index 00000000..664d8bb8 --- /dev/null +++ b/python/ctem/ctem/io.py @@ -0,0 +1,87 @@ +import numpy +import os + +def real2float(realstr): + """ + Converts a Fortran-generated ASCII string of a real value into a numpy float type. Handles cases where double precision + numbers in exponential notation use 'd' or 'D' instead of 'e' or 'E' + + Parameters + ---------- + realstr : string + Fortran-generated ASCII string of a real value. + + Returns + ------- + : float + The converted floating point value of the input string + """ + return float(realstr.replace('d', 'E').replace('D', 'E')) + +def read_param(parameters): + # Read and parse ctem.in file + inputfile = os.path.join(parameters['workingdir'], parameters['ctemfile']) + + # Read ctem.in file + print('Reading input file ' + parameters['ctemfile']) + with open(inputfile, 'r') as fp: + lines = fp.readlines() + + # Process file text + for line in lines: + fields = line.split() + if len(fields) > 0: + if ('pix' == fields[0].lower()): parameters['pix'] = real2float(fields[1]) + if ('gridsize' == fields[0].lower()): parameters['gridsize'] = int(fields[1]) + if ('seed' == fields[0].lower()): parameters['seed'] = int(fields[1]) + if ('sfdfile' == fields[0].lower()): parameters['sfdfile'] = fields[1] + if ('impfile' == fields[0].lower()): parameters['impfile'] = fields[1] + if ('maxcrat' == fields[0].lower()): parameters['maxcrat'] = real2float(fields[1]) + if ('sfdcompare' == fields[0].lower()): parameters['sfdcompare'] = fields[1] + if ('interval' == fields[0].lower()): parameters['interval'] = real2float(fields[1]) + if ('numintervals' == fields[0].lower()): parameters['numintervals'] = int(fields[1]) + if ('popupconsole' == fields[0].lower()): parameters['popupconsole'] = fields[1] + if ('saveshaded' == fields[0].lower()): parameters['saveshaded'] = fields[1] + if ('saverego' == fields[0].lower()): parameters['saverego'] = fields[1] + if ('savepres' == fields[0].lower()): parameters['savepres'] = fields[1] + if ('savetruelist' == fields[0].lower()): parameters['savetruelist'] = fields[1] + if ('runtype' == fields[0].lower()): parameters['runtype'] = fields[1] + if ('restart' == fields[0].lower()): parameters['restart'] = fields[1] + if ('shadedminh' == fields[0].lower()): + parameters['shadedminh'] = real2float(fields[1]) + parameters['shadedminhdefault'] = 0 + if ('shadedmaxh' == fields[0].lower()): + parameters['shadedmaxh'] = real2float(fields[1]) + parameters['shadedmaxhdefault'] = 0 + + # Test values for further processing + if (parameters['interval'] <= 0.0): + print('Invalid value for or missing variable INTERVAL in ' + inputfile) + if (parameters['numintervals'] <= 0): + print('Invalid value for or missing variable NUMINTERVALS in ' + inputfile) + if (parameters['pix'] <= 0.0): + print('Invalid value for or missing variable PIX in ' + inputfile) + if (parameters['gridsize'] <= 0): + print('Invalid value for or missing variable GRIDSIZE in ' + inputfile) + if (parameters['seed'] == 0): + print('Invalid value for or missing variable SEED in ' + inputfile) + if (parameters['sfdfile'] is None): + print('Invalid value for or missing variable SFDFILE in ' + inputfile) + if (parameters['impfile'] is None): + print('Invalid value for or missing variable IMPFILE in ' + inputfile) + if (parameters['popupconsole'] is None): + print('Invalid value for or missing variable POPUPCONSOLE in ' + inputfile) + if (parameters['saveshaded'] is None): + print('Invalid value for or missing variable SAVESHADED in ' + inputfile) + if (parameters['saverego'] is None): + print('Invalid value for or missing variable SAVEREGO in ' + inputfile) + if (parameters['savepres'] is None): + print('Invalid value for or missing variable SAVEPRES in ' + inputfile) + if (parameters['savetruelist'] is None): + print('Invalid value for or missing variable SAVETRUELIST in ' + inputfile) + if (parameters['runtype'] is None): + print('Invalid value for or missing variable RUNTYPE in ' + inputfile) + if (parameters['restart'] is None): + print('Invalid value for or missing variable RESTART in ' + inputfile) + + return parameters diff --git a/python/ctem/tests/testio/ctem.in b/python/ctem/tests/testio/ctem.in new file mode 100755 index 00000000..e57310dc --- /dev/null +++ b/python/ctem/tests/testio/ctem.in @@ -0,0 +1,70 @@ +! CTEM Input file + + +! Testing input. These are used to perform non-Monte Carlo tests. +testflag T ! Set to T to create a single crater with user-defined impactor properties +testimp 5.934e3 ! 93km crater ! Diameter of test impactor (m) +testvel 15.0e3 ! Velocity of test crater (m/s) +testang 90.0 ! Impact angle of test crater (deg) - Default 90.0 +testxoffset 0.0e0 ! x-axis offset of crater center from grid center (m) - Default 0.0 +testyoffset 0.0e0 ! y-axis offset of crater center from grid center (m) - Default 0.0 +tallyonly F ! Tally the craters without generating any craters +testtally F + +! IDL driver in uts +interval 1.0 +numintervals 1 ! Total number of intervals +restart F ! Restart a previous run +impfile NPFextrap.dat ! Impactor SFD rate file (col 1: Dimp (m), col 2: ! impactors > D (m**(-2) y**(-1)) +popupconsole F ! Pop up console window every output interval +saveshaded T ! Output shaded relief images +saverego F ! Output regolith map images +savepres F ! Output simplified console display images (presentation-compatible images) +savetruelist T ! Save the true cumulative crater distribution for each interval (large file size) +runtype statistical ! single: craters accumulate in successive intervals + ! statistical: surface is reset between intervals +sfdcompare FassettCounts.txt + +! CTEM required inputs +seed 23790 ! Random number generator seed +gridsize 1000 ! Size of grid in pixels +numlayers 10 ! Number of perched layers +pix 0.3468773099817e3 ! Pixel size (m) - Copernicus DEM compariso +mat rock ! Material (rock or ice) +! Bedrock scaling parameters +mu_b 0.55e0 ! Experimentally derived parameter for bedrock crater scaling law +kv_b 0.20e0 ! Experimentally derived parameter for bedrock crater scaling law +trho_b 2250.0e0 ! Target density (bedrock) (kg/m**3) +ybar_b 0.0e6 ! Bedrock strength (Pa) +! Regolith scaling parameters +mu_r 0.55e0 ! Experimentally derived parameter for regolith crater scaling law +kv_r 0.20e0 ! Experimentally derived parameter for regolith crater scaling law +trho_r 2250.0e0 ! Target density (regolith) (kg/m**3) +ybar_r 0.00e6 ! Regolith strength (Pa) +! Body parameters +gaccel 1.62e0 ! Gravitational acceleration at target (m/s**2) +trad 1737.35e3 ! Target radius (m) +prho 2500.0e0 ! Projectile density (kg/m**3) +sfdfile production.dat ! Impactor SFD file (col 1: Dimp (m), col 2: ! impactors > D +velfile lunar-MBA-impactor-velocities.dat ! Impactor velocity dist file + +! Seismic shaking input (required if seismic shaking is set to T, otherwise ignored) +doseismic F ! Perform seismic shaking calculations with each impact - Default F + +! Optional inputF These have internally set default values that work reasonable well. Comment them out with +deplimit 9e99 ! Depth limit for craters (m) - Default is to ignore. +maxcrat 1.00e0 ! Fraction of gridsize that maximum crater can be - Default 1.0 +killatmaxcrater F ! Stop the run if a crater larger than the maximum is produced - Default F +basinimp 35.0e3 ! Size of impactor to switch to lunar basin scaling law - Default is to ignore +docollapse T ! Do slope collapse - Default T +dosoftening F ! Do ejecta softening - Default T +doangle T ! Vary the impact angle. Set to F to have only vertical impacts - Default T +Kd1 0.0001 +psi 2.000 +fe 2.00 +ejecta_truncation 2.0 +dorays F +superdomain F + +dorealistic T + diff --git a/python/ctem/tests/testio/testio.py b/python/ctem/tests/testio/testio.py new file mode 100644 index 00000000..ae329313 --- /dev/null +++ b/python/ctem/tests/testio/testio.py @@ -0,0 +1,7 @@ +import ctem + +d = dir(ctem) + +sim = ctem.Simulation() +for k, v in sim.parameters.items(): + print(k, v) \ No newline at end of file From aa12ad186eb82229d0f1cff44ddccceb705066c5 Mon Sep 17 00:00:00 2001 From: David Minton Date: Mon, 21 Feb 2022 11:33:10 -0500 Subject: [PATCH 08/19] Updated frontend with new io routines --- python/ctem/ctem/frontend.py | 83 ++++++++- python/ctem/ctem/io.py | 327 ++++++++++++++++++++++++++++++++++- 2 files changed, 406 insertions(+), 4 deletions(-) diff --git a/python/ctem/ctem/frontend.py b/python/ctem/ctem/frontend.py index 6340d5ec..69327c37 100644 --- a/python/ctem/ctem/frontend.py +++ b/python/ctem/ctem/frontend.py @@ -1,4 +1,4 @@ -import numpy +import numpy as np import os import subprocess import shutil @@ -43,4 +43,83 @@ def __init__(self, param_file="ctem.in"): } self.parameters = io.read_param(self.parameters) - return \ No newline at end of file + # Set up data arrays + self.seedarr = np.zeros(100, dtype=int) + self.seedarr[0] = self.parameters['seed'] + self.odist = np.zeros([1, 6]) + self.pdist = np.zeros([1, 6]) + self.tdist = np.zeros([1, 6]) + self.surface_dem = np.zeros([self.parameters['gridsize'], self.parameters['gridsize']], dtype=np.float) + self.regolith = np.zeros([self.parameters['gridsize'], self.parameters['gridsize']], dtype=np.float) + + if self.parameters['sfdcompare'] is not None: + # Read sfdcompare file + sfdfile = os.path.join(self.parameters['workingdir'], self.parameters['sfdcompare']) + self.ph1 = io.read_formatted_ascii(sfdfile, skip_lines=0) + + # Read production function file + impfile = os.path.join(self.parameters['workingdir'], parameters['impfile']) + prodfunction = io.read_formatted_ascii(impfile, skip_lines=0) + + # Create impactor production population + area = (self.parameters['gridsize'] * self.parameters['pix']) ** 2 + self.production = np.copy(prodfunction) + self.production[:, 1] = self.production[:, 1] * area * self.parameters['interval'] + + # Write corrected production function to file + io.write_production(self.parameters, self.production) + + # Starting new or old run? + if (self.parameters['restart'].upper() == 'F'): + print('Starting a new run') + + if (self.parameters['runtype'].upper() == 'STATISTICAL'): + self.parameters['ncount'] = 1 + + # Write ctem.dat file + io.write_ctemdat(self.parameters, self.seedarr) + + else: + self.parameters['ncount'] = 0 + + # Delete tdistribution file, if it exists + tdist_file = os.path.join(self.parameters['workingdir'], 'tdistribution.dat') + if os.path.isfile(tdist_file): + os.remove(tdist_file) + else: + print('Continuing a previous run') + + # Read surface dem(shaded relief) and ejecta data files + dem_file = os.path.join(self.parameters['workingdir'], 'surface_dem.dat') + self.surface_dem = io.read_unformatted_binary(dem_file, self.parameters['gridsize']) + ejecta_file = os.path.join(self.parameters['workingdir'], 'surface_ejc.dat') + self.regolith = io.read_unformatted_binary(ejecta_file, self.parameters['gridsize']) + + # Read odistribution, tdistribution, and pdistribution files + ofile = os.path.join(self.parameters['workingdir'],'odistribution.dat') + odist = io.read_formatted_ascii(ofile, skip_lines=1) + tfile = os.path.join(self.parameters['workingdir'], 'tdistribution.dat') + tdist = io.read_formatted_ascii(tfile, skip_lines=1) + pfile = os.path.join(self.parameters['workingdir'], 'pdistribution.dat') + pdist = io.read_formatted_ascii(pfile, skip_lines=1) + + # Read impact mass from file + massfile = os.path.join(self.parameters['workingdir'], 'impactmass.dat') + impact_mass = io.read_impact_mass(massfile) + + # Read ctem.dat file + io.read_ctemdat(self.parameters, self.seedarr) + + # Open fracdonefile and regodepthfile for writing + # filename = os.path.join(self.parameters['workingdir'], 'fracdone.dat') + # fp_frac = open(filename, 'w') + # filename = os.path.join(self.parameters['workingdir'], 'regolithdepth.dat') + # fp_reg = open(filename, 'w') + + io.create_dir_structure(self.parameters) + + + return + + + def \ No newline at end of file diff --git a/python/ctem/ctem/io.py b/python/ctem/ctem/io.py index 664d8bb8..3786d87c 100644 --- a/python/ctem/ctem/io.py +++ b/python/ctem/ctem/io.py @@ -1,9 +1,17 @@ -import numpy +import numpy as np import os +import matplotlib.pyplot as plt +import shutil +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 def real2float(realstr): """ - Converts a Fortran-generated ASCII string of a real value into a numpy float type. Handles cases where double precision + Converts a Fortran-generated ASCII string of a real value into a np float type. Handles cases where double precision numbers in exponential notation use 'd' or 'D' instead of 'e' or 'E' Parameters @@ -85,3 +93,318 @@ def read_param(parameters): print('Invalid value for or missing variable RESTART in ' + inputfile) return parameters + + +def read_formatted_ascii(filename, skip_lines): + # Generalized ascii text reader + # For use with sfdcompare, production, odist, tdist, pdist data files + data = np.genfromtxt(filename, skip_header=skip_lines) + return data + + +def read_unformatted_binary(filename, gridsize): + # Read unformatted binary files created by Fortran + # For use with surface ejecta and surface dem data files + dt = np.float + data = np.fromfile(filename, dtype=dt) + data.shape = (gridsize, gridsize) + + return data + + +def read_ctemdat(parameters, seedarr): + # Read and parse ctem.dat file + datfile = parameters['workingdir'] + 'ctem.dat' + + # Read ctem.dat file + print('Reading input file ' + parameters['datfile']) + fp = open(datfile, 'r') + lines = fp.readlines() + fp.close() + + # Parse file lines and update parameter fields + fields = lines[0].split() + if len(fields) > 0: + parameters['totalimpacts'] = real2float(fields[0]) + parameters['ncount'] = int(fields[1]) + parameters['curyear'] = real2float(fields[2]) + parameters['restart'] = fields[3] + parameters['fracdone'] = real2float(fields[4]) + parameters['masstot'] = real2float(fields[5]) + + # Parse remainder of file to build seed array + nlines = len(lines) + index = 1 + while (index < nlines): + fields = lines[index].split() + seedarr[index - 1] = real2float(fields[0]) + index += 1 + + parameters['seedn'] = index - 1 + + return + + +def read_impact_mass(filename): + # Read impact mass file + + fp = open(filename, 'r') + line = fp.readlines() + fp.close() + + fields = line[0].split() + if (len(fields) > 0): + mass = real2float(fields[0]) + else: + mass = 0 + + return mass + + +# Write production function to file production.dat +# This file format does not exactly match that generated from IDL. Does it work? +def write_production(parameters, production): + filename = parameters['workingdir'] + parameters['sfdfile'] + np.savetxt(filename, production, fmt='%1.8e', delimiter=' ') + + return + + +def create_dir_structure(parameters): + # Create directories for various output files if they do not already exist + directories = ['dist', 'misc', 'rego', 'rplot', 'surf', 'shaded'] + + for directory in directories: + dir_test = parameters['workingdir'] + directory + if not os.path.isdir(dir_test): + os.makedirs(dir_test) + + return + + +def write_ctemdat(parameters, seedarr): + # Write various parameters and random number seeds into ctem.dat file + filename = parameters['workingdir'] + parameters['datfile'] + fp = open(filename, 'w') + + template = "%(totalimpacts)17d %(ncount)12d %(curyear)19.12E %(restart)s %(fracdone)9.6f %(masstot)19.12E\n" + fp.write(template % parameters) + + # Write random number seeds to the file + for index in range(parameters['seedn']): + fp.write("%12d\n" % seedarr[index]) + + fp.close() + + return + + +def copy_dists(parameters): + # Save copies of distribution files + + orig_list = ['odistribution', 'ocumulative', 'pdistribution', 'tdistribution'] + dest_list = ['odist', 'ocum', 'pdist', 'tdist'] + + for index in range(len(orig_list)): + forig = parameters['workingdir'] + orig_list[index] + '.dat' + fdest = parameters['workingdir'] + 'dist' + os.sep + dest_list[index] + "_%06d.dat" % parameters['ncount'] + shutil.copy2(forig, fdest) + + forig = parameters['workingdir'] + 'impactmass.dat' + fdest = parameters['workingdir'] + 'misc' + os.sep + "mass_%06d.dat" % parameters['ncount'] + shutil.copy2(forig, fdest) + + if (parameters['savetruelist'].upper() == 'T'): + forig = parameters['workingdir'] + 'tcumulative.dat' + fdest = parameters['workingdir'] + 'dist' + os.sep + "tcum_%06d.dat" % parameters['ncount'] + shutil.copy2(forig, fdest) + + return + + +# 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, 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 = plt.figure(figsize=(width, height), dpi=dpi) + ax = plt.axes([0, 0, 1, 1]) + ax.imshow(dem_img, interpolation="nearest", cmap='gray', vmin=0.0, vmax=1.0) + plt.axis('off') + # Save image to file + filename = parameters['workingdir'] + 'surf' + os.sep + "surf%06d.png" % parameters['ncount'] + plt.savefig(filename, dpi=dpi, bbox_inches=0) + + return + + +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'] + + ls = LightSource(azdeg=azimuth, altdeg=solar_angle) + cmap = cm.cividis + + # If min and max appear to be reversed, then fix them + if (parameters['shadedminh'] > parameters['shadedmaxh']): + temp = parameters['shadedminh'] + 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 = np.amin(DEM) + else: + shadedminh = parameters['shadedminh'] + if (parameters['shadedmaxhdefault'] == 1): + shadedmaxh = np.amax(DEM) + else: + 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 = plt.figure(figsize=(width, height), dpi=dpi) + ax = plt.axes([0, 0, 1, 1]) + ax.imshow(dem_img, interpolation="nearest", vmin=0.0, vmax=1.0) + plt.axis('off') + # Save image to file + filename = parameters['workingdir'] + 'shaded' + os.sep + "shaded%06d.png" % parameters['ncount'] + plt.savefig(filename, dpi=dpi, bbox_inches=0) + return parameters + + +def image_regolith(parameters, regolith): + # Create scaled regolith image + minref = parameters['pix'] * 1.0e-4 + maxreg = np.amax(regolith) + minreg = np.amin(regolith) + if (minreg < minref): minreg = minref + if (maxreg < minref): maxreg = (minref + 1.0e3) + regolith_scaled = np.copy(regolith) + np.place(regolith_scaled, regolith_scaled < minref, minref) + regolith_scaled = 254.0 * ( + (np.log(regolith_scaled) - np.log(minreg)) / (np.log(maxreg) - np.log(minreg))) + + # Save image to file + filename = parameters['workingdir'] + 'rego' + os.sep + "rego%06d.png" % parameters['ncount'] + height = parameters['gridsize'] / dpi + width = height + fig = plt.figure(figsize=(width, height), dpi=dpi) + fig.figimage(regolith_scaled, cmap=cm.nipy_spectral, origin='lower') + plt.savefig(filename) + + return + + +def create_rplot(parameters, odist, pdist, tdist, ph1): + # Parameters: empirical saturation limit and dfrac + satlimit = 3.12636 + dfrac = 2 ** (1. / 4) * 1.0e-3 + + # Calculate geometric saturation + minx = (parameters['pix'] / 3.0) * 1.0e-3 + maxx = 3 * parameters['pix'] * parameters['gridsize'] * 1.0e-3 + geomem = np.array([[minx, satlimit / 20.0], [maxx, satlimit / 20.0]]) + geomep = np.array([[minx, satlimit / 10.0], [maxx, satlimit / 10.0]]) + + # Create distribution arrays without zeros for plotting on log scale + idx = np.nonzero(odist[:, 5]) + odistnz = odist[idx] + odistnz = odistnz[:, [2, 5]] + + idx = np.nonzero(pdist[:, 5]) + pdistnz = pdist[idx] + pdistnz = pdistnz[:, [2, 5]] + + idx = np.nonzero(tdist[:, 5]) + tdistnz = tdist[idx] + tdistnz = tdistnz[:, [2, 5]] + + # Correct pdist + pdistnz[:, 1] = pdistnz[:, 1] * parameters['curyear'] / parameters['interval'] + + # Create sdist bin factors, which contain one crater per bin + area = (parameters['gridsize'] * parameters['pix'] * 1.0e-3) ** 2. + plo = 1 + sq2 = 2 ** (1. / 2) + while (sq2 ** plo > minx): + plo = plo - 1 + phi = plo + 1 + while (sq2 ** phi < maxx): + phi = phi + 1 + n = phi - plo + 1 + sdist = np.zeros([n, 2]) + p = plo + for index in range(n): + sdist[index, 0] = sq2 ** p + sdist[index, 1] = sq2 ** (2.0 * p + 1.5) / (area * (sq2 - 1)) + p = p + 1 + + # Create time label + tlabel = "%5.4e" % parameters['curyear'] + tlabel = tlabel.split('e') + texp = str(int(tlabel[1])) + timelabel = 'Time = ' + r'${}$ x 10$^{}$'.format(tlabel[0], texp) + ' yrs' + + # Save image to file + filename = parameters['workingdir'] + 'rplot' + os.sep + "rplot%06d.png" % parameters['ncount'] + height = parameters['gridsize'] / dpi + width = height + fig = plt.figure(figsize=(width, height), dpi=dpi) + + # Alter background color to be black, and change axis colors accordingly + plt.style.use('dark_background') + plt.rcParams['axes.prop_cycle'] + + # Plot data + plt.plot(odistnz[:, 0] * 1.0e-3, odistnz[:, 1], linewidth=3.0, color='blue') + plt.plot(pdistnz[:, 0] * 1.0e-3, pdistnz[:, 1], linewidth=2.0, linestyle='dashdot', color='white') + plt.plot(tdistnz[:, 0] * 1.0e-3, tdistnz[:, 1], linewidth=2.0, color='red') + + plt.plot(geomem[:, 0], geomem[:, 1], linewidth=2.0, linestyle=':', color='yellow') + plt.plot(geomep[:, 0], geomep[:, 1], linewidth=2.0, linestyle=':', color='yellow') + plt.plot(sdist[:, 0], sdist[:, 1], linewidth=2.0, linestyle=':', color='yellow') + + plt.plot(ph1[:, 0] * dfrac, ph1[:, 1], 'wo') + + # Create plot labels + plt.title('Crater Distribution R-Plot', fontsize=22) + plt.xlim([minx, maxx]) + plt.xscale('log') + plt.xlabel('Crater Diameter (km)', fontsize=18) + plt.ylim([5.0e-4, 5.0]) + plt.yscale('log') + plt.ylabel('R Value', fontsize=18) + plt.text(1.0e-2, 1.0, timelabel, fontsize=18) + + plt.tick_params(axis='both', which='major', labelsize=14) + plt.tick_params(axis='both', which='minor', labelsize=12) + + plt.savefig(filename) + + return From 58ee60823a62e6f78b2c7f53f372efd94ba2a99c Mon Sep 17 00:00:00 2001 From: David Minton Date: Mon, 21 Feb 2022 11:41:38 -0500 Subject: [PATCH 09/19] Fixed bug --- python/ctem/ctem/frontend.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/ctem/ctem/frontend.py b/python/ctem/ctem/frontend.py index 69327c37..f7ee95ea 100644 --- a/python/ctem/ctem/frontend.py +++ b/python/ctem/ctem/frontend.py @@ -122,4 +122,4 @@ def __init__(self, param_file="ctem.in"): return - def \ No newline at end of file + # def \ No newline at end of file From a462529ad2e2ee223f636aac0640114e8df6cb24 Mon Sep 17 00:00:00 2001 From: David Minton Date: Tue, 22 Feb 2022 08:44:06 -0500 Subject: [PATCH 10/19] Installed a nifty tool in Pycharm called Code Blocks Sorter that let me sort the methods in alphabetical order. --- python/ctem/ctem/io.py | 632 +++++++++++++++++++++-------------------- 1 file changed, 317 insertions(+), 315 deletions(-) diff --git a/python/ctem/ctem/io.py b/python/ctem/ctem/io.py index 3786d87c..d4842f8c 100644 --- a/python/ctem/ctem/io.py +++ b/python/ctem/ctem/io.py @@ -1,40 +1,289 @@ -import numpy as np -import os -import matplotlib.pyplot as plt -import shutil -from matplotlib.colors import LightSource -import matplotlib.cm as cm +def copy_dists(parameters): + # Save copies of distribution files + + orig_list = ['odistribution', 'ocumulative', 'pdistribution', 'tdistribution'] + dest_list = ['odist', 'ocum', 'pdist', 'tdist'] + + for index in range(len(orig_list)): + forig = parameters['workingdir'] + orig_list[index] + '.dat' + fdest = parameters['workingdir'] + 'dist' + os.sep + dest_list[index] + "_%06d.dat" % parameters['ncount'] + shutil.copy2(forig, fdest) + + forig = parameters['workingdir'] + 'impactmass.dat' + fdest = parameters['workingdir'] + 'misc' + os.sep + "mass_%06d.dat" % parameters['ncount'] + shutil.copy2(forig, fdest) + + if (parameters['savetruelist'].upper() == 'T'): + forig = parameters['workingdir'] + 'tcumulative.dat' + fdest = parameters['workingdir'] + 'dist' + os.sep + "tcum_%06d.dat" % parameters['ncount'] + shutil.copy2(forig, fdest) + + return -#Set pixel scaling common for image writing, at 1 pixel/ array element -dpi = 72.0 +def create_dir_structure(parameters): + # Create directories for various output files if they do not already exist + directories = ['dist', 'misc', 'rego', 'rplot', 'surf', 'shaded'] + + for directory in directories: + dir_test = parameters['workingdir'] + directory + if not os.path.isdir(dir_test): + os.makedirs(dir_test) + + return -def real2float(realstr): - """ - Converts a Fortran-generated ASCII string of a real value into a np float type. Handles cases where double precision - numbers in exponential notation use 'd' or 'D' instead of 'e' or 'E' - Parameters - ---------- - realstr : string - Fortran-generated ASCII string of a real value. +def create_rplot(parameters, odist, pdist, tdist, ph1): + # Parameters: empirical saturation limit and dfrac + satlimit = 3.12636 + dfrac = 2 ** (1. / 4) * 1.0e-3 + + # Calculate geometric saturation + minx = (parameters['pix'] / 3.0) * 1.0e-3 + maxx = 3 * parameters['pix'] * parameters['gridsize'] * 1.0e-3 + geomem = np.array([[minx, satlimit / 20.0], [maxx, satlimit / 20.0]]) + geomep = np.array([[minx, satlimit / 10.0], [maxx, satlimit / 10.0]]) + + # Create distribution arrays without zeros for plotting on log scale + idx = np.nonzero(odist[:, 5]) + odistnz = odist[idx] + odistnz = odistnz[:, [2, 5]] + + idx = np.nonzero(pdist[:, 5]) + pdistnz = pdist[idx] + pdistnz = pdistnz[:, [2, 5]] + + idx = np.nonzero(tdist[:, 5]) + tdistnz = tdist[idx] + tdistnz = tdistnz[:, [2, 5]] + + # Correct pdist + pdistnz[:, 1] = pdistnz[:, 1] * parameters['curyear'] / parameters['interval'] + + # Create sdist bin factors, which contain one crater per bin + area = (parameters['gridsize'] * parameters['pix'] * 1.0e-3) ** 2. + plo = 1 + sq2 = 2 ** (1. / 2) + while (sq2 ** plo > minx): + plo = plo - 1 + phi = plo + 1 + while (sq2 ** phi < maxx): + phi = phi + 1 + n = phi - plo + 1 + sdist = np.zeros([n, 2]) + p = plo + for index in range(n): + sdist[index, 0] = sq2 ** p + sdist[index, 1] = sq2 ** (2.0 * p + 1.5) / (area * (sq2 - 1)) + p = p + 1 + + # Create time label + tlabel = "%5.4e" % parameters['curyear'] + tlabel = tlabel.split('e') + texp = str(int(tlabel[1])) + timelabel = 'Time = ' + r'${}$ x 10$^{}$'.format(tlabel[0], texp) + ' yrs' + + # Save image to file + filename = parameters['workingdir'] + 'rplot' + os.sep + "rplot%06d.png" % parameters['ncount'] + height = parameters['gridsize'] / dpi + width = height + fig = plt.figure(figsize=(width, height), dpi=dpi) + + # Alter background color to be black, and change axis colors accordingly + plt.style.use('dark_background') + plt.rcParams['axes.prop_cycle'] + + # Plot data + plt.plot(odistnz[:, 0] * 1.0e-3, odistnz[:, 1], linewidth=3.0, color='blue') + plt.plot(pdistnz[:, 0] * 1.0e-3, pdistnz[:, 1], linewidth=2.0, linestyle='dashdot', color='white') + plt.plot(tdistnz[:, 0] * 1.0e-3, tdistnz[:, 1], linewidth=2.0, color='red') + + plt.plot(geomem[:, 0], geomem[:, 1], linewidth=2.0, linestyle=':', color='yellow') + plt.plot(geomep[:, 0], geomep[:, 1], linewidth=2.0, linestyle=':', color='yellow') + plt.plot(sdist[:, 0], sdist[:, 1], linewidth=2.0, linestyle=':', color='yellow') + + plt.plot(ph1[:, 0] * dfrac, ph1[:, 1], 'wo') + + # Create plot labels + plt.title('Crater Distribution R-Plot', fontsize=22) + plt.xlim([minx, maxx]) + plt.xscale('log') + plt.xlabel('Crater Diameter (km)', fontsize=18) + plt.ylim([5.0e-4, 5.0]) + plt.yscale('log') + plt.ylabel('R Value', fontsize=18) + plt.text(1.0e-2, 1.0, timelabel, fontsize=18) + + plt.tick_params(axis='both', which='major', labelsize=14) + plt.tick_params(axis='both', which='minor', labelsize=12) + + plt.savefig(filename) + + return + +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 = plt.figure(figsize=(width, height), dpi=dpi) + ax = plt.axes([0, 0, 1, 1]) + ax.imshow(dem_img, interpolation="nearest", cmap='gray', vmin=0.0, vmax=1.0) + plt.axis('off') + # Save image to file + filename = parameters['workingdir'] + 'surf' + os.sep + "surf%06d.png" % parameters['ncount'] + plt.savefig(filename, dpi=dpi, bbox_inches=0) + + return + + +def image_regolith(parameters, regolith): + # Create scaled regolith image + minref = parameters['pix'] * 1.0e-4 + maxreg = np.amax(regolith) + minreg = np.amin(regolith) + if (minreg < minref): minreg = minref + if (maxreg < minref): maxreg = (minref + 1.0e3) + regolith_scaled = np.copy(regolith) + np.place(regolith_scaled, regolith_scaled < minref, minref) + regolith_scaled = 254.0 * ( + (np.log(regolith_scaled) - np.log(minreg)) / (np.log(maxreg) - np.log(minreg))) + + # Save image to file + filename = parameters['workingdir'] + 'rego' + os.sep + "rego%06d.png" % parameters['ncount'] + height = parameters['gridsize'] / dpi + width = height + fig = plt.figure(figsize=(width, height), dpi=dpi) + fig.figimage(regolith_scaled, cmap=cm.nipy_spectral, origin='lower') + plt.savefig(filename) + + return - Returns - ------- - : float - The converted floating point value of the input string - """ - return float(realstr.replace('d', 'E').replace('D', 'E')) +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'] + + ls = LightSource(azdeg=azimuth, altdeg=solar_angle) + cmap = cm.cividis + + # If min and max appear to be reversed, then fix them + if (parameters['shadedminh'] > parameters['shadedmaxh']): + temp = parameters['shadedminh'] + 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 = np.amin(DEM) + else: + shadedminh = parameters['shadedminh'] + if (parameters['shadedmaxhdefault'] == 1): + shadedmaxh = np.amax(DEM) + else: + 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 = plt.figure(figsize=(width, height), dpi=dpi) + ax = plt.axes([0, 0, 1, 1]) + ax.imshow(dem_img, interpolation="nearest", vmin=0.0, vmax=1.0) + plt.axis('off') + # Save image to file + filename = parameters['workingdir'] + 'shaded' + os.sep + "shaded%06d.png" % parameters['ncount'] + plt.savefig(filename, dpi=dpi, bbox_inches=0) + return parameters + + +def read_ctemdat(parameters, seedarr): + # Read and parse ctem.dat file + datfile = parameters['workingdir'] + 'ctem.dat' + + # Read ctem.dat file + print('Reading input file ' + parameters['datfile']) + fp = open(datfile, 'r') + lines = fp.readlines() + fp.close() + + # Parse file lines and update parameter fields + fields = lines[0].split() + if len(fields) > 0: + parameters['totalimpacts'] = real2float(fields[0]) + parameters['ncount'] = int(fields[1]) + parameters['curyear'] = real2float(fields[2]) + parameters['restart'] = fields[3] + parameters['fracdone'] = real2float(fields[4]) + parameters['masstot'] = real2float(fields[5]) + + # Parse remainder of file to build seed array + nlines = len(lines) + index = 1 + while (index < nlines): + fields = lines[index].split() + seedarr[index - 1] = real2float(fields[0]) + index += 1 + + parameters['seedn'] = index - 1 + + return + + +def read_formatted_ascii(filename, skip_lines): + # Generalized ascii text reader + # For use with sfdcompare, production, odist, tdist, pdist data files + data = np.genfromtxt(filename, skip_header=skip_lines) + return data + + +def read_impact_mass(filename): + # Read impact mass file + + fp = open(filename, 'r') + line = fp.readlines() + fp.close() + + fields = line[0].split() + if (len(fields) > 0): + mass = real2float(fields[0]) + else: + mass = 0 + + return mass + + +# Write production function to file production.dat +# This file format does not exactly match that generated from IDL. Does it work? def read_param(parameters): # Read and parse ctem.in file inputfile = os.path.join(parameters['workingdir'], parameters['ctemfile']) - + # Read ctem.in file print('Reading input file ' + parameters['ctemfile']) with open(inputfile, 'r') as fp: lines = fp.readlines() - + # Process file text for line in lines: fields = line.split() @@ -61,7 +310,7 @@ def read_param(parameters): if ('shadedmaxh' == fields[0].lower()): parameters['shadedmaxh'] = real2float(fields[1]) parameters['shadedmaxhdefault'] = 0 - + # Test values for further processing if (parameters['interval'] <= 0.0): print('Invalid value for or missing variable INTERVAL in ' + inputfile) @@ -85,101 +334,42 @@ def read_param(parameters): print('Invalid value for or missing variable SAVEREGO in ' + inputfile) if (parameters['savepres'] is None): print('Invalid value for or missing variable SAVEPRES in ' + inputfile) - if (parameters['savetruelist'] is None): - print('Invalid value for or missing variable SAVETRUELIST in ' + inputfile) - if (parameters['runtype'] is None): - print('Invalid value for or missing variable RUNTYPE in ' + inputfile) - if (parameters['restart'] is None): - print('Invalid value for or missing variable RESTART in ' + inputfile) - - return parameters - - -def read_formatted_ascii(filename, skip_lines): - # Generalized ascii text reader - # For use with sfdcompare, production, odist, tdist, pdist data files - data = np.genfromtxt(filename, skip_header=skip_lines) - return data - - -def read_unformatted_binary(filename, gridsize): - # Read unformatted binary files created by Fortran - # For use with surface ejecta and surface dem data files - dt = np.float - data = np.fromfile(filename, dtype=dt) - data.shape = (gridsize, gridsize) - - return data - - -def read_ctemdat(parameters, seedarr): - # Read and parse ctem.dat file - datfile = parameters['workingdir'] + 'ctem.dat' - - # Read ctem.dat file - print('Reading input file ' + parameters['datfile']) - fp = open(datfile, 'r') - lines = fp.readlines() - fp.close() - - # Parse file lines and update parameter fields - fields = lines[0].split() - if len(fields) > 0: - parameters['totalimpacts'] = real2float(fields[0]) - parameters['ncount'] = int(fields[1]) - parameters['curyear'] = real2float(fields[2]) - parameters['restart'] = fields[3] - parameters['fracdone'] = real2float(fields[4]) - parameters['masstot'] = real2float(fields[5]) - - # Parse remainder of file to build seed array - nlines = len(lines) - index = 1 - while (index < nlines): - fields = lines[index].split() - seedarr[index - 1] = real2float(fields[0]) - index += 1 - - parameters['seedn'] = index - 1 - - return - - -def read_impact_mass(filename): - # Read impact mass file - - fp = open(filename, 'r') - line = fp.readlines() - fp.close() - - fields = line[0].split() - if (len(fields) > 0): - mass = real2float(fields[0]) - else: - mass = 0 + if (parameters['savetruelist'] is None): + print('Invalid value for or missing variable SAVETRUELIST in ' + inputfile) + if (parameters['runtype'] is None): + print('Invalid value for or missing variable RUNTYPE in ' + inputfile) + if (parameters['restart'] is None): + print('Invalid value for or missing variable RESTART in ' + inputfile) - return mass + return parameters -# Write production function to file production.dat -# This file format does not exactly match that generated from IDL. Does it work? -def write_production(parameters, production): - filename = parameters['workingdir'] + parameters['sfdfile'] - np.savetxt(filename, production, fmt='%1.8e', delimiter=' ') +def read_unformatted_binary(filename, gridsize): + # Read unformatted binary files created by Fortran + # For use with surface ejecta and surface dem data files + dt = np.float + data = np.fromfile(filename, dtype=dt) + data.shape = (gridsize, gridsize) - return + return data -def create_dir_structure(parameters): - # Create directories for various output files if they do not already exist - directories = ['dist', 'misc', 'rego', 'rplot', 'surf', 'shaded'] - - for directory in directories: - dir_test = parameters['workingdir'] + directory - if not os.path.isdir(dir_test): - os.makedirs(dir_test) - - return +def real2float(realstr): + """ + Converts a Fortran-generated ASCII string of a real value into a np float type. Handles cases where double precision + numbers in exponential notation use 'd' or 'D' instead of 'e' or 'E' + + Parameters + ---------- + realstr : string + Fortran-generated ASCII string of a real value. + + Returns + ------- + : float + The converted floating point value of the input string + """ + return float(realstr.replace('d', 'E').replace('D', 'E')) def write_ctemdat(parameters, seedarr): @@ -189,7 +379,7 @@ def write_ctemdat(parameters, seedarr): template = "%(totalimpacts)17d %(ncount)12d %(curyear)19.12E %(restart)s %(fracdone)9.6f %(masstot)19.12E\n" fp.write(template % parameters) - + # Write random number seeds to the file for index in range(parameters['seedn']): fp.write("%12d\n" % seedarr[index]) @@ -199,212 +389,24 @@ def write_ctemdat(parameters, seedarr): return -def copy_dists(parameters): - # Save copies of distribution files - - orig_list = ['odistribution', 'ocumulative', 'pdistribution', 'tdistribution'] - dest_list = ['odist', 'ocum', 'pdist', 'tdist'] - - for index in range(len(orig_list)): - forig = parameters['workingdir'] + orig_list[index] + '.dat' - fdest = parameters['workingdir'] + 'dist' + os.sep + dest_list[index] + "_%06d.dat" % parameters['ncount'] - shutil.copy2(forig, fdest) - - forig = parameters['workingdir'] + 'impactmass.dat' - fdest = parameters['workingdir'] + 'misc' + os.sep + "mass_%06d.dat" % parameters['ncount'] - shutil.copy2(forig, fdest) - - if (parameters['savetruelist'].upper() == 'T'): - forig = parameters['workingdir'] + 'tcumulative.dat' - fdest = parameters['workingdir'] + 'dist' + os.sep + "tcum_%06d.dat" % parameters['ncount'] - shutil.copy2(forig, fdest) - - return - - # 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, 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 = plt.figure(figsize=(width, height), dpi=dpi) - ax = plt.axes([0, 0, 1, 1]) - ax.imshow(dem_img, interpolation="nearest", cmap='gray', vmin=0.0, vmax=1.0) - plt.axis('off') - # Save image to file - filename = parameters['workingdir'] + 'surf' + os.sep + "surf%06d.png" % parameters['ncount'] - plt.savefig(filename, dpi=dpi, bbox_inches=0) - - return - - -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'] - - ls = LightSource(azdeg=azimuth, altdeg=solar_angle) - cmap = cm.cividis - - # If min and max appear to be reversed, then fix them - if (parameters['shadedminh'] > parameters['shadedmaxh']): - temp = parameters['shadedminh'] - 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 = np.amin(DEM) - else: - shadedminh = parameters['shadedminh'] - if (parameters['shadedmaxhdefault'] == 1): - shadedmaxh = np.amax(DEM) - else: - 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 = plt.figure(figsize=(width, height), dpi=dpi) - ax = plt.axes([0, 0, 1, 1]) - ax.imshow(dem_img, interpolation="nearest", vmin=0.0, vmax=1.0) - plt.axis('off') - # Save image to file - filename = parameters['workingdir'] + 'shaded' + os.sep + "shaded%06d.png" % parameters['ncount'] - plt.savefig(filename, dpi=dpi, bbox_inches=0) - return parameters +# http://nbviewer.jupyter.org/github/ThomasLecocq/geophysique.be/blob/master/2014-02-25%20Shaded%20Relief%20Map%20in%20Python.ipynb - -def image_regolith(parameters, regolith): - # Create scaled regolith image - minref = parameters['pix'] * 1.0e-4 - maxreg = np.amax(regolith) - minreg = np.amin(regolith) - if (minreg < minref): minreg = minref - if (maxreg < minref): maxreg = (minref + 1.0e3) - regolith_scaled = np.copy(regolith) - np.place(regolith_scaled, regolith_scaled < minref, minref) - regolith_scaled = 254.0 * ( - (np.log(regolith_scaled) - np.log(minreg)) / (np.log(maxreg) - np.log(minreg))) - - # Save image to file - filename = parameters['workingdir'] + 'rego' + os.sep + "rego%06d.png" % parameters['ncount'] - height = parameters['gridsize'] / dpi - width = height - fig = plt.figure(figsize=(width, height), dpi=dpi) - fig.figimage(regolith_scaled, cmap=cm.nipy_spectral, origin='lower') - plt.savefig(filename) +def write_production(parameters, production): + filename = parameters['workingdir'] + parameters['sfdfile'] + np.savetxt(filename, production, fmt='%1.8e', delimiter=' ') return -def create_rplot(parameters, odist, pdist, tdist, ph1): - # Parameters: empirical saturation limit and dfrac - satlimit = 3.12636 - dfrac = 2 ** (1. / 4) * 1.0e-3 - - # Calculate geometric saturation - minx = (parameters['pix'] / 3.0) * 1.0e-3 - maxx = 3 * parameters['pix'] * parameters['gridsize'] * 1.0e-3 - geomem = np.array([[minx, satlimit / 20.0], [maxx, satlimit / 20.0]]) - geomep = np.array([[minx, satlimit / 10.0], [maxx, satlimit / 10.0]]) - - # Create distribution arrays without zeros for plotting on log scale - idx = np.nonzero(odist[:, 5]) - odistnz = odist[idx] - odistnz = odistnz[:, [2, 5]] - - idx = np.nonzero(pdist[:, 5]) - pdistnz = pdist[idx] - pdistnz = pdistnz[:, [2, 5]] - - idx = np.nonzero(tdist[:, 5]) - tdistnz = tdist[idx] - tdistnz = tdistnz[:, [2, 5]] - - # Correct pdist - pdistnz[:, 1] = pdistnz[:, 1] * parameters['curyear'] / parameters['interval'] - - # Create sdist bin factors, which contain one crater per bin - area = (parameters['gridsize'] * parameters['pix'] * 1.0e-3) ** 2. - plo = 1 - sq2 = 2 ** (1. / 2) - while (sq2 ** plo > minx): - plo = plo - 1 - phi = plo + 1 - while (sq2 ** phi < maxx): - phi = phi + 1 - n = phi - plo + 1 - sdist = np.zeros([n, 2]) - p = plo - for index in range(n): - sdist[index, 0] = sq2 ** p - sdist[index, 1] = sq2 ** (2.0 * p + 1.5) / (area * (sq2 - 1)) - p = p + 1 - - # Create time label - tlabel = "%5.4e" % parameters['curyear'] - tlabel = tlabel.split('e') - texp = str(int(tlabel[1])) - timelabel = 'Time = ' + r'${}$ x 10$^{}$'.format(tlabel[0], texp) + ' yrs' - - # Save image to file - filename = parameters['workingdir'] + 'rplot' + os.sep + "rplot%06d.png" % parameters['ncount'] - height = parameters['gridsize'] / dpi - width = height - fig = plt.figure(figsize=(width, height), dpi=dpi) +dpi = 72.0 - # Alter background color to be black, and change axis colors accordingly - plt.style.use('dark_background') - plt.rcParams['axes.prop_cycle'] - # Plot data - plt.plot(odistnz[:, 0] * 1.0e-3, odistnz[:, 1], linewidth=3.0, color='blue') - plt.plot(pdistnz[:, 0] * 1.0e-3, pdistnz[:, 1], linewidth=2.0, linestyle='dashdot', color='white') - plt.plot(tdistnz[:, 0] * 1.0e-3, tdistnz[:, 1], linewidth=2.0, color='red') - - plt.plot(geomem[:, 0], geomem[:, 1], linewidth=2.0, linestyle=':', color='yellow') - plt.plot(geomep[:, 0], geomep[:, 1], linewidth=2.0, linestyle=':', color='yellow') - plt.plot(sdist[:, 0], sdist[:, 1], linewidth=2.0, linestyle=':', color='yellow') - - plt.plot(ph1[:, 0] * dfrac, ph1[:, 1], 'wo') +from matplotlib.colors import LightSource +import matplotlib.cm as cm - # Create plot labels - plt.title('Crater Distribution R-Plot', fontsize=22) - plt.xlim([minx, maxx]) - plt.xscale('log') - plt.xlabel('Crater Diameter (km)', fontsize=18) - plt.ylim([5.0e-4, 5.0]) - plt.yscale('log') - plt.ylabel('R Value', fontsize=18) - plt.text(1.0e-2, 1.0, timelabel, fontsize=18) - - plt.tick_params(axis='both', which='major', labelsize=14) - plt.tick_params(axis='both', which='minor', labelsize=12) - - plt.savefig(filename) - - return +# Set pixel scaling common for image writing, at 1 pixel/ array element +import matplotlib.pyplot as plt +import numpy as np +import os +import shutil From 56bc4c40995b046c79ec7a6e60304d674c7fc20d Mon Sep 17 00:00:00 2001 From: David Minton Date: Tue, 22 Feb 2022 09:52:12 -0500 Subject: [PATCH 11/19] Fixed position of import statements --- python/ctem/ctem/io.py | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/python/ctem/ctem/io.py b/python/ctem/ctem/io.py index d4842f8c..60ae499e 100644 --- a/python/ctem/ctem/io.py +++ b/python/ctem/ctem/io.py @@ -1,3 +1,13 @@ +import numpy as np +import os +import shutil +from matplotlib.colors import LightSource +import matplotlib.cm as cm +import matplotlib.pyplot as plt + +# Set pixel scaling common for image writing, at 1 pixel/ array element +dpi = 72.0 + def copy_dists(parameters): # Save copies of distribution files @@ -398,15 +408,3 @@ def write_production(parameters, production): return - -dpi = 72.0 - - -from matplotlib.colors import LightSource -import matplotlib.cm as cm - -# Set pixel scaling common for image writing, at 1 pixel/ array element -import matplotlib.pyplot as plt -import numpy as np -import os -import shutil From cee986627317452a41fe8bb383a00effce5e5457 Mon Sep 17 00:00:00 2001 From: David Minton Date: Tue, 22 Feb 2022 11:54:27 -0500 Subject: [PATCH 12/19] Restructured the driver to be object oriented and streamlined --- python/ctem/ctem/__init__.py | 2 +- python/ctem/ctem/ctem_draw_surf.py | 69 ----- python/ctem/ctem/ctem_driver.py | 229 --------------- python/ctem/ctem/ctem_io_readers.py | 163 ----------- python/ctem/ctem/ctem_io_writers.py | 264 ------------------ python/ctem/ctem/dist_reader.py | 29 -- python/ctem/ctem/driver.py | 263 +++++++++++++++++ python/ctem/ctem/frontend.py | 125 --------- python/ctem/ctem/io.py | 205 +++++++------- .../ctem/{ctem_viewer_3d.py => viewer3d.py} | 2 +- python/ctem/tests/testio/testio.py | 2 +- 11 files changed, 367 insertions(+), 986 deletions(-) delete mode 100644 python/ctem/ctem/ctem_draw_surf.py delete mode 100644 python/ctem/ctem/ctem_driver.py delete mode 100644 python/ctem/ctem/ctem_io_readers.py delete mode 100644 python/ctem/ctem/ctem_io_writers.py delete mode 100644 python/ctem/ctem/dist_reader.py create mode 100644 python/ctem/ctem/driver.py delete mode 100644 python/ctem/ctem/frontend.py rename python/ctem/ctem/{ctem_viewer_3d.py => viewer3d.py} (98%) diff --git a/python/ctem/ctem/__init__.py b/python/ctem/ctem/__init__.py index b9aa95ac..165a1bb0 100644 --- a/python/ctem/ctem/__init__.py +++ b/python/ctem/ctem/__init__.py @@ -1,4 +1,4 @@ # from ctem.ctem_io_readers import * # from ctem.ctem_io_writers import * # from ctem.ctem_driver import * -from ctem.frontend import * \ No newline at end of file +from ctem.driver import * \ No newline at end of file diff --git a/python/ctem/ctem/ctem_draw_surf.py b/python/ctem/ctem/ctem_draw_surf.py deleted file mode 100644 index 74bf81e1..00000000 --- a/python/ctem/ctem/ctem_draw_surf.py +++ /dev/null @@ -1,69 +0,0 @@ -#!/usr/local/bin/python -# -#Cratered Terrain Evolution Model driver -# -#Original IDL design: Jim Richardson, Arecibo Observatory -#Revised IDL design: David Minton, Purdue University -#Re-engineered design and Python implementation: Matthew Route, Purdue University -#August 2016 - -#Import general purpose modules -import numpy -import os -import subprocess -import shutil - -#Import CTEM modules -import ctem_io_readers -import ctem_io_writers - -#Create and initialize data dictionaries for parameters and options from CTEM.in -notset = '-NOTSET-' -currentdir = os.getcwd() + os.sep - -parameters={'restart': notset, - 'runtype': notset, - 'popupconsole': notset, - 'saveshaded': notset, - 'saverego': notset, - 'savepres': notset, - 'savetruelist': notset, - 'seedn': 1, - 'totalimpacts': 0, - 'ncount': 0, - 'curyear': 0.0, - 'fracdone': 1.0, - 'masstot': 0.0, - 'interval': 0.0, - 'numintervals': 0, - 'pix': -1.0, - 'gridsize': -1, - 'seed': 0, - 'maxcrat': 1.0, - 'shadedminhdefault': 1, - 'shadedmaxhdefault': 1, - 'shadedminh': 0.0, - 'shadedmaxh': 0.0, - 'workingdir': currentdir, - 'ctemfile': 'ctem.in', - 'datfile': 'ctem.dat', - 'impfile': notset, - 'sfdcompare': notset, - 'sfdfile': notset} - -#Read ctem.in to initialize parameter values based on user input -ctem_io_readers.read_ctemin(parameters,notset) -ctem_io_writers.create_dir_structure(parameters) - -dem_file = parameters['workingdir'] + 'surface_dem.dat' -surface_dem = ctem_io_readers.read_unformatted_binary(dem_file, parameters['gridsize']) - -#Set up data arrays -seedarr = numpy.zeros(100, dtype = numpy.int) -seedarr[0] = parameters['seed'] - -#Read ctem.dat file -ctem_io_readers.read_ctemdat(parameters, seedarr) - -#Write surface dem -ctem_io_writers.image_dem(parameters, surface_dem) \ No newline at end of file diff --git a/python/ctem/ctem/ctem_driver.py b/python/ctem/ctem/ctem_driver.py deleted file mode 100644 index 29057053..00000000 --- a/python/ctem/ctem/ctem_driver.py +++ /dev/null @@ -1,229 +0,0 @@ -#!/usr/local/bin/python -# -#Cratered Terrain Evolution Model driver -# -#Original IDL design: Jim Richardson, Arecibo Observatory -#Revised IDL design: David Minton, Purdue University -#Re-engineered design and Python implementation: Matthew Route, Purdue University -#August 2016 - -#Import general purpose modules -import numpy -import os -import subprocess -import shutil - -#Import CTEM modules -import ctem_io_readers -import ctem_io_writers - -#Create and initialize data dictionaries for parameters and options from CTEM.in -notset = '-NOTSET-' -currentdir = os.getcwd() + os.sep - -parameters={'restart': notset, - 'runtype': notset, - 'popupconsole': notset, - 'saveshaded': notset, - 'saverego': notset, - 'savepres': notset, - 'savetruelist': notset, - 'seedn': 1, - 'totalimpacts': 0, - 'ncount': 0, - 'curyear': 0.0, - 'fracdone': 1.0, - 'masstot': 0.0, - 'interval': 0.0, - 'numintervals': 0, - 'pix': -1.0, - 'gridsize': -1, - 'seed': 0, - 'maxcrat': 1.0, - 'shadedminhdefault': 1, - 'shadedmaxhdefault': 1, - 'shadedminh': 0.0, - 'shadedmaxh': 0.0, - 'workingdir': currentdir, - 'ctemfile': 'ctem.in', - 'datfile': 'ctem.dat', - 'impfile': notset, - 'sfdcompare': notset, - 'sfdfile': notset} - -#Read ctem.in to initialize parameter values based on user input -ctem_io_readers.read_ctemin(parameters,notset) - -#Read sfdcompare file -sfdfile = parameters['workingdir'] + parameters['sfdcompare'] -ph1 = ctem_io_readers.read_formatted_ascii(sfdfile, skip_lines = 0) - -#Set up data arrays -seedarr = numpy.zeros(100, dtype = numpy.int) -seedarr[0] = parameters['seed'] -odist = numpy.zeros([1, 6]) -pdist = numpy.zeros([1, 6]) -tdist = numpy.zeros([1, 6]) -surface_dem = numpy.zeros([parameters['gridsize'], parameters['gridsize']], dtype = numpy.float) -regolith = numpy.zeros([parameters['gridsize'], parameters['gridsize']], dtype =numpy.float) - -#Read production function file -impfile = parameters['workingdir'] + parameters['impfile'] -prodfunction = ctem_io_readers.read_formatted_ascii(impfile, skip_lines = 0) - -#Create impactor production population -area = (parameters['gridsize'] * parameters['pix'])**2 -production = numpy.copy(prodfunction) -production[:,1] = production[:,1] * area * parameters['interval'] - -#Write corrected production function to file -ctem_io_writers.write_production(parameters, production) - -#Starting new or old run? -if (parameters['restart'].upper() == 'F'): - print('Starting a new run') - - if (parameters['runtype'].upper() == 'STATISTICAL'): - parameters['ncount'] = 1 - - #Write ctem.dat file - ctem_io_writers.write_ctemdat(parameters, seedarr) - - else: - parameters['ncount'] = 0 - - #Delete tdistribution file, if it exists - tdist_file = parameters['workingdir'] + 'tdistribution.dat' - if os.path.isfile(tdist_file): - os.remove(tdist_file) - -else: - print('Continuing a previous run') - - #Read surface dem(shaded relief) and ejecta data files - dem_file = parameters['workingdir'] + 'surface_dem.dat' - surface_dem = ctem_io_readers.read_unformatted_binary(dem_file, parameters['gridsize']) - ejecta_file = parameters['workingdir'] + 'surface_ejc.dat' - regolith = ctem_io_readers.read_unformatted_binary(ejecta_file, parameters['gridsize']) - - #Read odistribution, tdistribution, and pdistribution files - ofile = parameters['workingdir'] + 'odistribution.dat' - odist = ctem_io_readers.read_formatted_ascii(ofile, skip_lines = 1) - tfile = parameters['workingdir'] + 'tdistribution.dat' - tdist = ctem_io_readers.read_formatted_ascii(tfile, skip_lines = 1) - pfile = parameters['workingdir'] + 'pdistribution.dat' - pdist = ctem_io_readers.read_formatted_ascii(pfile, skip_lines = 1) - - #Read impact mass from file - massfile = parameters['workingdir'] + 'impactmass.dat' - impact_mass = ctem_io_readers.read_impact_mass(massfile) - - #Read ctem.dat file - ctem_io_readers.read_ctemdat(parameters, seedarr) - -#Open fracdonefile and regodepthfile for writing -filename = parameters['workingdir'] + 'fracdone.dat' -fp_frac = open(filename,'w') -filename = parameters['workingdir'] + 'regolithdepth.dat' -fp_reg = open(filename,'w') - -#Begin CTEM processing loops -print('Beginning loops') - -ctem_io_writers.create_dir_structure(parameters) - -while (parameters['ncount'] <= parameters['numintervals']): - - #Create crater population - if (parameters['ncount'] > 0): - - #Move ctem.dat - forig = parameters['workingdir'] + 'ctem.dat' - fdest = parameters['workingdir'] + 'misc' + os.sep + "ctem%06d.dat" % parameters['ncount'] - shutil.copy2(forig, fdest) - - #Create crater population and display CTEM progress on screen - print(parameters['ncount'], ' Calling FORTRAN routine') - with subprocess.Popen([parameters['workingdir']+'CTEM'], - stdout=subprocess.PIPE, - universal_newlines=True) as p: - for line in p.stdout: - print(line, end='') - - - #Optional: do not pipe CTEM progress to the screen - #subprocess.check_output([parameters['workingdir']+'CTEM']) - - #Read Fortran output - print(parameters['ncount'], ' Reading Fortran output') - - #Read surface dem(shaded relief) and ejecta data files - dem_file = parameters['workingdir'] + 'surface_dem.dat' - surface_dem = ctem_io_readers.read_unformatted_binary(dem_file, parameters['gridsize']) - ejecta_file = parameters['workingdir'] + 'surface_ejc.dat' - regolith = ctem_io_readers.read_unformatted_binary(ejecta_file, parameters['gridsize']) - - #Read odistribution, tdistribution, and pdistribution files - ofile = parameters['workingdir'] + 'odistribution.dat' - odist = ctem_io_readers.read_formatted_ascii(ofile, skip_lines = 1) - tfile = parameters['workingdir'] + 'tdistribution.dat' - tdist = ctem_io_readers.read_formatted_ascii(tfile, skip_lines = 1) - pfile = parameters['workingdir'] + 'pdistribution.dat' - pdist = ctem_io_readers.read_formatted_ascii(pfile, skip_lines = 1) - - #Read impact mass from file - massfile = parameters['workingdir'] + 'impactmass.dat' - impact_mass = ctem_io_readers.read_impact_mass(massfile) - - #Read ctem.dat file - ctem_io_readers.read_ctemdat(parameters, seedarr) - - #Update parameters: mass, curyear, regolith properties - parameters['masstot'] = parameters['masstot'] + impact_mass - - parameters['curyear'] = parameters['curyear'] + parameters['fracdone'] * parameters['interval'] - template = "%(fracdone)9.6f %(curyear)19.12E\n" - fp_frac.write(template % parameters) - - reg_text = "%19.12E %19.12E %19.12E %19.12E\n" % (parameters['curyear'], - numpy.mean(regolith), numpy.amax(regolith), numpy.amin(regolith)) - fp_reg.write(reg_text) - - #Save copy of crater distribution files - ctem_io_writers.copy_dists(parameters) - - #Display results - print(parameters['ncount'], ' Displaying results') - - #Write surface dem, surface ejecta, shaded relief, and rplot data - ctem_io_writers.image_dem(parameters, surface_dem) - if (parameters['saverego'].upper() == 'T'): - ctem_io_writers.image_regolith(parameters, regolith) - if (parameters['saveshaded'].upper() == 'T'): - ctem_io_writers.image_shaded_relief(parameters, surface_dem) - if (parameters['savepres'].upper() == 'T'): - ctem_io_writers.create_rplot(parameters,odist,pdist,tdist,ph1) - - #Update ncount - parameters['ncount'] = parameters['ncount'] + 1 - - if ((parameters['runtype'].upper() == 'STATISTICAL') or (parameters['ncount'] == 1)): - parameters['restart'] = 'F' - parameters['curyear'] = 0.0 - parameters['totalimpacts'] = 0 - parameters['masstot'] = 0.0 - - #Delete tdistribution file, if it exists - tdist_file = parameters['workingdir'] + 'tdistribution.dat' - if os.path.isfile(tdist_file): - os.remove(tdist_file) - - else: - parameters['restart'] = 'T' - - #Write ctem.dat file - ctem_io_writers.write_ctemdat(parameters, seedarr) - -#Close updateable fracdonefile and regodepthfile files -fp_frac.close() -fp_reg.close() diff --git a/python/ctem/ctem/ctem_io_readers.py b/python/ctem/ctem/ctem_io_readers.py deleted file mode 100644 index 0935fac9..00000000 --- a/python/ctem/ctem/ctem_io_readers.py +++ /dev/null @@ -1,163 +0,0 @@ -#!/usr/local/bin/python -# -#Cratered Terrain Evolution Model file reading utilities -# -#Original IDL design: Jim Richardson, Arecibo Observatory -#Revised IDL design: David Minton, Purdue University -#Re-engineered design and Python implementation: Matthew Route, Purdue University -#August 2016 -# -#Known issues for operation with CTEM -#1) ctem.in has 1.0d0 value for maxcrat which is not readable by Python - -import numpy - - -def real2float(realstr): - """ - Converts a Fortran-generated ASCII string of a real value into a numpy float type. Handles cases where double precision - numbers in exponential notation use 'd' or 'D' instead of 'e' or 'E' - - Parameters - ---------- - realstr : string - Fortran-generated ASCII string of a real value. - - Returns - ------- - : float - The converted floating point value of the input string - """ - return float(realstr.replace('d', 'E').replace('D', 'E')) - - -def read_ctemin(parameters,notset): - #Read and parse ctem.in file - inputfile = parameters['workingdir'] + parameters['ctemfile'] - - #Read ctem.in file - print('Reading input file '+ parameters['ctemfile']) - fp = open(inputfile,'r') - lines = fp.readlines() - fp.close() - - #Process file text - for line in lines: - fields = line.split() - if len(fields) > 0: - if ('pix' == fields[0].lower()): parameters['pix']=real2float(fields[1]) - if ('gridsize' == fields[0].lower()): parameters['gridsize']=int(fields[1]) - if ('seed' == fields[0].lower()): parameters['seed']=int(fields[1]) - if ('sfdfile' == fields[0].lower()): parameters['sfdfile']=fields[1] - if ('impfile' == fields[0].lower()): parameters['impfile']=fields[1] - if ('maxcrat' == fields[0].lower()): parameters['maxcrat']=real2float(fields[1]) - if ('sfdcompare' == fields[0].lower()): parameters['sfdcompare']=fields[1] - if ('interval' == fields[0].lower()): parameters['interval']=real2float(fields[1]) - if ('numintervals' == fields[0].lower()): parameters['numintervals']=int(fields[1]) - if ('popupconsole' == fields[0].lower()): parameters['popupconsole']=fields[1] - if ('saveshaded' == fields[0].lower()): parameters['saveshaded']=fields[1] - if ('saverego' == fields[0].lower()): parameters['saverego']=fields[1] - if ('savepres' == fields[0].lower()): parameters['savepres']=fields[1] - if ('savetruelist' == fields[0].lower()): parameters['savetruelist']=fields[1] - if ('runtype' == fields[0].lower()): parameters['runtype']=fields[1] - if ('restart' == fields[0].lower()): parameters['restart']=fields[1] - if ('shadedminh' == fields[0].lower()): - parameters['shadedminh'] = real2float(fields[1]) - parameters['shadedminhdefault'] = 0 - if ('shadedmaxh' == fields[0].lower()): - parameters['shadedmaxh'] = real2float(fields[1]) - parameters['shadedmaxhdefault'] = 0 - - #Test values for further processing - if (parameters['interval'] <= 0.0): - print('Invalid value for or missing variable INTERVAL in '+ inputfile) - if (parameters['numintervals'] <= 0): - print('Invalid value for or missing variable NUMINTERVALS in '+ inputfile) - if (parameters['pix'] <= 0.0): - print('Invalid value for or missing variable PIX in '+ inputfile) - if (parameters['gridsize'] <= 0): - print('Invalid value for or missing variable GRIDSIZE in '+ inputfile) - if (parameters['seed'] == 0): - print('Invalid value for or missing variable SEED in '+ inputfile) - if (parameters['sfdfile'] == notset): - print('Invalid value for or missing variable SFDFILE in '+ inputfile) - if (parameters['impfile'] == notset): - print('Invalid value for or missing variable IMPFILE in '+ inputfile) - if (parameters['popupconsole'] == notset): - print('Invalid value for or missing variable POPUPCONSOLE in '+ inputfile) - if (parameters['saveshaded'] == notset): - print('Invalid value for or missing variable SAVESHADED in '+ inputfile) - if (parameters['saverego'] == notset): - print('Invalid value for or missing variable SAVEREGO in '+ inputfile) - if (parameters['savepres'] == notset): - print('Invalid value for or missing variable SAVEPRES in '+ inputfile) - if (parameters['savetruelist'] == notset): - print('Invalid value for or missing variable SAVETRUELIST in '+ inputfile) - if (parameters['runtype'] == notset): - print('Invalid value for or missing variable RUNTYPE in '+ inputfile) - if (parameters['restart'] == notset): - print('Invalid value for or missing variable RESTART in '+ inputfile) - - return - -def read_formatted_ascii(filename, skip_lines): - #Generalized ascii text reader - #For use with sfdcompare, production, odist, tdist, pdist data files - data = numpy.genfromtxt(filename, skip_header = skip_lines) - return data - -def read_unformatted_binary(filename, gridsize): - #Read unformatted binary files created by Fortran - #For use with surface ejecta and surface dem data files - dt = numpy.float - data = numpy.fromfile(filename, dtype = dt) - data.shape = (gridsize,gridsize) - - return data - -def read_ctemdat(parameters, seedarr): - #Read and parse ctem.dat file - datfile = parameters['workingdir'] + 'ctem.dat' - - #Read ctem.dat file - print('Reading input file '+ parameters['datfile']) - fp = open(datfile,'r') - lines = fp.readlines() - fp.close() - - #Parse file lines and update parameter fields - fields = lines[0].split() - if len(fields) > 0: - parameters['totalimpacts'] = real2float(fields[0]) - parameters['ncount'] = int(fields[1]) - parameters['curyear'] = real2float(fields[2]) - parameters['restart'] = fields[3] - parameters['fracdone'] = real2float(fields[4]) - parameters['masstot'] = real2float(fields[5]) - - #Parse remainder of file to build seed array - nlines = len(lines) - index = 1 - while (index < nlines): - fields = lines[index].split() - seedarr[index - 1] = real2float(fields[0]) - index += 1 - - parameters['seedn'] = index - 1 - - return - -def read_impact_mass(filename): - #Read impact mass file - - fp=open(filename,'r') - line=fp.readlines() - fp.close() - - fields = line[0].split() - if (len(fields) > 0): - mass = real2float(fields[0]) - else: - mass = 0 - - return mass \ No newline at end of file diff --git a/python/ctem/ctem/ctem_io_writers.py b/python/ctem/ctem/ctem_io_writers.py deleted file mode 100644 index 7bf21ce4..00000000 --- a/python/ctem/ctem/ctem_io_writers.py +++ /dev/null @@ -1,264 +0,0 @@ -#!/usr/local/bin/python -# -#Cratered Terrain Evolution Model file and graphical writing utilities -# -#Original IDL design: Jim Richardson, Arecibo Observatory -#Revised IDL design: David Minton, Purdue University -#Re-engineered design and Python implementation: Matthew Route, Purdue University -#August 2016 -# - -import matplotlib -from matplotlib import pyplot -import numpy -import os -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 - -#Write production function to file production.dat -#This file format does not exactly match that generated from IDL. Does it work? -def write_production(parameters, production): - filename = parameters['workingdir'] + parameters['sfdfile'] - numpy.savetxt(filename, production, fmt='%1.8e', delimiter=' ') - - return - -def create_dir_structure(parameters): - #Create directories for various output files if they do not already exist - directories=['dist','misc','rego','rplot','surf','shaded'] - - for directory in directories: - dir_test = parameters['workingdir'] + directory - if not os.path.isdir(dir_test): - os.makedirs(dir_test) - - return - -def write_ctemdat(parameters, seedarr): - #Write various parameters and random number seeds into ctem.dat file - filename = parameters['workingdir'] + parameters['datfile'] - fp = open(filename,'w') - - template = "%(totalimpacts)17d %(ncount)12d %(curyear)19.12E %(restart)s %(fracdone)9.6f %(masstot)19.12E\n" - fp.write(template % parameters) - - #Write random number seeds to the file - for index in range(parameters['seedn']): - fp.write("%12d\n" % seedarr[index]) - - fp.close() - - return - -def copy_dists(parameters): - #Save copies of distribution files - - orig_list = ['odistribution', 'ocumulative', 'pdistribution', 'tdistribution'] - dest_list = ['odist', 'ocum', 'pdist', 'tdist'] - - for index in range(len(orig_list)): - forig = parameters['workingdir'] + orig_list[index] + '.dat' - fdest = parameters['workingdir'] + 'dist' + os.sep + dest_list[index] + "_%06d.dat" % parameters['ncount'] - shutil.copy2(forig, fdest) - - forig = parameters['workingdir'] + 'impactmass.dat' - fdest = parameters['workingdir'] + 'misc' + os.sep + "mass_%06d.dat" % parameters['ncount'] - shutil.copy2(forig, fdest) - - if (parameters['savetruelist'].upper() == 'T'): - forig = parameters['workingdir'] + 'tcumulative.dat' - fdest = parameters['workingdir'] + 'dist' + os.sep + "tcum_%06d.dat" % parameters['ncount'] - shutil.copy2(forig, fdest) - - return - -#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, 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'] - matplotlib.pyplot.savefig(filename,dpi=dpi,bbox_inches=0) - - return - -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'] - - ls = LightSource(azdeg=azimuth, altdeg=solar_angle) - cmap = cm.cividis - - # If min and max appear to be reversed, then fix them - if (parameters['shadedminh'] > parameters['shadedmaxh']): - temp = parameters['shadedminh'] - 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'] - if (parameters['shadedmaxhdefault'] == 1): - shadedmaxh = numpy.amax(DEM) - else: - 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(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 - dfrac = 2**(1./4) * 1.0e-3 - - #Calculate geometric saturation - minx = (parameters['pix'] / 3.0) * 1.0e-3 - maxx = 3 * parameters['pix'] * parameters['gridsize'] * 1.0e-3 - geomem = numpy.array([[minx, satlimit / 20.0], [maxx, satlimit / 20.0]]) - geomep = numpy.array([[minx, satlimit / 10.0], [maxx, satlimit / 10.0]]) - - #Create distribution arrays without zeros for plotting on log scale - idx = numpy.nonzero(odist[:,5]) - odistnz = odist[idx] - odistnz = odistnz[:,[2,5]] - - idx = numpy.nonzero(pdist[:,5]) - pdistnz = pdist[idx] - pdistnz = pdistnz[:,[2,5]] - - idx = numpy.nonzero(tdist[:,5]) - tdistnz = tdist[idx] - tdistnz = tdistnz[:,[2,5]] - - #Correct pdist - pdistnz[:,1] = pdistnz[:,1] * parameters['curyear'] / parameters['interval'] - - #Create sdist bin factors, which contain one crater per bin - area = (parameters['gridsize'] * parameters['pix'] * 1.0e-3)**2. - plo = 1 - sq2 = 2**(1./2) - while (sq2**plo > minx): - plo = plo - 1 - phi = plo + 1 - while (sq2**phi < maxx): - phi = phi + 1 - n = phi - plo + 1 - sdist = numpy.zeros([n , 2]) - p = plo - for index in range(n): - sdist[index, 0] = sq2**p - sdist[index, 1] = sq2**(2.0*p + 1.5) / (area * (sq2 - 1)) - p = p + 1 - - #Create time label - tlabel = "%5.4e" % parameters['curyear'] - tlabel = tlabel.split('e') - texp = str(int(tlabel[1])) - timelabel = 'Time = '+ r'${}$ x 10$^{}$'.format(tlabel[0], texp) + ' yrs' - - #Save image to file - filename = parameters['workingdir'] + 'rplot' + os.sep + "rplot%06d.png" % parameters['ncount'] - height = parameters['gridsize'] / dpi - width = height - fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi) - - #Alter background color to be black, and change axis colors accordingly - matplotlib.pyplot.style.use('dark_background') - matplotlib.pyplot.rcParams['axes.prop_cycle'] - - #Plot data - matplotlib.pyplot.plot(odistnz[:,0]*1.0e-3, odistnz[:,1], linewidth=3.0, color = 'blue') - matplotlib.pyplot.plot(pdistnz[:,0]*1.0e-3, pdistnz[:,1], linewidth=2.0, linestyle='dashdot', color = 'white') - matplotlib.pyplot.plot(tdistnz[:,0]*1.0e-3, tdistnz[:,1], linewidth=2.0, color = 'red') - - matplotlib.pyplot.plot(geomem[:,0], geomem[:,1], linewidth=2.0, linestyle =':', color = 'yellow') - matplotlib.pyplot.plot(geomep[:,0], geomep[:,1], linewidth=2.0, linestyle =':', color = 'yellow') - matplotlib.pyplot.plot(sdist[:,0], sdist[:,1], linewidth=2.0, linestyle =':', color = 'yellow') - - matplotlib.pyplot.plot(ph1[:,0] * dfrac, ph1[:,1], 'wo') - - #Create plot labels - matplotlib.pyplot.title('Crater Distribution R-Plot',fontsize=22) - matplotlib.pyplot.xlim([minx, maxx]) - matplotlib.pyplot.xscale('log') - matplotlib.pyplot.xlabel('Crater Diameter (km)',fontsize=18) - matplotlib.pyplot.ylim([5.0e-4, 5.0]) - matplotlib.pyplot.yscale('log') - matplotlib.pyplot.ylabel('R Value', fontsize=18) - matplotlib.pyplot.text(1.0e-2, 1.0, timelabel, fontsize=18) - - matplotlib.pyplot.tick_params(axis='both', which='major', labelsize=14) - matplotlib.pyplot.tick_params(axis='both', which='minor', labelsize=12) - - matplotlib.pyplot.savefig(filename) - - return \ No newline at end of file diff --git a/python/ctem/ctem/dist_reader.py b/python/ctem/ctem/dist_reader.py deleted file mode 100644 index 31793793..00000000 --- a/python/ctem/ctem/dist_reader.py +++ /dev/null @@ -1,29 +0,0 @@ -import numpy as np -import pandas as pd - -gridsize = 2000 -pix = 3.08e3 -side = (gridsize * pix) -area = side**2 -iside = [-side, 0.0, side] -basin_crater_cutoff = 300000.0 #diameter in meters - -craters = pd.read_csv('NPF-global-Kd0.0001/dist/ocum_000100.dat',delim_whitespace=True) - -basins = craters.loc[craters['#Dcrat(m)'] >= basin_crater_cutoff] - -for bind, b in basins.iterrows(): - print(f"Basin{bind:02}: ",b['#Dcrat(m)'],b['time(y)']) - overlaps = open(f"NPF-global-Kd0.0001/basin_overlap/overlap{bind:02}.dat", 'w') - print(f"#BASIN NUMBER {bind}: Diameter = {b['#Dcrat(m)']}, Time {b['time(y)']}", file=overlaps) - for cind, c in craters.iterrows(): - if c['time(y)'] > b['time(y)']: - for i in iside: - for j in iside: - disx = b['xpos(m)'] - c['xpos(m)'] + i - disy = b['ypos(m)'] - c['ypos(m)'] + j - distance = np.sqrt(disx ** 2 + disy ** 2) - if distance < 0.5 * (b['#Dcrat(m)'] + c['#Dcrat(m)']): - print(c['#Dcrat(m)'],c['xpos(m)'],c['ypos(m)'],c['time(y)'], file=overlaps) - - diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py new file mode 100644 index 00000000..764759cc --- /dev/null +++ b/python/ctem/ctem/driver.py @@ -0,0 +1,263 @@ +import numpy as np +import os +import subprocess +import shutil +from ctem import io + +class Simulation: + """ + This is a class that defines the basic CTEM simulation object. It initializes a dictionary of user and reads + it in from a file (default name is ctem.in). It also creates the directory structure needed to store simulation + files if necessary. + """ + def __init__(self, param_file="ctem.in"): + currentdir = os.getcwd() + self.user = { + 'restart': None, + 'runtype': None, + 'popupconsole': None, + 'saveshaded': None, + 'saverego': None, + 'savepres': None, + 'savetruelist': None, + 'seedn': 1, + 'totalimpacts': 0, + 'ncount': 0, + 'curyear': 0.0, + 'fracdone': 1.0, + 'masstot': 0.0, + 'interval': 0.0, + 'numintervals': 0, + 'pix': -1.0, + 'gridsize': -1, + 'seed': 0, + 'maxcrat': 1.0, + 'shadedminhdefault': 1, + 'shadedmaxhdefault': 1, + 'shadedminh': 0.0, + 'shadedmaxh': 0.0, + 'workingdir': currentdir, + 'ctemfile': os.path.join(currentdir, param_file), + 'impfile': None, + 'sfdcompare': None, + 'sfdfile': None + } + self.user = io.read_user_input(self.user) + + self.output_filenames = { + 'dat': 'ctem.dat', + 'dem': 'surface_dem.dat', + 'diam': 'surface_diam.dat', + 'ejc': 'surface_ejc.dat', + 'pos': 'surface_pos.dat', + 'time': 'surface_time.dat', + 'ocum': 'ocumulative.dat', + 'odist': 'odistribution.dat', + 'pdist': 'pdistribution.dat', + 'tcum' : 'tcumulative.dat', + 'tdist' : 'tdistribution.dat', + 'impmass' : 'impactmass.dat', + 'fracdone' : 'fracdone.dat', + 'regodepth' : 'regolithdepth.dat', + } + + for k, v in self.output_filenames.items(): + self.output_filenames[k] = os.path.join(currentdir, v) + + # Set up data arrays + self.seedarr = np.zeros(100, dtype=int) + self.seedarr[0] = self.user['seed'] + self.odist = np.zeros([1, 6]) + self.pdist = np.zeros([1, 6]) + self.tdist = np.zeros([1, 6]) + self.surface_dem = np.zeros([self.user['gridsize'], self.user['gridsize']], dtype=float) + self.regolith = np.zeros([self.user['gridsize'], self.user['gridsize']], dtype=float) + + if self.user['sfdcompare'] is not None: + # Read sfdcompare file + sfdfile = os.path.join(self.user['workingdir'], self.user['sfdcompare']) + self.ph1 = io.read_formatted_ascii(sfdfile, skip_lines=0) + + # Scale the production function to the simulation domain + self.scale_production() + + # Starting new or old run? + if (self.user['restart'].upper() == 'F'): + print('Starting a new run') + + io.create_dir_structure(self.user) + + if (self.user['runtype'].upper() == 'STATISTICAL'): + self.user['ncount'] = 1 + + # Write ctem.dat file + io.write_datfile(self.user, self.seedarr) + else: + self.user['ncount'] = 0 + + # Delete any old output files + for k, v in self.output_filenames.items(): + if os.path.isfile(v): + os.remove(v) + else: + print('Continuing a previous run') + self.process_interval(isnew=False) + + + return + + def scale_production(self): + """ + Scales the production function to the simulation domain and saves the scaled production function to a file that + will be read in by the main Fortran program. + """ + + # Read production function file + impfile = os.path.join(self.user['workingdir'], self.user['impfile']) + prodfunction = io.read_formatted_ascii(impfile, skip_lines=0) + + # Create impactor production population + area = (self.user['gridsize'] * self.user['pix']) ** 2 + self.production = np.copy(prodfunction) + self.production[:, 1] = self.production[:, 1] * area * self.user['interval'] + + # Write the scaled production function to file + io.write_production(self.user, self.production) + + return + + def run(self): + """ + Runs a complete simulation over all intervals specified in the input user. This method loops over all + intervals and process outputs into images, and redirect output files to their apporpriate folders. + + This method replaces most of the functionality of the original ctem_driver.py + """ + + print('Beginning loops') + while (self.user['ncount'] <= self.user['numintervals']): + if (self.user['ncount'] > 0): + # Move ctem.dat + self.redirect_outputs(['dat'], 'misc') + + #Execute Fortran code + self.compute_one_interval() + + # Process output files + self.process_interval() + + # Display results + print(self.user['ncount'], ' Displaying results') + + # Write surface dem, surface ejecta, shaded relief, and rplot data + io.image_dem(self.user, self.surface_dem) + if (self.user['saverego'].upper() == 'T'): + io.image_regolith(self.user, self.surface_ejc) + if (self.user['saveshaded'].upper() == 'T'): + io.image_shaded_relief(self.user, self.surface_dem) + if (self.user['savepres'].upper() == 'T'): + io.create_rplot(self.user, self.odist, self.pdist, self.tdist, self.ph1) + + # Update ncount + self.user['ncount'] = self.user['ncount'] + 1 + + if ((self.user['runtype'].upper() == 'STATISTICAL') or (self.user['ncount'] == 1)): + self.user['restart'] = 'F' + self.user['curyear'] = 0.0 + self.user['totalimpacts'] = 0 + self.user['masstot'] = 0.0 + + # Delete tdistribution file, if it exists + tdist_file = self.user['workingdir'] + 'tdistribution.dat' + if os.path.isfile(tdist_file): + os.remove(tdist_file) + + else: + self.user['restart'] = 'T' + + # Write ctem.dat file + io.write_datfile(self.user, self.output_filenames['dat'],self.seedarr) + + return + + def compute_one_interval(self): + """ + Executes the Fortran code to generate one interval of simulation output. + """ + # Create crater population and display CTEM progress on screen + print(self.user['ncount'], ' Calling FORTRAN routine') + with subprocess.Popen([os.path.join(self.user['workingdir'], 'CTEM')], + stdout=subprocess.PIPE, + universal_newlines=True) as p: + for line in p.stdout: + print(line, end='') + + return + + def process_interval(self, isnew=True): + """ + Reads in all of the files generated by the Fortran code after one interval and redirects them to appropriate + folders for storage after appending the interval number to the filename. + """ + # Read Fortran output + print(self.user['ncount'], ' Reading Fortran output') + + # Read surface dem(shaded relief) and ejecta data files + self.surface_dem = io.read_unformatted_binary(self.output_filenames['dem'], self.user['gridsize']) + self.surface_ejc = io.read_unformatted_binary(self.output_filenames['ejc'], self.user['gridsize']) + + # Read odistribution, tdistribution, and pdistribution files + self.odist = io.read_formatted_ascii(self.output_filenames['odist'], skip_lines=1) + self.tdist = io.read_formatted_ascii(self.output_filenames['tdist'], skip_lines=1) + self.pdist = io.read_formatted_ascii(self.output_filenames['pdist'], skip_lines=1) + + # Read impact mass from file + impact_mass = io.read_impact_mass(self.output_filenames['massfile']) + + # Read ctem.dat file + io.read_ctemdat(self.user, self.seedarr) + + # Save copy of crater distribution files + if isnew: + + # Update user: mass, curyear, regolith properties + self.user['masstot'] = self.user['masstot'] + impact_mass + + self.user['curyear'] = self.user['curyear'] + self.user['fracdone'] * self.user['interval'] + template = "%(fracdone)9.6f %(curyear)19.12E\n" + with open(self.output_filenames['fracdone']) as fp_frac: + fp_frac.write(template % self.user) + + reg_text = "%19.12E %19.12E %19.12E %19.12E\n" % (self.user['curyear'], + np.mean(self.surface_ejc), np.amax(self.surface_ejc), + np.amin(self.surface_ejc)) + with open(self.output_filenames['regodepth']) as fp_reg: + fp_reg.write(reg_text) + + self.redirect_outputs(['odist', 'ocum', 'pdist', 'tdist'], 'dist') + if (self.user['savetruelist'].upper() == 'T'): + self.redirect_outputs(['tcum'], 'dist') + self.redirect_outputs(['impmass'], 'misc') + return + + def redirect_outputs(self, filekeys, foldername): + """ + Copies a set of output files from the working directory into a subfolder. + Takes as input the dictionaries of user parameters and output file names, the dictionary keys to the file names + you wish to redirect, and the redirection destination folder name. The new name will be the key name + zero padded + interval number. + + Example: + calling redirect_outputs(['odist','tdist'], 'dist') when user['ncount'] is 1 + should do the following: + copy 'odistribution.dat' to 'dist/odist_000001.dat' + copy 'tdistribution.dat' to 'dist/tdist_000001.dat' + """ + + for k in filekeys: + forig = self.output_filenames[k] + fdist = os.path.join(self.user['workingdir'], foldername, f"{k}_{self.user['ncount']:06d}.dat") + +if __name__ == '__main__': + sim = Simulation() + sim.run() \ No newline at end of file diff --git a/python/ctem/ctem/frontend.py b/python/ctem/ctem/frontend.py deleted file mode 100644 index f7ee95ea..00000000 --- a/python/ctem/ctem/frontend.py +++ /dev/null @@ -1,125 +0,0 @@ -import numpy as np -import os -import subprocess -import shutil -from ctem import io - -class Simulation: - """ - This is a class that defines the basic CTEM simulation object - """ - def __init__(self, param_file="ctem.in"): - currentdir = os.getcwd() - self.parameters = { - 'restart': None, - 'runtype': None, - 'popupconsole': None, - 'saveshaded': None, - 'saverego': None, - 'savepres': None, - 'savetruelist': None, - 'seedn': 1, - 'totalimpacts': 0, - 'ncount': 0, - 'curyear': 0.0, - 'fracdone': 1.0, - 'masstot': 0.0, - 'interval': 0.0, - 'numintervals': 0, - 'pix': -1.0, - 'gridsize': -1, - 'seed': 0, - 'maxcrat': 1.0, - 'shadedminhdefault': 1, - 'shadedmaxhdefault': 1, - 'shadedminh': 0.0, - 'shadedmaxh': 0.0, - 'workingdir': currentdir, - 'ctemfile': param_file, - 'datfile': 'ctem.dat', - 'impfile': None, - 'sfdcompare': None, - 'sfdfile': None - } - self.parameters = io.read_param(self.parameters) - - # Set up data arrays - self.seedarr = np.zeros(100, dtype=int) - self.seedarr[0] = self.parameters['seed'] - self.odist = np.zeros([1, 6]) - self.pdist = np.zeros([1, 6]) - self.tdist = np.zeros([1, 6]) - self.surface_dem = np.zeros([self.parameters['gridsize'], self.parameters['gridsize']], dtype=np.float) - self.regolith = np.zeros([self.parameters['gridsize'], self.parameters['gridsize']], dtype=np.float) - - if self.parameters['sfdcompare'] is not None: - # Read sfdcompare file - sfdfile = os.path.join(self.parameters['workingdir'], self.parameters['sfdcompare']) - self.ph1 = io.read_formatted_ascii(sfdfile, skip_lines=0) - - # Read production function file - impfile = os.path.join(self.parameters['workingdir'], parameters['impfile']) - prodfunction = io.read_formatted_ascii(impfile, skip_lines=0) - - # Create impactor production population - area = (self.parameters['gridsize'] * self.parameters['pix']) ** 2 - self.production = np.copy(prodfunction) - self.production[:, 1] = self.production[:, 1] * area * self.parameters['interval'] - - # Write corrected production function to file - io.write_production(self.parameters, self.production) - - # Starting new or old run? - if (self.parameters['restart'].upper() == 'F'): - print('Starting a new run') - - if (self.parameters['runtype'].upper() == 'STATISTICAL'): - self.parameters['ncount'] = 1 - - # Write ctem.dat file - io.write_ctemdat(self.parameters, self.seedarr) - - else: - self.parameters['ncount'] = 0 - - # Delete tdistribution file, if it exists - tdist_file = os.path.join(self.parameters['workingdir'], 'tdistribution.dat') - if os.path.isfile(tdist_file): - os.remove(tdist_file) - else: - print('Continuing a previous run') - - # Read surface dem(shaded relief) and ejecta data files - dem_file = os.path.join(self.parameters['workingdir'], 'surface_dem.dat') - self.surface_dem = io.read_unformatted_binary(dem_file, self.parameters['gridsize']) - ejecta_file = os.path.join(self.parameters['workingdir'], 'surface_ejc.dat') - self.regolith = io.read_unformatted_binary(ejecta_file, self.parameters['gridsize']) - - # Read odistribution, tdistribution, and pdistribution files - ofile = os.path.join(self.parameters['workingdir'],'odistribution.dat') - odist = io.read_formatted_ascii(ofile, skip_lines=1) - tfile = os.path.join(self.parameters['workingdir'], 'tdistribution.dat') - tdist = io.read_formatted_ascii(tfile, skip_lines=1) - pfile = os.path.join(self.parameters['workingdir'], 'pdistribution.dat') - pdist = io.read_formatted_ascii(pfile, skip_lines=1) - - # Read impact mass from file - massfile = os.path.join(self.parameters['workingdir'], 'impactmass.dat') - impact_mass = io.read_impact_mass(massfile) - - # Read ctem.dat file - io.read_ctemdat(self.parameters, self.seedarr) - - # Open fracdonefile and regodepthfile for writing - # filename = os.path.join(self.parameters['workingdir'], 'fracdone.dat') - # fp_frac = open(filename, 'w') - # filename = os.path.join(self.parameters['workingdir'], 'regolithdepth.dat') - # fp_reg = open(filename, 'w') - - io.create_dir_structure(self.parameters) - - - return - - - # def \ No newline at end of file diff --git a/python/ctem/ctem/io.py b/python/ctem/ctem/io.py index 60ae499e..0b587234 100644 --- a/python/ctem/ctem/io.py +++ b/python/ctem/ctem/io.py @@ -8,49 +8,50 @@ # Set pixel scaling common for image writing, at 1 pixel/ array element dpi = 72.0 -def copy_dists(parameters): +def copy_dists(user, output_filenames, distlist): # Save copies of distribution files orig_list = ['odistribution', 'ocumulative', 'pdistribution', 'tdistribution'] dest_list = ['odist', 'ocum', 'pdist', 'tdist'] for index in range(len(orig_list)): - forig = parameters['workingdir'] + orig_list[index] + '.dat' - fdest = parameters['workingdir'] + 'dist' + os.sep + dest_list[index] + "_%06d.dat" % parameters['ncount'] + forig = user['workingdir'] + orig_list[index] + '.dat' + fdest = user['workingdir'] + 'dist' + os.sep + dest_list[index] + "_%06d.dat" % user['ncount'] shutil.copy2(forig, fdest) - forig = parameters['workingdir'] + 'impactmass.dat' - fdest = parameters['workingdir'] + 'misc' + os.sep + "mass_%06d.dat" % parameters['ncount'] + forig = user['workingdir'] + 'impactmass.dat' + fdest = user['workingdir'] + 'misc' + os.sep + "mass_%06d.dat" % user['ncount'] shutil.copy2(forig, fdest) - if (parameters['savetruelist'].upper() == 'T'): - forig = parameters['workingdir'] + 'tcumulative.dat' - fdest = parameters['workingdir'] + 'dist' + os.sep + "tcum_%06d.dat" % parameters['ncount'] + if (user['savetruelist'].upper() == 'T'): + forig = user['workingdir'] + 'tcumulative.dat' + fdest = user['workingdir'] + 'dist' + os.sep + "tcum_%06d.dat" % user['ncount'] shutil.copy2(forig, fdest) return -def create_dir_structure(parameters): + +def create_dir_structure(user): # Create directories for various output files if they do not already exist directories = ['dist', 'misc', 'rego', 'rplot', 'surf', 'shaded'] for directory in directories: - dir_test = parameters['workingdir'] + directory + dir_test = user['workingdir'] + directory if not os.path.isdir(dir_test): os.makedirs(dir_test) return -def create_rplot(parameters, odist, pdist, tdist, ph1): +def create_rplot(user, odist, pdist, tdist, ph1): # Parameters: empirical saturation limit and dfrac satlimit = 3.12636 dfrac = 2 ** (1. / 4) * 1.0e-3 # Calculate geometric saturation - minx = (parameters['pix'] / 3.0) * 1.0e-3 - maxx = 3 * parameters['pix'] * parameters['gridsize'] * 1.0e-3 + minx = (user['pix'] / 3.0) * 1.0e-3 + maxx = 3 * user['pix'] * user['gridsize'] * 1.0e-3 geomem = np.array([[minx, satlimit / 20.0], [maxx, satlimit / 20.0]]) geomep = np.array([[minx, satlimit / 10.0], [maxx, satlimit / 10.0]]) @@ -68,10 +69,10 @@ def create_rplot(parameters, odist, pdist, tdist, ph1): tdistnz = tdistnz[:, [2, 5]] # Correct pdist - pdistnz[:, 1] = pdistnz[:, 1] * parameters['curyear'] / parameters['interval'] + pdistnz[:, 1] = pdistnz[:, 1] * user['curyear'] / user['interval'] # Create sdist bin factors, which contain one crater per bin - area = (parameters['gridsize'] * parameters['pix'] * 1.0e-3) ** 2. + area = (user['gridsize'] * user['pix'] * 1.0e-3) ** 2. plo = 1 sq2 = 2 ** (1. / 2) while (sq2 ** plo > minx): @@ -88,14 +89,14 @@ def create_rplot(parameters, odist, pdist, tdist, ph1): p = p + 1 # Create time label - tlabel = "%5.4e" % parameters['curyear'] + tlabel = "%5.4e" % user['curyear'] tlabel = tlabel.split('e') texp = str(int(tlabel[1])) timelabel = 'Time = ' + r'${}$ x 10$^{}$'.format(tlabel[0], texp) + ' yrs' # Save image to file - filename = parameters['workingdir'] + 'rplot' + os.sep + "rplot%06d.png" % parameters['ncount'] - height = parameters['gridsize'] / dpi + filename = user['workingdir'] + 'rplot' + os.sep + "rplot%06d.png" % user['ncount'] + height = user['gridsize'] / dpi width = height fig = plt.figure(figsize=(width, height), dpi=dpi) @@ -131,13 +132,13 @@ def create_rplot(parameters, odist, pdist, tdist, ph1): return -def image_dem(parameters, DEM): +def image_dem(user, DEM): dpi = 300.0 # 72.0 - pix = parameters['pix'] - gridsize = parameters['gridsize'] + pix = user['pix'] + gridsize = user['gridsize'] ve = 1.0 - azimuth = 300.0 # parameters['azimuth'] - solar_angle = 20.0 # parameters['solar_angle'] + azimuth = 300.0 # user['azimuth'] + 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) @@ -150,15 +151,15 @@ def image_dem(parameters, DEM): ax.imshow(dem_img, interpolation="nearest", cmap='gray', vmin=0.0, vmax=1.0) plt.axis('off') # Save image to file - filename = parameters['workingdir'] + 'surf' + os.sep + "surf%06d.png" % parameters['ncount'] + filename = user['workingdir'] + 'surf' + os.sep + "surf%06d.png" % user['ncount'] plt.savefig(filename, dpi=dpi, bbox_inches=0) return -def image_regolith(parameters, regolith): +def image_regolith(user, regolith): # Create scaled regolith image - minref = parameters['pix'] * 1.0e-4 + minref = user['pix'] * 1.0e-4 maxreg = np.amax(regolith) minreg = np.amin(regolith) if (minreg < minref): minreg = minref @@ -169,8 +170,8 @@ def image_regolith(parameters, regolith): (np.log(regolith_scaled) - np.log(minreg)) / (np.log(maxreg) - np.log(minreg))) # Save image to file - filename = parameters['workingdir'] + 'rego' + os.sep + "rego%06d.png" % parameters['ncount'] - height = parameters['gridsize'] / dpi + filename = user['workingdir'] + 'rego' + os.sep + "rego%06d.png" % user['ncount'] + height = user['gridsize'] / dpi width = height fig = plt.figure(figsize=(width, height), dpi=dpi) fig.figimage(regolith_scaled, cmap=cm.nipy_spectral, origin='lower') @@ -179,36 +180,36 @@ def image_regolith(parameters, regolith): return -def image_shaded_relief(parameters, DEM): +def image_shaded_relief(user, DEM): dpi = 300.0 # 72.0 - pix = parameters['pix'] - gridsize = parameters['gridsize'] + pix = user['pix'] + gridsize = user['gridsize'] ve = 1.0 mode = 'overlay' - azimuth = 300.0 # parameters['azimuth'] - solar_angle = 20.0 # parameters['solar_angle'] + azimuth = 300.0 # user['azimuth'] + solar_angle = 20.0 # user['solar_angle'] ls = LightSource(azdeg=azimuth, altdeg=solar_angle) cmap = cm.cividis # If min and max appear to be reversed, then fix them - if (parameters['shadedminh'] > parameters['shadedmaxh']): - temp = parameters['shadedminh'] - parameters['shadedminh'] = parameters['shadedmaxh'] - parameters['shadedmaxh'] = temp + if (user['shadedminh'] > user['shadedmaxh']): + temp = user['shadedminh'] + user['shadedminh'] = user['shadedmaxh'] + user['shadedmaxh'] = temp else: - parameters['shadedminh'] = parameters['shadedminh'] - parameters['shadedmaxh'] = parameters['shadedmaxh'] + user['shadedminh'] = user['shadedminh'] + user['shadedmaxh'] = user['shadedmaxh'] - # If no shadedmin/max parameters are read in from ctem.dat, determine the values from the data - if (parameters['shadedminhdefault'] == 1): + # If no shadedmin/max user are read in from ctem.dat, determine the values from the data + if (user['shadedminhdefault'] == 1): shadedminh = np.amin(DEM) else: - shadedminh = parameters['shadedminh'] - if (parameters['shadedmaxhdefault'] == 1): + shadedminh = user['shadedminh'] + if (user['shadedmaxhdefault'] == 1): shadedmaxh = np.amax(DEM) else: - shadedmaxh = parameters['shadedmaxh'] + shadedmaxh = user['shadedmaxh'] dem_img = ls.shade(DEM, cmap=cmap, blend_mode=mode, fraction=1.0, vert_exag=ve, dx=pix, dy=pix, @@ -222,17 +223,17 @@ def image_shaded_relief(parameters, DEM): ax.imshow(dem_img, interpolation="nearest", vmin=0.0, vmax=1.0) plt.axis('off') # Save image to file - filename = parameters['workingdir'] + 'shaded' + os.sep + "shaded%06d.png" % parameters['ncount'] + filename = user['workingdir'] + 'shaded' + os.sep + "shaded%06d.png" % user['ncount'] plt.savefig(filename, dpi=dpi, bbox_inches=0) - return parameters + return user -def read_ctemdat(parameters, seedarr): +def read_datfile(user, seedarr): # Read and parse ctem.dat file - datfile = parameters['workingdir'] + 'ctem.dat' + datfile = user['workingdir'] + 'ctem.dat' # Read ctem.dat file - print('Reading input file ' + parameters['datfile']) + print('Reading input file ' + user['datfile']) fp = open(datfile, 'r') lines = fp.readlines() fp.close() @@ -240,12 +241,12 @@ def read_ctemdat(parameters, seedarr): # Parse file lines and update parameter fields fields = lines[0].split() if len(fields) > 0: - parameters['totalimpacts'] = real2float(fields[0]) - parameters['ncount'] = int(fields[1]) - parameters['curyear'] = real2float(fields[2]) - parameters['restart'] = fields[3] - parameters['fracdone'] = real2float(fields[4]) - parameters['masstot'] = real2float(fields[5]) + user['totalimpacts'] = real2float(fields[0]) + user['ncount'] = int(fields[1]) + user['curyear'] = real2float(fields[2]) + user['restart'] = fields[3] + user['fracdone'] = real2float(fields[4]) + user['masstot'] = real2float(fields[5]) # Parse remainder of file to build seed array nlines = len(lines) @@ -255,7 +256,7 @@ def read_ctemdat(parameters, seedarr): seedarr[index - 1] = real2float(fields[0]) index += 1 - parameters['seedn'] = index - 1 + user['seedn'] = index - 1 return @@ -285,12 +286,12 @@ def read_impact_mass(filename): # Write production function to file production.dat # This file format does not exactly match that generated from IDL. Does it work? -def read_param(parameters): +def read_user_input(user): # Read and parse ctem.in file - inputfile = os.path.join(parameters['workingdir'], parameters['ctemfile']) + inputfile = user['ctemfile'] # Read ctem.in file - print('Reading input file ' + parameters['ctemfile']) + print('Reading input file ' + user['ctemfile']) with open(inputfile, 'r') as fp: lines = fp.readlines() @@ -298,60 +299,60 @@ def read_param(parameters): for line in lines: fields = line.split() if len(fields) > 0: - if ('pix' == fields[0].lower()): parameters['pix'] = real2float(fields[1]) - if ('gridsize' == fields[0].lower()): parameters['gridsize'] = int(fields[1]) - if ('seed' == fields[0].lower()): parameters['seed'] = int(fields[1]) - if ('sfdfile' == fields[0].lower()): parameters['sfdfile'] = fields[1] - if ('impfile' == fields[0].lower()): parameters['impfile'] = fields[1] - if ('maxcrat' == fields[0].lower()): parameters['maxcrat'] = real2float(fields[1]) - if ('sfdcompare' == fields[0].lower()): parameters['sfdcompare'] = fields[1] - if ('interval' == fields[0].lower()): parameters['interval'] = real2float(fields[1]) - if ('numintervals' == fields[0].lower()): parameters['numintervals'] = int(fields[1]) - if ('popupconsole' == fields[0].lower()): parameters['popupconsole'] = fields[1] - if ('saveshaded' == fields[0].lower()): parameters['saveshaded'] = fields[1] - if ('saverego' == fields[0].lower()): parameters['saverego'] = fields[1] - if ('savepres' == fields[0].lower()): parameters['savepres'] = fields[1] - if ('savetruelist' == fields[0].lower()): parameters['savetruelist'] = fields[1] - if ('runtype' == fields[0].lower()): parameters['runtype'] = fields[1] - if ('restart' == fields[0].lower()): parameters['restart'] = fields[1] + if ('pix' == fields[0].lower()): user['pix'] = real2float(fields[1]) + if ('gridsize' == fields[0].lower()): user['gridsize'] = int(fields[1]) + if ('seed' == fields[0].lower()): user['seed'] = int(fields[1]) + if ('sfdfile' == fields[0].lower()): user['sfdfile'] = os.path.join(user['workingdir'],fields[1]) + if ('impfile' == fields[0].lower()): user['impfile'] = os.path.join(user['workingdir'],fields[1]) + if ('maxcrat' == fields[0].lower()): user['maxcrat'] = real2float(fields[1]) + if ('sfdcompare' == fields[0].lower()): user['sfdcompare'] = os.path.join(user['workingdir'], fields[1]) + if ('interval' == fields[0].lower()): user['interval'] = real2float(fields[1]) + if ('numintervals' == fields[0].lower()): user['numintervals'] = int(fields[1]) + if ('popupconsole' == fields[0].lower()): user['popupconsole'] = fields[1] + if ('saveshaded' == fields[0].lower()): user['saveshaded'] = fields[1] + if ('saverego' == fields[0].lower()): user['saverego'] = fields[1] + if ('savepres' == fields[0].lower()): user['savepres'] = fields[1] + if ('savetruelist' == fields[0].lower()): user['savetruelist'] = fields[1] + if ('runtype' == fields[0].lower()): user['runtype'] = fields[1] + if ('restart' == fields[0].lower()): user['restart'] = fields[1] if ('shadedminh' == fields[0].lower()): - parameters['shadedminh'] = real2float(fields[1]) - parameters['shadedminhdefault'] = 0 + user['shadedminh'] = real2float(fields[1]) + user['shadedminhdefault'] = 0 if ('shadedmaxh' == fields[0].lower()): - parameters['shadedmaxh'] = real2float(fields[1]) - parameters['shadedmaxhdefault'] = 0 + user['shadedmaxh'] = real2float(fields[1]) + user['shadedmaxhdefault'] = 0 # Test values for further processing - if (parameters['interval'] <= 0.0): + if (user['interval'] <= 0.0): print('Invalid value for or missing variable INTERVAL in ' + inputfile) - if (parameters['numintervals'] <= 0): + if (user['numintervals'] <= 0): print('Invalid value for or missing variable NUMINTERVALS in ' + inputfile) - if (parameters['pix'] <= 0.0): + if (user['pix'] <= 0.0): print('Invalid value for or missing variable PIX in ' + inputfile) - if (parameters['gridsize'] <= 0): + if (user['gridsize'] <= 0): print('Invalid value for or missing variable GRIDSIZE in ' + inputfile) - if (parameters['seed'] == 0): + if (user['seed'] == 0): print('Invalid value for or missing variable SEED in ' + inputfile) - if (parameters['sfdfile'] is None): + if (user['sfdfile'] is None): print('Invalid value for or missing variable SFDFILE in ' + inputfile) - if (parameters['impfile'] is None): + if (user['impfile'] is None): print('Invalid value for or missing variable IMPFILE in ' + inputfile) - if (parameters['popupconsole'] is None): + if (user['popupconsole'] is None): print('Invalid value for or missing variable POPUPCONSOLE in ' + inputfile) - if (parameters['saveshaded'] is None): + if (user['saveshaded'] is None): print('Invalid value for or missing variable SAVESHADED in ' + inputfile) - if (parameters['saverego'] is None): + if (user['saverego'] is None): print('Invalid value for or missing variable SAVEREGO in ' + inputfile) - if (parameters['savepres'] is None): + if (user['savepres'] is None): print('Invalid value for or missing variable SAVEPRES in ' + inputfile) - if (parameters['savetruelist'] is None): + if (user['savetruelist'] is None): print('Invalid value for or missing variable SAVETRUELIST in ' + inputfile) - if (parameters['runtype'] is None): + if (user['runtype'] is None): print('Invalid value for or missing variable RUNTYPE in ' + inputfile) - if (parameters['restart'] is None): + if (user['restart'] is None): print('Invalid value for or missing variable RESTART in ' + inputfile) - return parameters + return user def read_unformatted_binary(filename, gridsize): @@ -382,16 +383,15 @@ def real2float(realstr): return float(realstr.replace('d', 'E').replace('D', 'E')) -def write_ctemdat(parameters, seedarr): - # Write various parameters and random number seeds into ctem.dat file - filename = parameters['workingdir'] + parameters['datfile'] +def write_datfile(user, filename, seedarr): + # Write various user and random number seeds into ctem.dat file fp = open(filename, 'w') template = "%(totalimpacts)17d %(ncount)12d %(curyear)19.12E %(restart)s %(fracdone)9.6f %(masstot)19.12E\n" - fp.write(template % parameters) + fp.write(template % user) # Write random number seeds to the file - for index in range(parameters['seedn']): + for index in range(user['seedn']): fp.write("%12d\n" % seedarr[index]) fp.close() @@ -399,11 +399,8 @@ def write_ctemdat(parameters, seedarr): return -# Possible references -# http://nbviewer.jupyter.org/github/ThomasLecocq/geophysique.be/blob/master/2014-02-25%20Shaded%20Relief%20Map%20in%20Python.ipynb - -def write_production(parameters, production): - filename = parameters['workingdir'] + parameters['sfdfile'] +def write_production(user, production): + filename = user['sfdfile'] np.savetxt(filename, production, fmt='%1.8e', delimiter=' ') return diff --git a/python/ctem/ctem/ctem_viewer_3d.py b/python/ctem/ctem/viewer3d.py similarity index 98% rename from python/ctem/ctem/ctem_viewer_3d.py rename to python/ctem/ctem/viewer3d.py index b4517054..684aed85 100644 --- a/python/ctem/ctem/ctem_viewer_3d.py +++ b/python/ctem/ctem/viewer3d.py @@ -14,7 +14,7 @@ def __init__(self): return def readctem(self, surface_file): - # Create and initialize data dictionaries for parameters and options from CTEM.in + # Create and initialize data dictionaries for user and options from CTEM.in notset = '-NOTSET-' currentdir = os.getcwd() + os.sep diff --git a/python/ctem/tests/testio/testio.py b/python/ctem/tests/testio/testio.py index ae329313..1212193e 100644 --- a/python/ctem/tests/testio/testio.py +++ b/python/ctem/tests/testio/testio.py @@ -3,5 +3,5 @@ d = dir(ctem) sim = ctem.Simulation() -for k, v in sim.parameters.items(): +for k, v in sim.user.items(): print(k, v) \ No newline at end of file From c5c6c74d33753dda1629628cc4f5564898ac047d Mon Sep 17 00:00:00 2001 From: David Minton Date: Tue, 22 Feb 2022 12:01:21 -0500 Subject: [PATCH 13/19] Fixed a bunch of typos --- examples/global-lunar-bombardment/ctem.in | 4 ++-- python/ctem/ctem/driver.py | 6 +++--- python/ctem/ctem/io.py | 5 ++--- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/examples/global-lunar-bombardment/ctem.in b/examples/global-lunar-bombardment/ctem.in index 470aaee4..ca50a7fa 100755 --- a/examples/global-lunar-bombardment/ctem.in +++ b/examples/global-lunar-bombardment/ctem.in @@ -34,9 +34,9 @@ runtype single ! Run type: options are normal / ! CTEM required inputs seed 76535 ! Random number generator seed -gridsize 2000 ! Size of grid in pixels +gridsize 200 ! Size of grid in pixels numlayers 10 ! Number of perched layers -pix 3.08e3 ! Pixel size (m) +pix 3.08e4 ! Pixel size (m) mat rock ! Material (rock or ice) ! Bedrock scaling parameters mu_b 0.55e0 ! Experimentally derived parameter for bedrock crater scaling law diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py index 764759cc..061e8098 100644 --- a/python/ctem/ctem/driver.py +++ b/python/ctem/ctem/driver.py @@ -91,7 +91,7 @@ def __init__(self, param_file="ctem.in"): self.user['ncount'] = 1 # Write ctem.dat file - io.write_datfile(self.user, self.seedarr) + io.write_datfile(self.user, self.output_filenames['dat'], self.seedarr) else: self.user['ncount'] = 0 @@ -212,10 +212,10 @@ def process_interval(self, isnew=True): self.pdist = io.read_formatted_ascii(self.output_filenames['pdist'], skip_lines=1) # Read impact mass from file - impact_mass = io.read_impact_mass(self.output_filenames['massfile']) + impact_mass = io.read_impact_mass(self.output_filenames['impmass']) # Read ctem.dat file - io.read_ctemdat(self.user, self.seedarr) + io.read_datfile(self.user, self.output_filenames['dat'], self.seedarr) # Save copy of crater distribution files if isnew: diff --git a/python/ctem/ctem/io.py b/python/ctem/ctem/io.py index 0b587234..92f5dc7a 100644 --- a/python/ctem/ctem/io.py +++ b/python/ctem/ctem/io.py @@ -228,12 +228,11 @@ def image_shaded_relief(user, DEM): return user -def read_datfile(user, seedarr): +def read_datfile(user, datfile, seedarr): # Read and parse ctem.dat file - datfile = user['workingdir'] + 'ctem.dat' # Read ctem.dat file - print('Reading input file ' + user['datfile']) + print('Reading input file ' + datfile) fp = open(datfile, 'r') lines = fp.readlines() fp.close() From 7557f79ea1e98a276cac46a4dc32705826ee97a1 Mon Sep 17 00:00:00 2001 From: David Minton Date: Tue, 22 Feb 2022 12:44:51 -0500 Subject: [PATCH 14/19] Fixed more typos. --- python/ctem/ctem/driver.py | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py index 061e8098..c09cf20a 100644 --- a/python/ctem/ctem/driver.py +++ b/python/ctem/ctem/driver.py @@ -218,22 +218,21 @@ def process_interval(self, isnew=True): io.read_datfile(self.user, self.output_filenames['dat'], self.seedarr) # Save copy of crater distribution files + # Update user: mass, curyear, regolith properties + self.user['masstot'] = self.user['masstot'] + impact_mass + + self.user['curyear'] = self.user['curyear'] + self.user['fracdone'] * self.user['interval'] + template = "%(fracdone)9.6f %(curyear)19.12E\n" + with open(self.output_filenames['fracdone'], 'w') as fp_frac: + fp_frac.write(template % self.user) + + reg_text = "%19.12E %19.12E %19.12E %19.12E\n" % (self.user['curyear'], + np.mean(self.surface_ejc), np.amax(self.surface_ejc), + np.amin(self.surface_ejc)) + with open(self.output_filenames['regodepth'], 'w') as fp_reg: + fp_reg.write(reg_text) + if isnew: - - # Update user: mass, curyear, regolith properties - self.user['masstot'] = self.user['masstot'] + impact_mass - - self.user['curyear'] = self.user['curyear'] + self.user['fracdone'] * self.user['interval'] - template = "%(fracdone)9.6f %(curyear)19.12E\n" - with open(self.output_filenames['fracdone']) as fp_frac: - fp_frac.write(template % self.user) - - reg_text = "%19.12E %19.12E %19.12E %19.12E\n" % (self.user['curyear'], - np.mean(self.surface_ejc), np.amax(self.surface_ejc), - np.amin(self.surface_ejc)) - with open(self.output_filenames['regodepth']) as fp_reg: - fp_reg.write(reg_text) - self.redirect_outputs(['odist', 'ocum', 'pdist', 'tdist'], 'dist') if (self.user['savetruelist'].upper() == 'T'): self.redirect_outputs(['tcum'], 'dist') From ffb03159b02854c497bfcb61c7b83a4fe655c91a Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 23 Feb 2022 12:58:58 -0500 Subject: [PATCH 15/19] Fixed some problems related to directory structure generation and accidentally deleting ctem.dat after generating it --- python/ctem/ctem/driver.py | 19 +++++++++---------- python/ctem/ctem/io.py | 2 +- 2 files changed, 10 insertions(+), 11 deletions(-) diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py index c09cf20a..9bda0aa8 100644 --- a/python/ctem/ctem/driver.py +++ b/python/ctem/ctem/driver.py @@ -10,7 +10,7 @@ class Simulation: it in from a file (default name is ctem.in). It also creates the directory structure needed to store simulation files if necessary. """ - def __init__(self, param_file="ctem.in"): + def __init__(self, param_file="ctem.in", isnew=True): currentdir = os.getcwd() self.user = { 'restart': None, @@ -78,14 +78,19 @@ def __init__(self, param_file="ctem.in"): sfdfile = os.path.join(self.user['workingdir'], self.user['sfdcompare']) self.ph1 = io.read_formatted_ascii(sfdfile, skip_lines=0) - # Scale the production function to the simulation domain - self.scale_production() # Starting new or old run? - if (self.user['restart'].upper() == 'F'): + if (self.user['restart'].upper() == 'F' and isnew): print('Starting a new run') io.create_dir_structure(self.user) + # Delete any old output files + for k, v in self.output_filenames.items(): + if os.path.isfile(v): + os.remove(v) + + # Scale the production function to the simulation domain + self.scale_production() if (self.user['runtype'].upper() == 'STATISTICAL'): self.user['ncount'] = 1 @@ -94,16 +99,10 @@ def __init__(self, param_file="ctem.in"): io.write_datfile(self.user, self.output_filenames['dat'], self.seedarr) else: self.user['ncount'] = 0 - - # Delete any old output files - for k, v in self.output_filenames.items(): - if os.path.isfile(v): - os.remove(v) else: print('Continuing a previous run') self.process_interval(isnew=False) - return def scale_production(self): diff --git a/python/ctem/ctem/io.py b/python/ctem/ctem/io.py index 92f5dc7a..4088c7ea 100644 --- a/python/ctem/ctem/io.py +++ b/python/ctem/ctem/io.py @@ -37,7 +37,7 @@ def create_dir_structure(user): directories = ['dist', 'misc', 'rego', 'rplot', 'surf', 'shaded'] for directory in directories: - dir_test = user['workingdir'] + directory + dir_test = os.path.join(user['workingdir'], directory) if not os.path.isdir(dir_test): os.makedirs(dir_test) From a3d78a9bb6becda9fba45ce255174362aae83140 Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 23 Feb 2022 12:59:19 -0500 Subject: [PATCH 16/19] Restructured the 3d viewer so that it uses a derived class of the main Simulation class --- python/ctem/ctem/viewer3d.py | 75 ++++++++---------------------------- 1 file changed, 15 insertions(+), 60 deletions(-) diff --git a/python/ctem/ctem/viewer3d.py b/python/ctem/ctem/viewer3d.py index 684aed85..7c68d995 100644 --- a/python/ctem/ctem/viewer3d.py +++ b/python/ctem/ctem/viewer3d.py @@ -1,64 +1,22 @@ import numpy as np import os import open3d -from ctem import ctem_io_readers +import ctem import os - -class polysurface: +class Polysurface(ctem.Simulation): """A model of a self-gravitating small body""" def __init__(self): - self.polyfile = "" + ctem.Simulation.__init__(self, isnew=False) # Initialize the data structures, but doesn't start a new run #used for Open3d module self.mesh = open3d.geometry.TriangleMesh() - return - - def readctem(self, surface_file): - # Create and initialize data dictionaries for user and options from CTEM.in - notset = '-NOTSET-' - currentdir = os.getcwd() + os.sep - - self.parameters = {'restart': notset, - 'runtype': notset, - 'popupconsole': notset, - 'saveshaded': notset, - 'saverego': notset, - 'savepres': notset, - 'savetruelist': notset, - 'seedn': 1, - 'totalimpacts': 0, - 'ncount': 0, - 'curyear': 0.0, - 'fracdone': 1.0, - 'masstot': 0.0, - 'interval': 0.0, - 'numintervals': 0, - 'pix': -1.0, - 'gridsize': -1, - 'seed': 0, - 'maxcrat': 1.0, - 'shadedminhdefault': 1, - 'shadedmaxhdefault': 1, - 'shadedminh': 0.0, - 'shadedmaxh': 0.0, - 'workingdir': currentdir, - 'ctemfile': 'ctem.in', - 'datfile': 'ctem.dat', - 'impfile': notset, - 'sfdcompare': notset, - 'sfdfile': notset} - - # Read ctem.in to initialize parameter values based on user input - ctem_io_readers.read_ctemin(self.parameters, notset) - # Read surface dem(shaded relief) and ejecta data files - dem_file = self.parameters['workingdir'] + surface_file - self.dem = ctem_io_readers.read_unformatted_binary(dem_file, self.parameters['gridsize']) + self.process_interval(isnew=False) return def ctem2blockmesh(self): # Build mesh grid - s = self.parameters['gridsize'] - pix = self.parameters['pix'] + s = self.user['gridsize'] + pix = self.user['pix'] ygrid, xgrid = np.mgrid[-s/2:s/2, -s/2:s/2] * pix y0, x0 = ygrid - pix/2, xgrid - pix/2 y1, x1 = ygrid - pix/2, xgrid + pix/2 @@ -124,13 +82,13 @@ def ctem2blockmesh(self): def ctem2trimesh(self): # Build mesh grid - s = self.parameters['gridsize'] - pix = self.parameters['pix'] + s = self.user['gridsize'] + pix = self.user['pix'] ygrid, xgrid = np.mgrid[-s/2:s/2, -s/2:s/2] * pix xvals = xgrid.flatten() yvals = ygrid.flatten() - zvals = self.dem.flatten() + zvals = self.surface_dem.flatten() verts = np.array((xvals, yvals, zvals)).T faces = np.full([2 * (s-1)**2, 3], -1, dtype=np.int64) for j in range(s - 1): @@ -154,9 +112,9 @@ def visualize(self): opt = vis.get_render_option() opt.background_color = np.asarray([0, 0, 0]) - # zmax = np.amax(surf.dem) - # zmin = np.amin(surf.dem) - # cnorm = (surf.dem.flatten() - zmin) / (zmax - zmin) + # zmax = np.amax(self.surface_dem) + # zmin = np.amin(self.surface_dem) + # cnorm = (self.surface_dem.flatten() - zmin) / (zmax - zmin) # cval = plt.cm.terrain(cnorm)[:, :3] self.mesh.paint_uniform_color([0.5, 0.5, 0.5]) @@ -167,12 +125,9 @@ def visualize(self): if __name__ == '__main__': import matplotlib.pyplot as plt - surf = polysurface() - surf.readctem(surface_file="surface_dem.dat") - surf.ctem2trimesh() - surf.visualize() - - + sim = Polysurface() + sim.ctem2trimesh() + sim.visualize() From 769ed1709c52ee534ff2bdf42bcd1f2996270cda Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 23 Feb 2022 13:04:24 -0500 Subject: [PATCH 17/19] Fixed directory structure naming using os.path.join instead of appending separation characters manually --- python/ctem/ctem/io.py | 32 ++++---------------------------- 1 file changed, 4 insertions(+), 28 deletions(-) diff --git a/python/ctem/ctem/io.py b/python/ctem/ctem/io.py index 4088c7ea..d0d383e5 100644 --- a/python/ctem/ctem/io.py +++ b/python/ctem/ctem/io.py @@ -8,30 +8,6 @@ # Set pixel scaling common for image writing, at 1 pixel/ array element dpi = 72.0 -def copy_dists(user, output_filenames, distlist): - # Save copies of distribution files - - orig_list = ['odistribution', 'ocumulative', 'pdistribution', 'tdistribution'] - dest_list = ['odist', 'ocum', 'pdist', 'tdist'] - - for index in range(len(orig_list)): - forig = user['workingdir'] + orig_list[index] + '.dat' - fdest = user['workingdir'] + 'dist' + os.sep + dest_list[index] + "_%06d.dat" % user['ncount'] - shutil.copy2(forig, fdest) - - forig = user['workingdir'] + 'impactmass.dat' - fdest = user['workingdir'] + 'misc' + os.sep + "mass_%06d.dat" % user['ncount'] - shutil.copy2(forig, fdest) - - if (user['savetruelist'].upper() == 'T'): - forig = user['workingdir'] + 'tcumulative.dat' - fdest = user['workingdir'] + 'dist' + os.sep + "tcum_%06d.dat" % user['ncount'] - shutil.copy2(forig, fdest) - - return - - - def create_dir_structure(user): # Create directories for various output files if they do not already exist directories = ['dist', 'misc', 'rego', 'rplot', 'surf', 'shaded'] @@ -95,7 +71,7 @@ def create_rplot(user, odist, pdist, tdist, ph1): timelabel = 'Time = ' + r'${}$ x 10$^{}$'.format(tlabel[0], texp) + ' yrs' # Save image to file - filename = user['workingdir'] + 'rplot' + os.sep + "rplot%06d.png" % user['ncount'] + filename = os.path.join(user['workingdir'],'rplot',"rplot%06d.png" % user['ncount']) height = user['gridsize'] / dpi width = height fig = plt.figure(figsize=(width, height), dpi=dpi) @@ -151,7 +127,7 @@ def image_dem(user, DEM): ax.imshow(dem_img, interpolation="nearest", cmap='gray', vmin=0.0, vmax=1.0) plt.axis('off') # Save image to file - filename = user['workingdir'] + 'surf' + os.sep + "surf%06d.png" % user['ncount'] + filename = os.path.join(user['workingdir'],'surf',"surf%06d.png" % user['ncount']) plt.savefig(filename, dpi=dpi, bbox_inches=0) return @@ -170,7 +146,7 @@ def image_regolith(user, regolith): (np.log(regolith_scaled) - np.log(minreg)) / (np.log(maxreg) - np.log(minreg))) # Save image to file - filename = user['workingdir'] + 'rego' + os.sep + "rego%06d.png" % user['ncount'] + filename = os.path.join(user['workingdir'],'rego', "rego%06d.png" % user['ncount']) height = user['gridsize'] / dpi width = height fig = plt.figure(figsize=(width, height), dpi=dpi) @@ -223,7 +199,7 @@ def image_shaded_relief(user, DEM): ax.imshow(dem_img, interpolation="nearest", vmin=0.0, vmax=1.0) plt.axis('off') # Save image to file - filename = user['workingdir'] + 'shaded' + os.sep + "shaded%06d.png" % user['ncount'] + filename = os.path.join(user['workingdir'],'shaded',"shaded%06d.png" % user['ncount']) plt.savefig(filename, dpi=dpi, bbox_inches=0) return user From bc52cba3630f374213ecdd062ac29892d6ab254b Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 23 Feb 2022 15:05:05 -0500 Subject: [PATCH 18/19] Finished redirect method to copy files to subfolders and append interval numbers to them --- python/ctem/ctem/driver.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py index 9bda0aa8..7dd8dfef 100644 --- a/python/ctem/ctem/driver.py +++ b/python/ctem/ctem/driver.py @@ -254,7 +254,8 @@ def redirect_outputs(self, filekeys, foldername): for k in filekeys: forig = self.output_filenames[k] - fdist = os.path.join(self.user['workingdir'], foldername, f"{k}_{self.user['ncount']:06d}.dat") + fdest = os.path.join(self.user['workingdir'], foldername, f"{k}_{self.user['ncount']:06d}.dat") + shutil.copy2(forig, fdest) if __name__ == '__main__': sim = Simulation() From ed8cac8b07f5f3f1b6b91a5df2cdc58338b076c0 Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 23 Feb 2022 15:07:59 -0500 Subject: [PATCH 19/19] Streamlined 3d viewer by removing extraneous code --- python/ctem/ctem/viewer3d.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/python/ctem/ctem/viewer3d.py b/python/ctem/ctem/viewer3d.py index 7c68d995..e3d66cb7 100644 --- a/python/ctem/ctem/viewer3d.py +++ b/python/ctem/ctem/viewer3d.py @@ -1,8 +1,6 @@ import numpy as np -import os import open3d import ctem -import os class Polysurface(ctem.Simulation): """A model of a self-gravitating small body""" @@ -10,7 +8,6 @@ def __init__(self): ctem.Simulation.__init__(self, isnew=False) # Initialize the data structures, but doesn't start a new run #used for Open3d module self.mesh = open3d.geometry.TriangleMesh() - self.process_interval(isnew=False) return def ctem2blockmesh(self): @@ -94,7 +91,6 @@ def ctem2trimesh(self): for j in range(s - 1): for i in range(s - 1): idx = (s - 1) * j + i - # faces[idx,:] = [ j*s+i, j*s+i+1, (j+1)*s+i+1, (j+1)*s+i ] faces[idx, :] = [j * s + i, j * s + i + 1, (j + 1) * s + i + 1] idx += (s - 1) ** 2 faces[idx, :] = [(j + 1) * s + i + 1, (j + 1) * s + i, j * s + i ] @@ -112,19 +108,12 @@ def visualize(self): opt = vis.get_render_option() opt.background_color = np.asarray([0, 0, 0]) - # zmax = np.amax(self.surface_dem) - # zmin = np.amin(self.surface_dem) - # cnorm = (self.surface_dem.flatten() - zmin) / (zmax - zmin) - # cval = plt.cm.terrain(cnorm)[:, :3] - self.mesh.paint_uniform_color([0.5, 0.5, 0.5]) vis.run() vis.destroy_window() if __name__ == '__main__': - import matplotlib.pyplot as plt - sim = Polysurface() sim.ctem2trimesh() sim.visualize()