From 37e880264319d2fc6d46d86d5ddec29f50859865 Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Tue, 13 Jun 2023 13:53:04 -0400 Subject: [PATCH] Added my single-run analysis script; currently separate from CTEM package (also CTEM works without fractions) --- python/analyze_single_run.py | 317 +++++++++++++++++++++++++++++++++++ 1 file changed, 317 insertions(+) create mode 100644 python/analyze_single_run.py diff --git a/python/analyze_single_run.py b/python/analyze_single_run.py new file mode 100644 index 00000000..8ba87db0 --- /dev/null +++ b/python/analyze_single_run.py @@ -0,0 +1,317 @@ +import numpy as np +import matplotlib.pyplot as plt +import ctem +import matplotlib.patches as patches +import matplotlib.colors as colors +import matplotlib.image as mpimg + + +def find_apollo_coords(gridsize, pix): + a14xm = -529810 + a14ym = -110559 + a15xm = 110165 + a15ym = 792410 + a16xm = 470012 + a16ym = -272122 + a17xm = 933110 + a17ym = 612289 + + a14xo = a14xm / pix + a14yo = a14ym / pix + a15xo = a15xm / pix + a15yo = a15ym / pix + a16xo = a16xm / pix + a16yo = a16ym / pix + a17xo = a17xm / pix + a17yo = a17ym / pix + + center = int(gridsize / 2) + + a14xp = center + int(a14xo) + a14yp = center + int(a14yo) + a15xp = center + int(a15xo) + a15yp = center + int(a15yo) + a16xp = center + int(a16xo) + a16yp = center + int(a16yo) + a17xp = center + int(a17xo) + a17yp = center + int(a17yo) + + axs = [a14xp, a15xp, a16xp, a17xp] + ays = [a14yp, a15yp, a16yp, a17yp] + + return axs, ays + +def traverse(x, y, depth, rego): + """finds number of layers to reach depth""" + n = 1 + nm = 0 + while nm < depth: + i = -n + layerdepth = rego[x,y][i] + nm = nm + layerdepth + n += 1 + n = n-1 + partial = nm - depth + thickness = layerdepth - partial + return n, thickness + +def find_melt_at_pixel(x,y,depth,rego,melt,meltdist): + nlayers, bothickness = traverse(y,x,depth,rego) + meltarray = np.zeros((nlayers, nlist)) + nm = np.zeros((nlayers, nlist)) + rt = np.zeros(nlayers) + rt[0] = bothickness + if nlayers > 1: + rt[1:] = rego[y,x][-nlayers+1:] + for i in range(1,nlayers+1): + if i == 1: + meltarray[i-1] = meltdist[y,x][-nlist:] + else: + meltarray[i-1] = meltdist[y,x][-(i * nlist):-(nlist * (i-1))] + nm[i-1] = meltarray[i-1] * rt[-i] + meltarray[nlayers-1] = meltarray[nlayers-1] * (rt[0]/rego[y,x][-nlayers]) + totmelt = np.zeros(nlist) + for k in range(nlist): + totmelt[k] = np.sum(meltarray[:,k]) + return totmelt, np.sum(rt) + + +def aggregate(x,y,depth,n,rego,melt,meltdist,filter=False,apollo=False): + if n==0: + meltarray,rtlist = find_melt_at_pixel(x,y,depth,rego,melt,meltdist) + elif n<0: + print("number of pixels must be >=0") + return + else: + meltarray = [] + rtlist = [] + k = 0 + for i in range(-n,n+1): + for j in range(-n,n+1): + meltarray.append(np.zeros(nlist)) + m,rt = find_melt_at_pixel(x+i,y+j,depth,rego,melt,meltdist) + meltarray[k] = m + rtlist.append(rt) + k += 1 + + meltarray = np.asarray(meltarray) + totmelt = np.sum(meltarray) + + norm = np.zeros(nlist) + totm = np.zeros(nlist) + for l in range(nlist): + if n !=0: + totm[l] = np.sum(meltarray[:,l]) + norm[l] = totm[l] / totmelt + else: + totm[l] = meltarray[l] + norm[l] = totm[l] / totmelt + mf = totmelt/(np.sum(rtlist)*pix*pix) + + #setup data to plot + nonzerocraters = [] + nonzeromelt = [] + nonzerocraternames = [] + sigcraters = [] + sigmelt = [] + nonsigmelt = [] + sigcraternames = [] + for l in range(nlist): + if filter: + if norm[l] > 1e-4: + nonzerocraters.append(str(l)) + nonzeromelt.append(norm[l]) + nonzerocraternames.append(craters[l]) + if norm[l] < 0.05: + nonsigmelt.append(norm[l]) + else: + sigcraters.append(str(l)) + sigmelt.append(norm[l]) + sigcraternames.append(craters[l]) + + else: + if norm[l] > 0: + nonzerocraters.append(str(l)) + nonzeromelt.append(norm[l]) + nonzerocraternames.append(craters[l]) + if norm[l] < 0.05: + nonsigmelt.append(norm[l]) + else: + sigcraters.append(str(l)) + sigmelt.append(norm[l]) + sigcraternames.append(craters[l]) + + if len(nonsigmelt) > 0: + sigcraters.append('%i Others'%len(nonsigmelt)) + sigmelt.append(np.sum(nonsigmelt)) + sigcraternames.append('%i Others'%len(nonsigmelt)) + + #plot a pie chart + + f1, (a1, a2) = plt.subplots(1,2) + f1.suptitle('Pixel: [' + '{:3d}'.format(x) + ',' + '{:3d}'.format(y) + '] (Aggregated to ' + str(n) + ' pixels) (Total melt fraction = ' + '{0:.4f}'.format(mf) + ') (depth = ' + '{0:.2f}'.format(depth) + ' m)') + a1.pie(sigmelt,labels=sigcraternames,startangle=90,autopct='%1.1f%%',shadow=False,colors = iter(plt.cm.Set3(np.linspace(0,1,len(sigmelt))))) + a1.axis('equal') + a2.imshow(img) + #a2.scatter(x, gridsize-y, color='blue', label='Selected Pixel') + a2.xaxis.set_ticklabels([]) + a2.yaxis.set_ticklabels([]) + if n == 0: + a2.scatter(x, gridsize-y, color='blue', label='Selected Region') + else: + rect = patches.Rectangle((x,gridsize-y),n,n,linewidth=1,edgecolor='b',facecolor='b',label='Selected Region') + a2.add_patch(rect) + if apollo: + a2.scatter(apollo_x, surfcoords, color='red', label='Apollo 14-17 Landing Sites') + a2.tick_params(left = False) + a2.tick_params(bottom = False) + a2.legend() + plt.show() + + plt.rcParams['figure.figsize'] = [13, 5] + + return nonzerocraters, nonzerocraternames, nonzeromelt + +def map_crater_to_depth(crater,depth,rego,meltdist,melt,apollo=True,log=False): + if crater > nlist: + print("Please specify a crater number between 1 and " + str(nlist)) + else: + cratermelt = np.zeros((gridsize,gridsize)) + cratermf = np.zeros((gridsize,gridsize)) + for x in range(gridsize): + for y in range(gridsize): + nlayers, bothickness = traverse(y,x,depth) + meltarray = np.zeros((nlayers, nlist)) + nm = np.zeros((nlayers, nlist)) + rt = np.zeros(nlayers) + rt[0] = bothickness + m = np.zeros(nlayers) + if nlayers > 1: + rt[1:] = rego[y,x][-nlayers+1:] + mt = melt[y,x][-nlayers:] + mt[0] = mt[0] * (depth/rego[y,x][-nlayers]) + for i in range(1,nlayers+1): + if i == 1: + meltarray[i-1] = meltdist[y,x][-nlist:] + else: + meltarray[i-1] = meltdist[y,x][-(i * nlist):-(nlist * (i-1))] + nm[i-1] = meltarray[i-1] * rt[-i] + meltarray[nlayers-1] = meltarray[nlayers-1] * (rt[0]/rego[y,x][-nlayers]) + tot = np.zeros(nlayers) + for k in range(nlayers): + m[k] = meltarray[k, crater] + cratermelt[y,x] = np.sum(m) + cratermf[y,x] = cratermelt[y,x] / (depth * pix * pix) + + f1 = plt.figure(1) + a1 = f1.add_subplot(111) + if crater == nlist-1: + a1.set_title('Local melt to a depth of ' + str(depth) + ' m' ) + else: + a1.set_title('Melt from crater ' + craters[crater] + ' to a depth of ' + str(depth) + ' m' ) + if not log: + md = a1.imshow(cratermelt,origin='lower')#,vmin=0.0,vmax=1.0) + else: + md = a1.imshow(cratermelt,origin='lower', norm=colors.LogNorm())#vmin=1e-3))#,vmax=1.0)) + if crater == nlist-1: + cbar = f1.colorbar(md, label='Volume of total melt from local craters') + else: + cbar = f1.colorbar(md, label='Volume of total melt from this crater') + if apollo: + a1.scatter(apollo_x, apollo_y, color='red', label='Apollo 14-17 Landing Sites') + plt.legend() + f2 = plt.figure(2) + a2 = f2.add_subplot(111) + if crater == nlist-1: + a2.set_title('Local melt to a depth of ' + str(depth) + ' m' ) + else: + a2.set_title('Melt from crater ' + craters[crater] + ' to a depth of ' + str(depth) + ' m' ) + if not log: + md = a2.imshow(cratermf,origin='lower')#,vmin=0.0,vmax=1.0) + else: + md = a2.imshow(cratermf,origin='lower', norm=colors.LogNorm())#vmin=1e-3))#,vmax=1.0)) + if crater == nlist-1: + cbar = f2.colorbar(md, label='Fraction of total melt from local craters') + else: + cbar = f2.colorbar(md, label='Fraction of total melt from this crater') + if apollo: + a2.scatter(apollo_x, apollo_y, color='red', label='Apollo 14-17 Landing Sites') + plt.legend() + + plt.show() + + return + +def map_meltfrac_to_depth(depth,apollo=True,log=False): + mf = np.zeros((gridsize,gridsize)) + for x in range(gridsize): + for y in range(gridsize): + nlayers, bothickness = traverse(y,x,depth) + mfarray = np.zeros(nlayers) + rt = np.zeros(nlayers) + if nlayers > 1: + rt[1:] = rego[y,x][-nlayers+1:] + rt[0] = bothickness + mt = melt[y,x][-nlayers:] + mt[0] = mt[0] * (rt[0]/rego[y,x][-nlayers]) + mf[y,x] = np.sum(mt) / (depth * pix * pix) + else: + mf[y,x] = ((melt[y,x][-1] * (depth/rego[y,x][-1])) / (depth * pix * pix)) + + f1 = plt.figure(1) + a1 = f1.add_subplot(111) + if not log: + md = a1.imshow(mf,origin='lower')#,vmin=0.0,vmax=1.0) + else: + md = a1.imshow(mf,origin='lower', norm=colors.LogNorm())#vmin=1e-3))#,vmax=1.0)) + cbar = f1.colorbar(md, label='Total melt fraction') + if apollo: + a1.scatter(apollo_x, apollo_y, color='red', label='Apollo 14-17 Landing Sites') + plt.legend() + plt.title('Total melt fraction to a depth of ' + str(depth) + ' m') + plt.show() + + return + +if __name__ == '__main__': + path = '/Users/owner/Documents/git/CTEM/examples/global-lunar-bombardment/' + #path = '/Users/owner/Desktop/debug/mftest/' + gridsize = 500 + pix = 12.32e3 + pixkm = pix / 1000 + + apollo_x, apollo_y = find_apollo_coords(gridsize,pix) + surfcoords = np.zeros(len(apollo_y)) + for i in range(len(apollo_y)): + surfcoords[i] = gridsize - apollo_y[i] + + rego = ctem.util.read_linked_list_binary(path + '/surface_rego.dat', gridsize) + melt = ctem.util.read_linked_list_binary(path + '/surface_melt.dat', gridsize) + stack = ctem.util.read_unformatted_binary(path + '/surface_stacknum.dat', gridsize, kind='I4B') + meltdist = ctem.util.read_linked_list_binary(path + '/surface_meltdist.dat', gridsize, kind='SP') + + surfimg = path + '/surf/surf000001.png' + img = mpimg.imread(surfimg) + + #nlist = 117 + # craters = ["South Pole-Aitken", "Fecunditatis", "Australe North", "TOPO-13", "Medii", "Bartels-Voskresenskiy", + # "Serenitatis North", "Lamont", "Vaporum", "Fowler-Charlier", "Mutus-Vlacq", "Asperitatis", "Aestuum", + # "Copernicus-H", "Balmer-Kapetyn", "Cruger-Sirsalis", "Dirichlet-Jackson", "Coloumb-Sarton", + # "Schiller-Zucchius", "Smythii", "Amundsen-Ganswindt", "Nubium", "Rupes Recta", "Deslandres", "Lorentz", + # "Pozcobutt", "Fitzgerald-Jackson", "Poincare", "Wegener-Winlock", "Orientale Southwest", "TOPO-22", + # "Galois", "Schickard", "Fermi", "Ingenii", "Gagarin", "Keeler West", "Birkhoff", "Harkhebi", "Campbell", + # "Serenitatis", "Freundlich-Sharanov", "Landau", "Von Karman M", "Leibnitz", "Apollo", "Milne", + # "Grimaldi", "Pasteur", "Nectaris", "Mendel-Rydberg", "Moscoviense North", "Moscoviense", "Korolev", + # "Mendeleev", "Planck", "Hertzsprung", "Humorum", "Crisisum", "Humboldtianum", "Sikorsky-Rittenhaus", + # "Bel'Kovich", "Bailly", "d'Alembert", "Oppenheimer", "Clavius", "Crisium East", "SChwarzschild", + # "La Condamine", "Maupertuis", "Imbrium", "Schrodinger", "Ukert", "Campanus", "Protagoras", "Calippus", + # "Babbage A", "Prinz", "Lambert R", "Egede", "Fontenelle", "T. Mayer", "Stadius", "Cassini", "Letronne", + # "Cruger", "Marius", "Arzachel", "T. Mayer W", "Wallace", "Kepler D", "Orientale", "Lansberg", "Humboldt", + # "Kies", "Timaeus", "Bonpland D", "Krieger", "Lansberg C", "Briggs", "Archimedes", "Thebit", "Hansteen", + # "Davy", "Markov", "Iridum", "Seleucus", "Cardanus", "Krafft", "Dalton", "St. George", "Bianchinni", + # "Herotodus", "Plato", "Damoiseau", "Mairan", "Local"] + + nlist = 8 + craters = ['South Pole-Aitken','Serenitatis','Nectaris','Crisium','Imbrium','Orientale','Iridum','Local'] + + aggregate(apollo_x[2],apollo_y[2],5,0,rego,melt,meltdist)