From cee986627317452a41fe8bb383a00effce5e5457 Mon Sep 17 00:00:00 2001 From: David Minton Date: Tue, 22 Feb 2022 11:54:27 -0500 Subject: [PATCH] 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