From 766d8a305e450ea92b1e317fc484104696b1c9af Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 17 Feb 2022 12:17:19 -0500 Subject: [PATCH] Moved Python code to correct new location --- python/{ => ctem/ctem}/ctem_draw_surf.py | 0 python/{ => ctem/ctem}/ctem_driver.py | 458 +++++++++---------- python/{ => ctem/ctem}/ctem_io_readers.py | 286 ++++++------ python/{ => ctem/ctem}/ctem_io_writers.py | 526 +++++++++++----------- python/ctem/ctem/ctem_viewer_3d.py | 18 +- python/{ => ctem/ctem}/dist_reader.py | 0 6 files changed, 652 insertions(+), 636 deletions(-) rename python/{ => ctem/ctem}/ctem_draw_surf.py (100%) rename python/{ => ctem/ctem}/ctem_driver.py (97%) rename python/{ => ctem/ctem}/ctem_io_readers.py (97%) rename python/{ => ctem/ctem}/ctem_io_writers.py (97%) rename python/{ => ctem/ctem}/dist_reader.py (100%) diff --git a/python/ctem_draw_surf.py b/python/ctem/ctem/ctem_draw_surf.py similarity index 100% rename from python/ctem_draw_surf.py rename to python/ctem/ctem/ctem_draw_surf.py diff --git a/python/ctem_driver.py b/python/ctem/ctem/ctem_driver.py similarity index 97% rename from python/ctem_driver.py rename to python/ctem/ctem/ctem_driver.py index 712a6157..29057053 100644 --- a/python/ctem_driver.py +++ b/python/ctem/ctem/ctem_driver.py @@ -1,229 +1,229 @@ -#!/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() +#!/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_io_readers.py b/python/ctem/ctem/ctem_io_readers.py similarity index 97% rename from python/ctem_io_readers.py rename to python/ctem/ctem/ctem_io_readers.py index a58ba9ab..84cf7e06 100644 --- a/python/ctem_io_readers.py +++ b/python/ctem/ctem/ctem_io_readers.py @@ -1,144 +1,144 @@ -#!/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 - +#!/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_io_writers.py b/python/ctem/ctem/ctem_io_writers.py similarity index 97% rename from python/ctem_io_writers.py rename to python/ctem/ctem/ctem_io_writers.py index 7d7934e8..7bf21ce4 100644 --- a/python/ctem_io_writers.py +++ b/python/ctem/ctem/ctem_io_writers.py @@ -1,264 +1,264 @@ -#!/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) - +#!/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 index a278239c..237fd99b 100644 --- a/python/ctem/ctem/ctem_viewer_3d.py +++ b/python/ctem/ctem/ctem_viewer_3d.py @@ -95,7 +95,23 @@ def convert_trimesh_to_o3d(self): dem_vec = np.array((xvals, yvals, zvals*10)).T pcd = open3d.geometry.PointCloud() pcd.points = open3d.utility.Vector3dVector(dem_vec) + pcd.estimate_normals() - o3d.visualization.draw_geometries([pcd]) + # estimate radius for rolling ball + distances = pcd.compute_nearest_neighbor_distance() + avg_dist = np.mean(distances) + radius = 1.5 * avg_dist + + mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_ball_pivoting( + pcd, + o3d.utility.DoubleVector([radius, radius * 2])) + + # create the triangular mesh with the vertices and faces from open3d + tri_mesh = trimesh.Trimesh(np.asarray(mesh.vertices), np.asarray(mesh.triangles), + vertex_normals=np.asarray(mesh.vertex_normals)) + + trimesh.convex.is_convex(tri_mesh) + + o3d.visualization.draw_geometries([trimesh]) diff --git a/python/dist_reader.py b/python/ctem/ctem/dist_reader.py similarity index 100% rename from python/dist_reader.py rename to python/ctem/ctem/dist_reader.py