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/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/__init__.py b/python/ctem/ctem/__init__.py index 4774c6af..165a1bb0 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.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 84cf7e06..00000000 --- a/python/ctem/ctem/ctem_io_readers.py +++ /dev/null @@ -1,144 +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 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']=float(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 ('sfdcompare' == fields[0].lower()): parameters['sfdcompare']=fields[1] - if ('interval' == fields[0].lower()): parameters['interval']=float(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'] = float(fields[1]) - parameters['shadedminhdefault'] = 0 - if ('shadedmaxh' == fields[0].lower()): - parameters['shadedmaxh'] = float(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'] = float(fields[0]) - parameters['ncount'] = int(fields[1]) - parameters['curyear'] = float(fields[2]) - parameters['restart'] = fields[3] - parameters['fracdone'] = float(fields[4]) - parameters['masstot'] = float(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]) - 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 = float(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/ctem_viewer_3d.py b/python/ctem/ctem/ctem_viewer_3d.py deleted file mode 100644 index be967d10..00000000 --- a/python/ctem/ctem/ctem_viewer_3d.py +++ /dev/null @@ -1,92 +0,0 @@ -import numpy as np -import os -import open3d -import ctem_io_readers -import os - - -class polysurface: - """A model of a self-gravitating small body""" - def __init__(self): - self.polyfile = "" - #used for Open3d module - self.mesh = open3d.geometry.TriangleMesh() - return - - def ctem2mesh(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, - '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 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']) - # Build mesh grid - s = parameters['gridsize'] - pix = 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.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] - idx += (s - 1) ** 2 - faces[idx, :] = [j * s + i + 1, (j + 1) * s + i + 1, (j + 1) * s + i] - self.mesh.vertices = open3d.utility.Vector3dVector(verts) - self.mesh.triangles = open3d.utility.Vector3iVector(faces) - self.mesh.compute_vertex_normals() - - -if __name__ == '__main__': - import matplotlib.pyplot as plt - - surf = polysurface() - surf.ctem2mesh(surface_file="surface_dem.dat") - - 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]) - - 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..7dd8dfef --- /dev/null +++ b/python/ctem/ctem/driver.py @@ -0,0 +1,262 @@ +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", isnew=True): + 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) + + + # Starting new or old run? + 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 + + # Write ctem.dat file + io.write_datfile(self.user, self.output_filenames['dat'], self.seedarr) + else: + self.user['ncount'] = 0 + 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['impmass']) + + # Read ctem.dat file + 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: + 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] + fdest = os.path.join(self.user['workingdir'], foldername, f"{k}_{self.user['ncount']:06d}.dat") + shutil.copy2(forig, fdest) + +if __name__ == '__main__': + sim = Simulation() + sim.run() \ 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..d0d383e5 --- /dev/null +++ b/python/ctem/ctem/io.py @@ -0,0 +1,382 @@ +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 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 = os.path.join(user['workingdir'], directory) + if not os.path.isdir(dir_test): + os.makedirs(dir_test) + + return + + +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 = (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]]) + + # 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] * user['curyear'] / user['interval'] + + # Create sdist bin factors, which contain one crater per bin + area = (user['gridsize'] * user['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" % 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 = 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) + + # 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(user, DEM): + dpi = 300.0 # 72.0 + pix = user['pix'] + gridsize = user['gridsize'] + ve = 1.0 + 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) + + # 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 = os.path.join(user['workingdir'],'surf',"surf%06d.png" % user['ncount']) + plt.savefig(filename, dpi=dpi, bbox_inches=0) + + return + + +def image_regolith(user, regolith): + # Create scaled regolith image + minref = user['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 = 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) + fig.figimage(regolith_scaled, cmap=cm.nipy_spectral, origin='lower') + plt.savefig(filename) + + return + + +def image_shaded_relief(user, DEM): + dpi = 300.0 # 72.0 + pix = user['pix'] + gridsize = user['gridsize'] + ve = 1.0 + mode = 'overlay' + 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 (user['shadedminh'] > user['shadedmaxh']): + temp = user['shadedminh'] + user['shadedminh'] = user['shadedmaxh'] + user['shadedmaxh'] = temp + else: + user['shadedminh'] = user['shadedminh'] + user['shadedmaxh'] = user['shadedmaxh'] + + # 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 = user['shadedminh'] + if (user['shadedmaxhdefault'] == 1): + shadedmaxh = np.amax(DEM) + else: + shadedmaxh = user['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 = os.path.join(user['workingdir'],'shaded',"shaded%06d.png" % user['ncount']) + plt.savefig(filename, dpi=dpi, bbox_inches=0) + return user + + +def read_datfile(user, datfile, seedarr): + # Read and parse ctem.dat file + + # Read ctem.dat file + print('Reading input file ' + 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: + 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) + index = 1 + while (index < nlines): + fields = lines[index].split() + seedarr[index - 1] = real2float(fields[0]) + index += 1 + + user['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_user_input(user): + # Read and parse ctem.in file + inputfile = user['ctemfile'] + + # Read ctem.in file + print('Reading input file ' + user['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()): 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()): + user['shadedminh'] = real2float(fields[1]) + user['shadedminhdefault'] = 0 + if ('shadedmaxh' == fields[0].lower()): + user['shadedmaxh'] = real2float(fields[1]) + user['shadedmaxhdefault'] = 0 + + # Test values for further processing + if (user['interval'] <= 0.0): + print('Invalid value for or missing variable INTERVAL in ' + inputfile) + if (user['numintervals'] <= 0): + print('Invalid value for or missing variable NUMINTERVALS in ' + inputfile) + if (user['pix'] <= 0.0): + print('Invalid value for or missing variable PIX in ' + inputfile) + if (user['gridsize'] <= 0): + print('Invalid value for or missing variable GRIDSIZE in ' + inputfile) + if (user['seed'] == 0): + print('Invalid value for or missing variable SEED in ' + inputfile) + if (user['sfdfile'] is None): + print('Invalid value for or missing variable SFDFILE in ' + inputfile) + if (user['impfile'] is None): + print('Invalid value for or missing variable IMPFILE in ' + inputfile) + if (user['popupconsole'] is None): + print('Invalid value for or missing variable POPUPCONSOLE in ' + inputfile) + if (user['saveshaded'] is None): + print('Invalid value for or missing variable SAVESHADED in ' + inputfile) + if (user['saverego'] is None): + print('Invalid value for or missing variable SAVEREGO in ' + inputfile) + if (user['savepres'] is None): + print('Invalid value for or missing variable SAVEPRES in ' + inputfile) + if (user['savetruelist'] is None): + print('Invalid value for or missing variable SAVETRUELIST in ' + inputfile) + if (user['runtype'] is None): + print('Invalid value for or missing variable RUNTYPE in ' + inputfile) + if (user['restart'] is None): + print('Invalid value for or missing variable RESTART in ' + inputfile) + + return user + + +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 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_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 % user) + + # Write random number seeds to the file + for index in range(user['seedn']): + fp.write("%12d\n" % seedarr[index]) + + fp.close() + + return + + +def write_production(user, production): + filename = user['sfdfile'] + np.savetxt(filename, production, fmt='%1.8e', delimiter=' ') + + return + diff --git a/python/ctem/ctem/viewer3d.py b/python/ctem/ctem/viewer3d.py new file mode 100644 index 00000000..e3d66cb7 --- /dev/null +++ b/python/ctem/ctem/viewer3d.py @@ -0,0 +1,122 @@ +import numpy as np +import open3d +import ctem + +class Polysurface(ctem.Simulation): + """A model of a self-gravitating small body""" + 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() + return + + def ctem2blockmesh(self): + # Build mesh grid + 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 + 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(), + 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 + 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[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 ] + 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(f2) + self.mesh.compute_vertex_normals() + self.mesh.compute_triangle_normals() + + return + + def ctem2trimesh(self): + # Build mesh grid + 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.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): + 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] + 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() + vis.add_geometry(self.mesh) + opt = vis.get_render_option() + opt.background_color = np.asarray([0, 0, 0]) + + self.mesh.paint_uniform_color([0.5, 0.5, 0.5]) + vis.run() + vis.destroy_window() + + +if __name__ == '__main__': + sim = Polysurface() + sim.ctem2trimesh() + sim.visualize() + + + 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..1212193e --- /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.user.items(): + print(k, v) \ No newline at end of file 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