#!/usr/bin/env python # -*- coding: utf-8 -*- # ========================================================================= # ========================================================================= ''' CarboCELL_5-1.py Experimental carbonate Cellular Automaton (CA) population simulator P Burgess RHUL, C Jenkins INSTAAR, D Potts UCSC Begun at CSDMS SanAntonio Conference 16 Oct 2010 Now version 5-1: 31 Oct 2011 at INSTAAR using Python v2.6 and matplotlib 1.0.0 This version gets links the program with events and adjusts the colour bars ''' # ====================== 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 ==================== # ========================================================================= import time,sys,os,glob import random import numpy as NPY import scipy.signal as SS from pylab import * import matplotlib.pyplot as PLT import matplotlib.path as PATH import matplotlib.patches as PATCHES import matplotlib.colors as COL # ========================================================================= def myMAIN(): global gridSZ,runDirS,prjS global tStep,tNow,tTotal,tShow global neighbDepth,centreN, cellScale global taskQ,numFacies,colorsA,hexColrsA global stratThicknsYXa, stratFaciesYXa #Final 3D matrices global winnersA,losersA #tallies of the successful processes each timestep global tallyA global caStatsDataA global mnmxNeighbr,mnmxTrigger global prvAgesA,prvReliefA #(Put in global because it rides through all modules) global nextFacies global winnersA; winnersS="staySame,newGrwth,newClone,newSpawn,newImmi" global losersA; losersS="genrlMort,competdOut,massMort,longevMort,stillBarren" global iNUL iNUL=-99 loudQ=False #Controls the amount of screen output for this module #Banner ------------------------------------------------------------------ print "*"*70;print print "\t"*3,"Program pyCarboCAT_*.py"; print print "\t","A Python program that holds a cellular population model "; print print "\t"," and a version of carboCAT by P Burgess"; print print "\t","Peter Burgess RHUL, Chris Jenkins INSTAAR, Donald Potts UCSC" print "\t","Questions to: chris.jenkins@colorado.edu" print "\t"*4,"Jan 2011" print;print "*"*70;print #Initialize some settings ----------- print; print "Setup ..."; print gridSZ=51 tStep=1; tTotal=100; tIni=-1 #Meant to be ~years stepTotal=tTotal/tStep step=0 #Step counts chronostratigraphic grid layers tA=[tIni] tShow=20 #Steps on which a graphic will be produced # (in addition to initial and final) #Create the project and the output sites -------------- prjS="su" #"Sunday" Clinic 2013 (use no underscore) prjDirS="_"+prjS+"/" try: os.mkdir("_"+prjS) except: pass #Probably just that it already exists runS=time.strftime("%d%b%Y_%H%M", time.localtime()) #runS=time.strftime("%a, %d %b %Y %H:%M:%S +0000", gmtime()) #Eg: 'Thu, 28 Jun 2001 14:17:15 +0000' runDirS=prjDirS+"_"+runS+"/" print "File destination: runDirS=",runDirS try: os.mkdir(runDirS) except: pass #Probably just that it already exists try: os.mkdir(runDirS+"_dumps/") except: pass try: os.mkdir(runDirS+"_correl/") except: pass try: os.mkdir(runDirS+"_pattern/") except: pass try: os.mkdir(runDirS+"_ages/") except: pass print "Outputs will go to... runDirS=",runDirS #Set some basic model configurations ------------------- taskQ="model" #taskQ="autom" #Choose whether carboPOP ("model") or carboCAT ("autom") cellScale=1. #Approx scale of cells (for rugosity only) neighbDepth=2 #0=> 1x1, 1=> 3x3, 2=> 5x5 neighbourhoods #Note: Total number of cells is always odd centreN=((neighbDepth*2+1)**2)/2 #(Index of centre cell in unravelled array) foo=getEvents() #Reads massMortA,massInflxA from file *_events.txt #Set the facies characteristics including build up ... #(Note: at ini this loads in ALL the facies characteristics) foo=caAccumulate(NPY.array([]),"ini") print; print "taskQ,numFacies=",taskQ,numFacies #Make discrete colour scale for numFacies... norm = COL.Normalize(0.,float(numFacies-1)) colorsA = [cm.jet(norm(nF)) for nF in map(float,range(numFacies))] colorsA[0]=(0.5, 0.5, 0.5) #Set (0) to grey... print "colorsA=",colorsA hexColrsA=[] for rgbA in colorsA: hexColrsS=COL.rgb2hex(rgbA[:3]) hexColrsA.append(hexColrsS) #print rgbA,hexColrsS print "colorsA,hexColrsA=",colorsA,hexColrsA #Set the CA rules... if taskQ=="autom": print 'This run is of pyCarboCAT (taskQ="autom")' #These are used if "autom"; "model" will set it's own #(Peter's standard settings are [4,10][6,10] with # 4 lithofacies incl empty) mnmxNeighbr=[5,7] mnmxTrigger=[5,8] elif taskQ=="model": print 'This run is of CarboCELL (taskQ="model")' mnmxNeighbr=[0,0]; mnmxTrigger=[0,0] #Unused else: raw_input("Bad task type ! ("+taskQ+")") #Make a neighbourhood cookie cutter r=NPY.arange(-neighbDepth,neighbDepth+1,+1) #Start at bottom left cutXa,cutYa=NPY.meshgrid(r,r) print "cutYa=","\t".join(map(str,cutYa)) print "cutXa=","\t".join(map(str,cutXa)) #Initial grid ------------------------ iniXYa=[]; tNow=tIni #tallyFaciesA=[0]*numFacies for nY in range(gridSZ): caAy=[] for nX in range(gridSZ): randFacies=int(random()*float(numFacies)) caAy.append(randFacies) #Tally... #tallyFaciesA[randFacies]+=1 iniXYa.append(caAy) iniXYa=NPY.array(iniXYa) #print "iniXYa.shape,tallyFaciesA=",iniXYa.shape,tallyFaciesA #caAge is a mapping of the current age of each cell's occupant prvAgesA=NPY.array([[0]*iniXYa.shape[0]]*iniXYa.shape[1]) #print prvAgesA.shape,prvAgesA newAgesA=prvAgesA[:] #Create operating ca grids ... prvGridA=iniXYa[:] newGridA=iniXYa[:] #Initialize the final 3D strat grids... stratFaciesYXa=NPY.array([[[0]*gridSZ]*gridSZ]*(stepTotal+1)) stratThicknsYXa=NPY.array([[[0.]*gridSZ]*gridSZ]*(stepTotal+1)) #print "Checking: step,stratFaciesYXa.shape=",step,stratFaciesYXa.shape #print " stratFaciesYXa[0]=",stratFaciesYXa[0] #Bathymetry (added by CJJ 31Oct2011) prvReliefA=NPY.array([[0.]*iniXYa.shape[0]]*iniXYa.shape[1]) newReliefA=prvReliefA[:] print; print "Setup completed !" kyS=raw_input(" proceed ? (.y/*n) (Will delete prior plots)") if kyS.lower()<>"y" and kyS<>"": print "Stopping and exiting !" sys.quit() print;print "-"*60 #Delete previous plots... globA=glob.glob(runDirS+"*.*") for globE in globA: os.remove(globE) globA=glob.glob(runDirS+"_dumps/*.*") for globE in globA: os.remove(globE) #Plot and write the initial state... foo=writeOutput(iniXYa,"facies_"+taskQ+"_"+prjS+"_"+str(tIni)) foo=plotFacies(iniXYa,"facies_"+taskQ+"_"+prjS+"_"+str(tIni),0) foo=plotAges(prvAgesA,"ages_"+taskQ+"_"+prjS+"_"+str(tIni),0) foo=plotRelief(prvReliefA,"relief_"+taskQ+"_"+prjS+"_"+str(tIni),0) #foo=plotFFT(iniXYa,"fft_"+taskQ+"_"+prjS+"_"+str(tIni),0) #foo=plotCorrel(iniXYa,"correl_"+taskQ+"_"+prjS+"_"+str(tIni),0) #Gather the statistics on the initial ca ----- caStatsDataA=[[0]]*numFacies statsXYa=caStatistics(iniXYa) for cfN in range(numFacies): caStatsDataA[cfN]=[statsXYa[cfN]] #Do the accumulation sums on the initial grid ----- foo=caAccumulate(newGridA,"dat") #thicknsYXa,faciesYXa come back in global for rY in range(gridSZ): pRAy=prvReliefA[rY] tYXa=thicknsYXa[rY] nRAy=[pRAy[rX]+tYXa[rX] for rX in range(gridSZ)] #print "nRAy=",nRAy newReliefA[rY]=nRAy[:] #print "newReliefA=",newReliefA #Populate with the initial grids ... stratFaciesYXa[step]=faciesYXa[:] #Step here is 0 stratThicknsYXa[step]=thicknsYXa[:] rugosty=rugosity(newReliefA) #Initial rugosity #--------------------------------------------------------- for tNow in range(0,tTotal,tStep+1): #The +1 ensures that the closing step is included if tNow % tShow==0: print "tNow=",tNow step+=1 tA.append(tNow) #Per map... winnersA=[0]*5; losersA=[0]*5 tallyA=[0]*numFacies nextFacies=-1 #ie Resets per map and includes empty ! for nY in range(gridSZ): for nX in range(gridSZ): #Working on cell iniXYa[nY,nX] neighbFaciesA=[] #Is a list - will be formed unravelled ! neighbAgesA=[] for cutI in range(len(cutYa)): for cutJ in range(len(cutXa)): cellY=nY+cutYa[cutI][cutJ] cellX=nX+cutXa[cutI][cutJ] #print "cutI,cutJ,cutYa[cutI][cutJ],cutXa[cutI][cutJ]=",\ # cutI,cutJ,cutYa[cutI][cutJ],cutXa[cutI][cutJ] #Allow for wrap-around... cellY0=cellY; cellX0=cellX if cellY>=gridSZ or cellY<0: cellY=cellY % gridSZ if cellX>=gridSZ or cellX<0: cellX=cellX % gridSZ if cellY<>cellY0 or cellX<>cellX0: #print "Wrapped:" ,cellY,cellX,"(",cellY0,cellX0,")" pass neighbFaciesA.append(prvGridA[cellY,cellX]) neighbAgesA.append(prvAgesA[cellY,cellX]) #print "cellY,cellX=",cellY,cellX #The CA rules implementation --- neighbFaciesA=NPY.array(neighbFaciesA) neighbAgesA=NPY.array(neighbAgesA) '''if neighbFaciesA[4]<>prvGridA[nY][nX]: print neighbFaciesA[4],\ prvGridA[nY][nX]," at [nY][nX]=",nY,nX print "neighbFaciesA=",neighbFaciesA raw_input("Mismatch !")''' if taskQ=="model": newGridA[nY,nX],plusAge=cellModelling(\ neighbFaciesA,neighbAgesA) newAgesA[nY,nX]=plusAge #Chris Jenkins' carboPOP elif taskQ=="autom": newGridA[nY,nX],plusAge=cellAutomaton(\ neighbFaciesA,neighbAgesA) #print "at assign nGa nYnX:",prvGridA[0] newAgesA[nY,nX]=plusAge #Peter Burgess' carboCAT else: raw_input("Unknown task !") if loudQ or tNow % tShow==0: print "winnersA=",winnersA,"(",winnersS,")" print "losersA=",losersA,"(",losersS,")" print "tallyA=",tallyA foo=writeOutput(newGridA,"facies_"+taskQ+"_"+prjS+"_"+str(tNow)) if tNow % tShow==0: print "Plotting at tNow=",tNow plotFacies(newGridA,"facies_"+taskQ+"_"+prjS+"_"+str(tNow),0) foo=plotAges(newAgesA,"ages_"+taskQ+"_"+prjS+"_"+str(tNow),0) #foo=plotFFT(newGridA,"fft_"+taskQ+"_"+prjS+"_"+str(tNow),0) #foo=plotCorrel(newGridA,"correl_"+taskQ+"_"+prjS+"_"+str(tNow),0) foo=plotRelief(newAgesA,"relief_"+taskQ+"_"+prjS+"_"+str(tNow),0) rugosty=rugosity(newReliefA) #Initial rugosity print "Rugosity: rugosty=",rugosty #Gather the statistics on the initial ca ----- statsXYa=caStatistics(newGridA) for cfN in range(numFacies): caStatsDataA[cfN].append(statsXYa[cfN]) #Do the accumulation sums -------------------- foo=caAccumulate(newGridA,"dat") #thicknsYXa,faciesYXa come back in global for rY in range(gridSZ): pRAy=prvReliefA[rY] tYXa=thicknsYXa[rY] nRAy=[pRAy[rX]+tYXa[rX] for rX in range(gridSZ)] #print "nRAy=",nRAy newReliefA[rY]=nRAy[:] stratFaciesYXa[step]=faciesYXa[:] stratThicknsYXa[step]=thicknsYXa[:] #Put the newly calculated grid to the 'previous' grid ... prvGridA=newGridA[::] prvAgesA=newAgesA[::] prvReliefA=newReliefA[::] print "-"*60;print "Finishing: make summary graphics" foo=caAccumulate(NPY.array([]),"end") #In figures, stratThicknsYXa will control block height # and stratFaciesYXa will control block colour #print "stratThicknsYXa=",stratThicknsYXa #print "stratFaciesYXa=",stratFaciesYXa #Now plot profile sections of the stratThicknsYXa matrix ----------------- sectionA=[[0.5,0.5],[0.0,1.0]] #(Sets the cut to make as fraction down/across the map area) foo=plotStratig(sectionA,"strat_"+taskQ+"_"+prjS+"_"+str(tNow),1) foo=plotHistories(tA,caStatsDataA,"stats_"+taskQ+"_"+prjS+"_"+str(tNow),1) #raw_input("End of processing") print "Ended" return # ========================================================================= def rugosity(rlifA): #Calc along Y then also along X #This is an approximation # (better is: "http://www.jennessent.com/downloads/WSB_32_3_Jenness.pdf") #This is written assuming the map and cells are equi-sided #Note also the role of cellScale (approx. absolute cell dimensions) #print "rlifA=",list(rlifA) stepsX=0.;stepsY=0. for rY in range(1,gridSZ): stepsX+=sum([abs(rlifA[rY][rX]-rlifA[rY][rX-1])*cellScale for rX in range(1,gridSZ)]) stepsY+=sum([abs(rlifA[rY][rX]-rlifA[rY-1][rX])*cellScale for rX in range(gridSZ)]) allArea=float((gridSZ*cellScale)**2) stepWiseArea=stepsX+stepsY+allArea rugos=stepWiseArea/allArea "<>writeOutput:" outF=open(runDirS+"_dumps/"+nameS+".txt","w") outF.write("pyCarboCAT run: "+\ ":".join(map(str,time.localtime()[:5]))+\ "at time step: "+str(tNow)+\ "\n") outS="\n".join(map(str,caOutA)) outF.write(outS+"\n") outF.close() return # ========================================================================= def caStatistics(caAnalyzA): loudQ=False if loudQ: print ">>caStatistics:" yFaciesCountA=[0]*numFacies for iY in range(len(caAnalyzA)): caAnalyzY=list(caAnalyzA[iY]) for fN in range(numFacies): yFaciesCountA[fN]+=caAnalyzY.count(fN) #print "<>caAccumulate: statusQ,nowXYa.shape=",\ statusQ,nowXYa.shape #print nowXYa if statusQ=="ini": if taskQ=="autom": numFacies=3+1 faciesNamesA=["Empty"] faciesNamesA.extend("Facies_"+str(n) for n in range(numFacies-1)) print "numFacies,faciesNamesA=",numFacies,faciesNamesA faciesAccumRates=[0.] faciesAccumRates.extend([1.]*(numFacies-1)) faciesSolidity=[0.] faciesSolidity.extend([1.]*(numFacies-1)) elif taskQ=="model": #Get the facies details --------------------------- faciesNamesA=[] faciesAccumRates=[] #ie Growth faciesLongevA=[]; faciesMatureA=[] faciesSolidity=[] #ie 1/collapsability faciesRivalsA=[] #Successful competitors fFacies=open("_facies_"+prjS+".txt","r") if loudQ: print "The facies: (number,name,accum_mm/yr,solidity_*)" while fFacies: fFaciesS=fFacies.readline() if fFaciesS=="": break fFaciesA=fFaciesS.strip().split("\t") if not "Index" in fFaciesS: #Don't do for header print "fFaciesA=",fFaciesA faciesAccumRates.append(float(fFaciesA[2])/1000.) #(/1000 is to convert to metres) faciesSolidity.append(float(fFaciesA[3])) faciesNamesA.append(fFaciesA[1]) faciesLongevA.append(int(fFaciesA[4])) faciesMatureA.append(int(fFaciesA[5])) if fFaciesA[6]<>"[]": rivalsA=map(int,fFaciesA[6].strip('"')\ .strip("[").strip("]").split(",")) else: rivalsA=[] faciesRivalsA.append(rivalsA) fFacies.close() numFacies=len(faciesNamesA) if loudQ: print "faciesNamesA,faciesAccumRates,faciesSolidity=",\ faciesNamesA,faciesAccumRates,faciesSolidity print "faciesMatureA,faciesLongevA,faciesRivalsA=",\ faciesMatureA,faciesLongevA,faciesRivalsA #-------------------------------------------------- elif statusQ=="dat": #Applies to tNow: change stock to carb thickness #print "faciesSolidity=",faciesSolidity thicknsYXa=NPY.array([[0.]*gridSZ]*gridSZ) faciesYXa=NPY.array([[0.]*gridSZ]*gridSZ) #print "thicknsYXa.shape=",thicknsYXa.shape for iY in range(gridSZ): for iX in range(gridSZ): #print "nowXYa[iY,iX]=",nowXYa[iY,iX] thicknsYXa[iY,iX]=float(tStep)\ *faciesAccumRates[nowXYa[iY,iX]]\ *faciesSolidity[nowXYa[iY,iX]] #ie time*growth*solidity faciesYXa[iY,iX]=nowXYa[iY,iX] if loudQ: print "thicknsYXa (mm)=",thicknsYXa print "faciesYXa (#code)=",faciesYXa elif statusQ=="end": pass return #Outputs passed via global # ========================================================================= def getEvents(): global massMortA,massInflxA massMortA=[]; massInflxA=[] '''massMortA=[["bleaching",7,0.3,[1,2,3,5]],\ ["cots",20,0.5,[1,3,5]], ["cyclone",14,0.1,[1,3,4,5]], ["disease",5,0.2,[2]] ] #ie [event-type,period,mortality/1,facies-affected]''' eventsF=open("_events_"+prjS+".txt","r") eventsS=eventsF.readline() #Header while eventsF: eventsS=eventsF.readline() if eventsS=="": break #EOF #print "eventsS=",eventsS evA=eventsS.strip().split("\t") #print "evA=",evA try: mortFaciesA=map(int,evA[3].strip("[").strip("]").split(",")) except: mortFaciesA=[] massMortA.append([evA[0],float(evA[1]),float(evA[2]),mortFaciesA]) #print evA[5].strip("[").strip("]").split(",") try: influxFaciesA=map(int,evA[5].strip("[").strip("]").split(",")) except: influxFaciesA=[] massInflxA.append([evA[0],float(evA[1]),float(evA[4]),influxFaciesA]) #Reading from a time series will be implemented later eventsF.close() print; print "Set up events..." print "massMortA="," ".join(map(str,massMortA)) print "massInflxA="," ".join(map(str,massInflxA)) print return # ========================================================================= def cellModelling(neighbFaciesA,neighbAgesA): global winnersA #staySame,newGrwth,newClone,newSpawn,newImmi global losersA #genMort,comptnMort,massMort,longevMort,stillBarren cellFaciesA=[-1,-1] loudQ=False #print "neighbFaciesA=",neighbFaciesA #Make some counts ----------------------------------------------- neighbFaciesA=list(neighbFaciesA) neighbAgesA=list(neighbAgesA) centreN=((neighbDepth*2+1)**2)/2 centreAge=neighbAgesA[centreN] centreFacies=neighbFaciesA[centreN] centreCount=neighbFaciesA.count(centreFacies) foo=neighbFaciesA.pop(centreN) #(Neighbour array neighbFaciesA is now -1 element, the centre) #Make handy neighbourhood statistics... neighbCountsA=[] for i in range(numFacies): neighbCountsA.append(list(neighbFaciesA).count(i)) if loudQ: print ">>cellModelling: centreFacies,centreAge,centreCount=",\ centreFacies,centreAge,centreCount totalNeighbours=sum(neighbCountsA[1:]) #Number of neighbourhood non-empties dominantNeighb=NPY.argmax(neighbCountsA[1:])+1 #from non-empties dominantCount=neighbCountsA[dominantNeighb] #(Ditto) randomCell=random() #Out of 1. if loudQ: print "cell Rules: tNow,centreFacies,randomCell,=",\ tNow,centreFacies,randomCell if centreFacies<>0: #Occupied cells - prone to mortalities #Is it prone to massMort ? First calc some params ... massMortQ=False for mM in massMortA: #print "mM=",mM,mM[2]*10. if tNow % mM[1]==0 and tNow<>0: if randomCell=faciesLongevA[centreFacies]: #Longevity mortality (occupied cells only) cellFaciesA=[0,0] losersA[3]+=1 if loudQ: print "longevityMort: centreAge,faciesLongevA[centreFacies]=",\ centreAge,faciesLongevA[centreFacies] elif numRivals>=totalNeighbours/2: #Competition ... # the facies is outcompeted when a majority of # neighbours are rivals ... cellFaciesA=[dominantNeighb,0] losersA[1]+=1 if loudQ: print "Competition: (type)",rival #Time-based events...? elif massMortQ: #The massMort events come in through getEvents() #print massMortA cellFaciesA=[0,0] losersA[2]+=1 if loudQ: print "massMort: (type)",mM[0] else: #Cell unchanged ... cellFaciesA=[centreFacies,centreAge+1] if centreFacies==0: losersA[4]+=1 else: winnersA[0]+=1 if loudQ: print "Unchanged: centreFacies,centreAge=",\ centreFacies,centreAge else: #Empty cell #Spawning (local,random) spawnFacies=int(random()*numFacies) if spawnFacies>0: cellFaciesA=[spawnFacies,0] #(Empty facies cant be spawned) if loudQ: print "Spawn: cellFaciesA(new)=",cellFaciesA winnersA[2]+=1 elif dominantCount>=totalNeighbours/2: #Opportunity for growth or dominant facies (organism) # (provided more than 2 of max one in neighb) #growth into empty cells cellFaciesA=[dominantNeighb,0] #print "Grwth: cellFaciesA,dominantCount=",\ # cellFaciesA,dominantCount if loudQ: print "Grwth instance" winnersA[1]+=1 else: #Cell unchanged ... cellFaciesA=[0,centreAge+1] if centreFacies==0: losersA[4]+=1 else: winnersA[0]+=1 if loudQ: print "Unchanged: centreFacies,centreAge=",\ centreFacies,centreAge if loudQ: print "cellFaciesA",cellFaciesA if -1 in cellFaciesA: raw_input("Problem at cellFaciesA") return cellFaciesA # ========================================================================= def cellAutomaton(neighbFaciesA,neighbAgesA): global winnersA #staySame,newGrwth,newClone,newSpawn,newImmi global losersA #genMort,comptnMort,massMort,longevMort,stillBarren global tallyA global nextFacies #Cycles through facies for recruitments loudQ=False #This is Peter Burgess' original carboCAT # (Note: Settings are made in MAIN) fateA=[""]*2 #Make some counts ----------------------------------------------- neighbFaciesA=list(neighbFaciesA.ravel()) neighbAgesA=list(neighbAgesA.ravel()) centreFacies=neighbFaciesA[centreN] #(The facies number of the centre cell) centreCount=neighbFaciesA.count(centreFacies) #(How many of the crentre facies are in the neighbood) centreAge=neighbAgesA[centreN] '''#Testing... #print "centreN,neighbFaciesA,centreFacies=",\ # centreN,neighbFaciesA,centreFacies cellFaciesA=[centreFacies,centreAge] return cellFaciesA''' foo=neighbFaciesA.pop(centreN) #Note: Neighbour array neighbFaciesA is now minus element, the centre #Make handy neighbourhood statistics... neighbCountsA=[] for i in range(numFacies): neighbCountsA.append(list(neighbFaciesA).count(i)) #print "neighbFaciesA,neighbCountsA=",neighbFaciesA,neighbCountsA if loudQ: print; print ">>cellAutomaton: centreFacies,neighbFaciesA,centreAge=",\ centreFacies,neighbFaciesA,centreAge if centreFacies==0: #Recruitment into empties ... #print "max(neighbCountsA[1:])=",max(neighbCountsA[1:]) neighbTestA=[0]*numFacies nT=-1 for nC in neighbCountsA: nT+=1 if nC>=mnmxTrigger[0] and nC=mnmxNeighbr[0] and \ neighbCountsA[centreFacies]>plotHistories: caStatsDataA=",caStatsDataA #Will use use colorsA computed at start... f4 = PLT.figure() #Set graph box to make space for legends ax1=f4.add_axes([0.1,0.3,0.8,0.6]) #l,b,w,h cfN=len(statsDatA) #Number of species, etc labelA=['']*cfN lineA=['']*cfN for cfL in range(cfN): #print "cfL,tA,colorsA[cfL],hexColrsA[cfL]=",\ # cfL,tA,colorsA[cfL],hexColrsA[cfL] lineA[cfL]=PLT.plot(tA, statsDatA[cfL], hexColrsA[cfL], \ label=tA[cfL]) PLT.plot(0.,1.,hexColrsA[-1]) #Sets constant graph extent #Legend ------------------------- leg=f4.legend(lineA,faciesNamesA,'lower left',labelsep=0) for t in leg.get_texts(): t.set_fontsize(8) # the legend text fontsize PLT.title(nameS+"_"+"".join(map(str,time.localtime()[:5]))) f4.savefig(runDirS+nameS+'.png') if drawQ==1: PLT.draw() if loudQ: print "<>plotStratig: sectionA=",sectionA #prepare the plot space... f2 = PLT.figure() ax = f2.add_subplot(111) #Cut the stratal matrices ------------------------------------------ #(stratThicknsYXa,stratFaciesYXa come in via global) sectnCut=int(sectionA[0][0]*gridSZ) faciesProfileA=stratFaciesYXa[:,sectnCut,:] thicknsProfileA=stratThicknsYXa[:,sectnCut,:] profileShapeA=faciesProfileA.shape #[T,x] if loudQ: print "Checking: sectnCut,profileShapeA=",sectnCut,profileShapeA #print " faciesProfileA=",faciesProfileA #print " thicknsProfileA=",thicknsProfileA #------------------------------------------------------------------- #Follows the pplot auto-colours maxT=0. faciesLeg=[[],[]] patches4legA=[] tChron=tShow #Will draw time line each tChron time steps (usu =tShow) chronVerts=[[] for i in range(tTotal/tChron)] #Note: chron=0 is the first (non-zero) one #Change the matrix into a structure of the box shapes... for spX in range(profileShapeA[1]): #for x-width (ie work on each column) fspX=float(spX) tBot=0. tTop=0. #thicknsProfileA[-1,spX] #Bottom cell in xth column #print "spX,fspX,tBot,tTop=",spX,fspX,tBot,tTop #print " spT,spTT,tBot,tTop,maxT,faciesNum,colorsA[faciesNum]=" for spT in range(profileShapeA[0]): #for y-height (and start from bottom) #print "fspT,fspX",fspT,spTT,fspX tBot=tTop #Previous value tTop+=thicknsProfileA[spT,spX] maxT=max(maxT,tTop) faciesNum=faciesProfileA[spT,spX] if tTop==tBot: #No thickness alteration pass else: #print " ",spT,spT,tBot,tTop,maxT,\ # faciesNum,colorsA[faciesNum] #Draw the main patchwork... verts = [ (fspX, tBot), # left, bottom (fspX, tTop), # left, top (fspX+1., tTop), # right, top (fspX+1., tBot), # right, bottom (fspX, tBot), # (Closure) ignored ] codes = [PATH.Path.MOVETO, PATH.Path.LINETO, PATH.Path.LINETO, PATH.Path.LINETO, PATH.Path.CLOSEPOLY, ] caPath = PATH.Path(verts, codes) #The drawing... patch=PATCHES.PathPatch(caPath, \ facecolor=colorsA[faciesNum], lw=1) ax.add_patch(patch) #Add to the legend lists if necessary... if not faciesNum in faciesLeg[0]: patches4legA.append(patch) faciesLeg[0].append(faciesNum) faciesLeg[1].append(faciesNamesA[faciesNum]) #Every so often add a white time line... if spT % tChron==0 and spT<>0: chron=spT/tChron-1 #Index for the chron #print "do tChron line: spT,tTop,chron=",spT,tTop,chron chronVerts[chron].append([fspX, tTop]) chronVerts[chron].append([fspX+1., tTop]) #Draw the chron lines... #print "chronVerts=",chronVerts chronVerts=NPY.array(chronVerts) for chronV in chronVerts: #print "chronV=",chronV PLT.plot(chronV[:,0], chronV[:,1], lw=1, color='w') #Fix the axes post-hoc... ax.set_xlim(0,profileShapeA[1]) ax.set_ylim(0,maxT*1.2) #Quikfix: The 1.2 makes room for the legend PLT.title(nameS+"".join(map(str,time.localtime()[:4]))) #Legend ------------------------- leg=PLT.legend(patches4legA,faciesLeg[1],loc='upper left',labelsep=0) for t in leg.get_texts(): t.set_fontsize(8) # the legend text fontsize f2.savefig(runDirS+nameS+'.png') if drawQ==1: PLT.draw() if loudQ: print "<>plotAges: nameS,drawQ=",nameS,drawQ #print plotXYa.shape,plotXYa plotXYa[0,0]=0 plotXYa[-1,-1]=tTotal #Start it... f1 = PLT.figure() #First make the coord and flattened data arrays... r=NPY.arange(0.,float(gridSZ)+.01,1.) X,Y=NPY.meshgrid(r,r) PLT.pcolor(X,Y,plotXYa) #PLT.scatter(nA, mA, c=nmA, s=symSz, label='n m', faceted=False) PLT.colorbar() PLT.title(nameS+" "+"".join(map(str,time.localtime()[:4]))) f1.savefig(runDirS+"_ages/"+nameS+'.png') if drawQ==1: PLT.draw() #print "<>plotRelief: nameS,drawQ=",nameS,drawQ #print plotXYa.shape,plotXYa plotXYa[0,0]=0 plotXYa[-1,-1]=tTotal #Start it... f1 = PLT.figure() #First make the coord and flattened data arrays... r=NPY.arange(0.,float(gridSZ)+.01,1.) X,Y=NPY.meshgrid(r,r) PLT.pcolor(X,Y,plotXYa) #PLT.scatter(nA, mA, c=nmA, s=symSz, label='n m', faceted=False) PLT.colorbar() PLT.title(nameS+" "+"".join(map(str,time.localtime()[:4]))) f1.savefig(runDirS+"_pattern/"+nameS+'.png') if drawQ==1: PLT.draw() #print "<>plotFacies: nameS,drawQ=",nameS,drawQ toPlotXYa[0,0]=numFacies-1 toPlotXYa[-1,-1]=0 #Start it... f1 = PLT.figure() #First make the coord and flattened data arrays... r=NPY.arange(0.,float(gridSZ)+.01,1.) X,Y=NPY.meshgrid(r,r) PLT.pcolor(X,Y,toPlotXYa,cmap=cm.jet) #PLT.scatter(nA, mA, c=nmA, s=symSz, label='n m', faceted=False) PLT.colorbar() PLT.title(nameS+" "+"".join(map(str,time.localtime()[:4]))) f1.savefig(runDirS+"_pattern/"+nameS+'.png') if drawQ==1: PLT.draw() #print "<