# -*- coding: utf-8 -*- # ============================================================================ # ============================================================================ ''' carboLOT_*.py Prototype of a population simulation of carbonate production based on Lotke-Volterra Competition methods Chris Jenkins INSTAAR Donald Potts UCSC Peter Burgess RHUL Started May 2010 as 'numericalCARBONATE' NSF supported EAGER award from Oct 2010 Version 16: 27 Oct 2011 CJJ at BoCO - Revised time base, burst and reporting - Now writing and pick-up from dump files - Changed folders for outputs, put units on graphics - Wave climate - but not completely (needs Snell's Law) - 1-cell selvage around the map - better handling of append to arrays - setup verification ''' # ========================= 16 16 16 16 16 16 16 16 ========================== # ============================================================================ import time,os,glob import sys #Tame the decimal points ? #sys.float_output_precision = 4 #import scipy as SCY import scipy.integrate as scint import numpy as NPY #from numpy import MSK #For masking import struct as STR import matplotlib as MPL import matplotlib.pyplot as PLT import matplotlib.path as PATH import matplotlib.patches as PATCHES import matplotlib.colors as COL #import userPORTS #Ready for users to put modules in global kbFileName kbFileName="_carboKB_7.okb" #Determines the KnowledgeBase eg "KB_3.txt" #(Concievably, users may substitute their own) loudQ=False #Determines whether run screen is verbose #(Not a global: different each module)iCr,kbWdHabittA # ============================================================================ # ============================================================================ ''' Structure: getKBdata() Reads the organism-knowledegebase into arrays reconcileValuesKB(kbDataA) Tidies up the KB data as it proceeds into arrays: brings human-friendly syntaxes into machine-syntaxes interpHabiRange(enviroPrpty,habitatS) Calculates a organism parameter (eg suitability) from an enviro value of a grid cell by interpolating to a position in a linear set from the KnowledgeBase. getELEVN(lolaA) Finds the bathymetry over the whole work area and loads it into elevnA. lolaA here is the central position of the square grid (presently modlSZx*modlSZy). get3x3elevn(lolaSEEKa) The code for interrogating the GEBCO 1nm global bathymetry data grid. Outputs elevn (interps 4 closest grid points), and slope and roughness (based on 3x3 local matrix) landINFLUNZ(lolaA) Calculates the nearness of a cell to land and how exposed to the ocean (waves) it is. getAQUA(lolaA) Interrogates the DAAC OceanColor grids from Goddard SC for Chlorophyll, PAR Sea-Surface Irradiance, Ocean absorbance (chlorophyll based, after Morell work. mapOutpParam(z,nameS,showQ) Makes matplotlib graphic of a parameter (z,nameS) over the area plotLinesXY(x,y,nameS,showQ) plotPointsXY(x,y,nameS,showQ) Plot parameters in various ways. ShowQ tells whether to display the result at run-time. calcMilieu(enviroA,iCr) How suited is the organism (based on it's KB data) to the environment of the gridcell ? At present outputs the maximum suitability, and autrophic potential (leading to growth). Presently based on WD,Temp,Light(PAR). crSuitab,crTrophic are outputs calcTrophics(irrad) Parameterizes carbonate skeletal extension using the benthic irradiation as the driver, following Gattuso and Kleypas. Light is the main source of metabolic energy. heterotrophic uptake is also allowed (eg at night, and at depths). calcLotStocks(creatrA,competnA) Lotke-Volterra calculation of population growth spawning(Creatr) How many spawn are produced ? (Sexual and asexual; at each cell) cohorts(Creatr,cohortsA) Re-calculates the population profile of a organism (at each cell) Some definitions... prodnA - carbonate production kg/m2 insituPetA - insitu materials description detriPetA - detrital materials description thknsA[0] - insitu material thickness (collapsed) thknsA[1] - detrital thickness calcSedProduction(creatrA) Calculate the amount of sediment produced per grid cell based on the population levels present and Calculate the thickness of accumulation per grid cell for the current time step based on sediment production Will also be dependent on sediment breakdown and transport - how much of the produced material is removed, but this can be added later... _myMAIN_ The main driver of the program, including the moving (by loops) over the yx map, tBench computing step and iC organism drivers. seaLEVEL(findKy,step) Varies sea-level and thereby, many other enviro parameters. (Uses the Sedflux 140ky sealevel array for the time being.) Other explanations: 1. Arrays are generally marked with an ending "a" or "A" 2. Strings are generally marked with an ending "s" or "S" 3. Data directly taken from the Knowledge Base is prefixed with "kb" 4. The parameter names of the KB have a web syntax eg creature.param.subparam 5. Usually, counters use "i" or "n", and elements of an array may use "ë" 6. Terms and Arrays are named this way forwhatWhatCoords and then the a/s/e/n/i suffix: eg "crStocksMapA" organism-stock at Yth row-Xth col in array 7. tStart/tStop are the start/stop times based on tNow tNow is the step the model is at, and t can advance ahead of tNow in tBurstLen ''' # ============================================================================ def physIn3D(toMapA,nameS,unitS,showQ): #After: # "http://hyry.dip.jp:8000/scipybook/default/file/05-matplotlib/mplot3d_surface.py" #and "http://matplotlib.sourceforge.net/trunk-docs/examples/mplot3d/surface3d_demo3.html" from mpl_toolkits.mplot3d import axes3d #{1} loudQ=False if loudQ: print ">>physIn3D:",nameS xA, yA = NPY.mgrid[0:modlSZx,0:modlSZy] #{2} zA=NPY.array(elevnA[:]) toMapA=NPY.array(toMapA) #Specify the colours... norm = COL.normalize(toMapA.min(),toMapA.max()) #cc = lambda arg: COL.colorConverter.to_rgba(arg, alpha=0.6) fColorsA=[] for y in range(modlSZy): fColorsA.append([]) for x in range(modlSZx): fColorsA[y].append(MPL.cm.jet(norm(toMapA[y][x]))) #fColorsA.append(norm(toMapA[y][x])) if loudQ: print "fColorsA=",fColorsA fig3d = PLT.figure() ax = axes3d.Axes3D(fig3d) #{3} #v1: ax.plot_surface(xA, yA, zA, rstride=1, cstride=1, cmap = PLT.cm.Blues_r) #{4} surf=ax.plot_surface(xA, yA, zA, rstride=1, cstride=1, facecolors=fColorsA) #{4} ax.set_xlabel("Y") ax.set_ylabel("X") #(Note: CJJ will need to get some flexibility into the graph orientations) ax.set_zlabel("m") #The elevation PLT.title(nameS+":"+runYMDHMs) #The color bar has the colours, not elevation #cb=fig3d.colorbar(surf) cbNorm = COL.Normalize(toMapA.min(),toMapA.max()) cax,kw=MPL.colorbar.make_axes(ax,shrink=.66) cb=MPL.colorbar.ColorbarBase(cax, cmap=MPL.cm.jet,norm=cbNorm) cb.set_label(unitS) fig3d.savefig(dirS+runDirS+"_3dPlots/"+nameS+'_3d.png') if showQ<>0: PLT.show() raw_input("Show") PLT.close() return # ============================================================================ def myMAIN(): global dirS,runDirS,runYMDHMs,lolaA global modlSZx,modlSZy global modlRESN,modlSZa,LoLfHtWdA global tNow,nY,nX,showA,reptCells global pi,d2r global iNUL,fNUL,stockNull global tStep,tRept,tStop,tCompute global tBurstLen,tBurstIni,tBurstEnd,tBurstRun global prjA,prjS global seaLEVL,elevnA #elevnA is adjusted in myMAIN global logPlotQ global prodnMapA,lithoMapA,thknsMapA global colorsA,hexColorsA global showCellQ #Says whether to display on screen for this cell global initClock,startClock,endClock #Program runtimes initClock=time.clock() global logF #Banner .... print "*"*70;print print "\t"*3,"Program carboLOT_*.py"; print print "\t","Chris Jenkins INSTAAR, Peter Burgess RHUL, Donald Potts UCSC" print "\t"*4,"Sep 2010 onwards" print;print "*"*70;print #Initializations ------- pi=NPY.pi; d2r=pi/180. #Use as deg*d2r iNUL=-99; fNUL=-99.; stockNull=10.**-4 seaLEVL=0.0 logPlotQ=True #True= plot stocks as logs; else False=linear loudQ=False #Define the project ==================================================== print; print "Project (location) "+"="*30 prjSelectA=["agin","chag","midw","lhis"] #All lower case prjS="-"; print prjS=raw_input("Please give a 4-char code from "+\ "|".join(prjSelectA)+" or new ...") prjS=prjS.strip().lower() while prjS not in prjSelectA: print prjS=raw_input("Please try again... ") else: print;print "Accepted ! Working on: " runYMDHMs="".join(map(str,time.localtime()[:5])) #The run folder uses project & date/time -------------------- dirS="_"+prjS.lower()+"/" #(Note this needs to be made manually for new project areas, # and have the Setup,Event files and Setup folder installed runDirS=runYMDHMs+"/" try: os.mkdir(dirS+runDirS) except: pass #May already exist... print "dirS+runDirS=",dirS+runDirS #Testing the setup files are present ... stopQ=False try: setupTryF=open(dirS+"_"+prjS+"_SETUP.txt","r") setupTryF.close() except: raw_input("You need to place a SETUP file in to: "+dirS) stopQ=True try: setupTryF=open(dirS+"_"+prjS+"_EVENTS.txt","r") setupTryF.close() except: raw_input("You need to place an EVENTS file in to: "+dirS) stopQ=True if stopQ: sys.exit("Can't continue without the files") #Log file... logF=open(dirS+runDirS+"_"+prjS+"_logFile.txt","w") try: os.mkdir(dirS+runDirS+"_StockDumps") os.mkdir(dirS+runDirS+"_StratDumps") os.mkdir(dirS+runDirS+"_GraphPlots") os.mkdir(dirS+runDirS+"_MapPlots") os.mkdir(dirS+runDirS+"_3dPlots") except: print "Problem making project/run folders" #Get project details --------------------------- outS="\n"+"Project setup variables "+"="*30 valuesA=getSetUp() prjA,lolaA,modlSZa,modlRESN,showA,sectnA,\ tStart,tStop,tStep,tRept,tCompute,tBurstLen,\ sStart=valuesA modlSZx,modlSZy=modlSZa reptCells=2 tNow=tStart-1 #(Begin one step behind actual start) #Location --------------------------------------------------- modlSZa=[modlSZx,modlSZy] outS="\n"+"Project location data "+"="*30+"\n" outS+=" project= "+";".join(map(str,prjA))+"\n" outS+=" lon/lat= "+";".join(map(str,lolaA))+"\n" outS+=" model dims,resoln= "+"".join(map(str,modlSZa))+","\ +str(modlRESN)+"\n" print outS; logF.write(outS) outS="\n"+"Map layout: modlSZa,modlRESN,showA="+\ ";".join(map(str,modlSZa))+","+str(modlRESN)+","\ +";".join(map(str,showA))+"\n" #showA is the cell for which a screen display is given during run LoLfHtWdA=[lolaA[0]-float(modlSZx)*modlRESN/2,\ lolaA[1]-float(modlSZy)*modlRESN/2,\ float(modlSZx)*modlRESN,\ float(modlSZy)*modlRESN] outS+="Map extents: LoLfHtWdA=..."+";".join(map(str,LoLfHtWdA))+"\n" print outS; logF.write(outS) #Biological values -------------------------------------- outS="\n"+"Interrogate organism KnowledgeBase "+"="*20+"\n" print outS; logF.write(outS) foo=getKBdata() #This loads creatrA, proptyA into global #Make discrete colour scale for numOrgnsms ---------------------------- norm = COL.Normalize(0.,float(numOrgnsms-1)) colorsA = [MPL.cm.jet(norm(nF)) for nF in map(float,range(numOrgnsms))] colorsA[0]=(0.5, 0.5, 0.5) #Set (0) to grey... #print "colorsA=",colorsA hexColorsA=[] for rgbA in colorsA: hexColrsS=COL.rgb2hex(rgbA[:3]) hexColorsA.append(hexColrsS) #print rgbA,hexColrsS #print "colorsA,hexColorsA=",colorsA,hexColorsA #Events ------------------------------------------------------ print; print "Setting up the events "+"-"*20 evTimesA=getEvents() #eventsA output is global, evTimesA is local #We use evTimesA to truncate the burst calculations # also potentially to plot lines on graphs #Environmental values ------------------------------------------------------ #foo=groundOverlay([],"","","","ini") #Prepares kml for ground overlays print; print "Set up environmental YX model grids "+"-"*40 kyS=raw_input("Use stored enviro grids (if available; else make new) ? (y/.n*)") print if kyS.lower()=="y": tempA,parA,chlA,irradA,hsA,tpA,dpA,zEnrgyDissiA=getDataStore() #We will assume this is successful... print;print "Note: no maps are drawn from this retrieval" print " but a mask is made" else: #Note: lolaA is the lower left; # the grid is then made from modlSZx;y & modlRESN foo=getELEVN(lolaA) #Get bathy based on a lat/lon #(elevnA is returned via global) tempA=getTEMP(lolaA) #Temps from WOA05 #print "tempA=",tempA parA,chlA=getAQUA(lolaA) #Get light distribution from MODIS AQUA irradA=calcLIGHT(parA,chlA) #(Initial values only, here) hsA,tpA,dpA=getWAVES(lolaA) #Deep water, undirected for the present zEnrgyDissiA=calcENERGY(hsA,tpA,dpA) #(Initial values only, here) #Make the grid dumps ready for the next runs... foo=putSetupStore(tempA,parA,chlA,irradA,hsA,tpA,dpA,zEnrgyDissiA) print " ... enviro grids re-computed" time.sleep(1) #THE MODELLING ========================================================== loudQ=False startClock=time.clock() tBurstIni=-1.; tBurstEnd=-1. #Announce... print;print;print "*"*60 print "Start modelling the carbonate ... " print #Map arrays (very similar to PeterB's) --- prodnMapA=[[0.]*modlSZx for i in xrange(modlSZy)] #4 near-seabed bins #(The 3 upper bins will compact to form the tickness) thknsMapA=[[0.]*modlSZx for i in xrange(modlSZy)] depoMapA=[[0.]*modlSZx for i in xrange(modlSZy)] sedThknsMapA=[[0.]*modlSZx for i in xrange(modlSZy)] #Lithology describes the lithology of the carbonate lithoMapA=[[""]*modlSZx for i in xrange(modlSZy)] detriMapA=[[""]*modlSZx for i in xrange(modlSZy)] #Initialize the main biol data accumulators ----------------- print ">>setIniStock: sStart=",sStart crStocksMapA=[[0.]*modlSZx for i in xrange(modlSZy)] #Will take total matrix of organism stocks # updated by xth-cell then yth-row (at initialization # set in one go) crDeadsMapA=[[0.]*modlSZx for i in xrange(modlSZy)] vacntMapA=[[0.]*modlSZx for i in xrange(modlSZy)] roughnsMapA=[[0.]*modlSZx for i in xrange(modlSZy)] for nY in range(modlSZy): for nX in range(modlSZx): showCellQ=False if [nY,nX]==showA: showCellQ=True #Find the cell's environment wd=seaLEVL-elevnA[nY][nX] #enviroA is the array of environment values for the site; # they are used to calculate habitat suitability enviroA=[['wd',wd],['temp',tempA[nY][nX]],\ ['par',irradA[nY][nX]],['chl',chlA[nY][nX]],\ ['enrgydens',zEnrgyDissiA[nY][nX]] ] #Initialize the stocks as necessary ... #(Will set according to options given in setUp) crStocksMapA[nY][nX]=setIniStock(sStart,nY,nX,enviroA) crDeadsMapA[nY][nX]=[0.]*numOrgnsms vacntMapA[nY][nX]=calcVacancy(crStocksMapA[nY][nX]) roughnsMapA[nY][nX]=calcRoughness(crStocksMapA[nY][nX],\ crDeadsMapA[nY][nX],vacntMapA[nY][nX]) #Note: stocksNowInCellA and deadNowInCellA are reserved # for post-initialization #Deads holds 'stock' that has died this step; it is used for # deposition ( as opposed to accumulation) if showCellQ and loudQ: print "On initialization: "+\ "crStocksMapA[nY][nX],vacntMapA[nY][nX]=",\ crStocksMapA[nY][nX],vacntMapA[nY][nX] foo=putPops(nY,nX,-1,crStocksMapA[nY][nX]) foo=putDead(nY,nX,-1,crDeadsMapA[nY][nX]) print "Initialization for t=-1 is complete" #For the total YX matrix of nCr organism stocks # updated by xth-cell then yth-row (at initialization # set in one go) inBurstQ=False #Initializing while tNow>putPops") #(Uses the stocks for one cell/timestep) #This can be applied as a check in different places #records this time's data for production, etc... stocksNowInCellA=stockSetsInCellA[0] deadNowInCellA=deadSetsInCellA[0] #(Note: bioerosion is rated to the amount of recently # dead, ie defenceless) else: if inBurstQ: #Note: the stockSetsInCellA array is # already prepared and written, so no action # needed on that #Pick up the existing live & dead stocks in cell... stocksNowInCellA,deadNowInCellA=\ getStocksAtT(nY,nX,tNow) #(ie for all cr at YXT) if showCellQ and loudQ: print "Still in burst: tNow,tNow-tBurstIni=",\ tNow,tNow-tBurstIni #print " stocksNowInCellA=",stocksNowInCellA #putPops and putDead not needed else: #Pick up the pre-existing live & dead stocks in cell # and extend them one more step... stocksNowInCellA,deadNowInCellA=\ getStocksAtT(nY,nX,tNow-1) #(ie array of 1-dim cr entries; for previous year) if showCellQ and loudQ: print "Drifting: tNow,=",tNow #print "stocksNowInCellA=",stocksNowInCellA,inBurstQ foo=putPops(nY,nX,tNow,stocksNowInCellA) foo=putDead(nY,nX,tNow,deadSetsInCellA[tGet]) #print "doBurstQ,inBurstQ=",doBurstQ,inBurstQ #Then get appropriate values in time for each cr #(This works for burst and annual array instances) #print stockSetsInCellA.shape,tNow if showCellQ and loudQ: print "tNow,tGet,inBurstQ,stocksNowInCellA=",\ tNow,tGet,inBurstQ,stocksNowInCellA,\ NPY.array(stocksNowInCellA).shape #Which is most abundant organism ? ------------- #(But not if they are at the floor of 10E-4) crDomnt=-99; crDomnc=-99. #In case of no clear winner crDomnc=max(stocksNowInCellA) if crDomnc>stockNull: if stocksNowInCellA.count(crDomnc)==1: crDomnt=stocksNowInCellA.index(crDomnc) crDomntRow.append(crDomnt) crDomncRow.append(crDomnc) if showCellQ and loudQ: print "crDomnt,crDomnc=",crDomnt,crDomnc #Put burst outputs to cellwise accumulators --------- # ready for the next time step... crStocksRow.append(stocksNowInCellA) #Total organism coverage... crStockSumRow.append(sum(stocksNowInCellA)) if showCellQ and loudQ: print "nY,nX,stocksNowInCellA,sum(stocksNowInCellA)=",\ nY,nX,stocksNowInCellA,sum(stocksNowInCellA) foo=show2D(crStockSumMapA,"nYnX crStockSumMapA",False) #Accumulation of carbonate... prodnA=calcSedProduction(stocksNowInCellA) #Living vacncy=calcVacancy(stocksNowInCellA) #Bare roughns=calcRoughness(crStocksMapA[nY][nX],\ deadNowInCellA,vacncy) thknsA,insituPetA,detriPetA=calcSedBuildup(deadNowInCellA) #Dead prodnRow.append(sum(prodnA)) lithoRow.append(insituPetA) detriRow.append(detriPetA) roughnsRow.append(roughns) vacntRow.append(vacncy) thicknRow.append(thknsA[0]) depoRow.append(thknsA[1]) #Note: Not the underground accumulations) elevnA[nY][nX]+=sum(thknsA) # elevation changes with accumulated thickness. # water depth will also change... if elevnA[nY][nX]>1.: print "Broaching: elevnA[nY][nX],tNow,nY,nX,thknsA=",\ elevnA[nY][nX],tNow,nY,nX,thknsA else: #print "On land !" numOnLand+=1 enviroA=[['wd',-1],['temp',-99],['par',-99],['chl',-99]] suitabsInCellA=[0.]*numOrgnsms trophicsInCellA=[[0.]*3 for nO in range(numOrgnsms)] stocksNowInCellA=[stockNull]*numOrgnsms stockSetsInCellA=[[stockNull]*numOrgnsms for tN in range(tNow+1)] deadNowInCellA=[stockNull]*numOrgnsms #Multiple through the burst foo=putPops(nY,nX,tNow,stocksNowInCellA) foo=putDead(nY,nX,tNow,deadNowInCellA) crStocksRow.append([0.]*numOrgnsms) #Total organism coverage... crStockSumRow.append(0.) #No most abundant organism crDomntRow.append(-1) crDomncRow.append(0.) #No additions to carbonate budgets prodnRow.append(0.) lithoRow.append("-") detriRow.append("-") roughnsRow.append(0.) vacntRow.append(0.) thicknRow.append(0.) depoRow.append(0.) #-------------------------------------------------- #End nX loop crDomntMapA[nY]=crDomntRow[:] crDomncMapA[nY]=crDomncRow[:] crStocksMapA[nY]=crStocksRow[:] crStockSumMapA[nY]=crStockSumRow[:] prodnMapA[nY]=prodnRow[:] lithoMapA[nY]=lithoRow[:] detriMapA[nY]=detriRow[:] roughnsMapA[nY]=roughnsRow[:] vacntMapA[nY]=vacntRow[:] thknsMapA[nY]=thicknRow[:] depoMapA[nY]=depoRow[:] #End nY loop if loudQ: foo=show2D(crStockSumMapA,"crStockSumMapA",False) #Redistribute any sediment that exists across the map... depoMapA=distribSED(depoMapA,roughnsMapA) #sedThknsMapA is the whole-of run deposition map for dY in range(modlSZy): sedThknsRow=[0.]*modlSZx for dX in range(modlSZx): sedThknsRow[dX]=sedThknsMapA[dY][dX]+depoMapA[dY][dX] sedThknsMapA[dY]=sedThknsRow[:] #Dump the map characteristics to a file for each time step... foo=dumpCharData("prodn",prodnMapA) foo=dumpCharData("litho",lithoMapA) foo=dumpCharData("detri",detriMapA) foo=dumpCharData("thkns",thknsMapA) foo=dumpCharData("depos",depoMapA) foo=dumpCharData("sedthk",sedThknsMapA) if tNow % tRept==0: #print "\t"+"Making maps" #Create new map mask to recognize new land areas foo=calcMASK() #Dump the start conditions (ie for tStart) foo=reportEnviro(tempA,parA,chlA,irradA,zEnrgyDissiA) #Plot the yx model cells of most abundant organism ---- # and put each grid out to KML... # ...(LoLfHtWdA,projectS,nameS,descnS,imageS) imgNameS=prjS+"_DomStock_t"+str(tNow) #print "crDomntMapA,crDomncMapA",crDomntMapA,crDomncMapA foo=mapDomOrgnsms(crDomntMapA,crDomncMapA,imgNameS,0) '''foo=groundOverlay(LoLfHtWdA,\ "DomntOrgnsms",\ "Dominant organisms of each cell (orgnm #)",\ "_MapPlots/"+imgNameS,"map")''' imgNameS=prjS+"_StockSum_t"+str(tNow) unitS='m^2 inhabited/m^2 area' if loudQ: print "crStockSumMapA=",crStockSumMapA,\ NPY.array(crStockSumMapA).shape foo=mapOutpParam(crStockSumMapA,imgNameS,unitS,0) '''foo=groundOverlay(LoLfHtWdA,\ "TotalStocks",\ "Total Stocks of carbonate organisms (%live areal cover)",\ "_MapPlots/"+imgNameS,"map")''' foo=physIn3D(crStockSumMapA,imgNameS,unitS,0) imgNameS=prjS+"_Prodn_t"+str(tNow) unitS='kg/m^2/yr' #print "prodnMapA=",prodnMapA,NPY.array(prodnMapA).shape foo=mapOutpParam(prodnMapA,imgNameS,unitS,0) '''foo=groundOverlay(LoLfHtWdA,\ "Prodn",\ "Carbonate production (kg/m2)",\ "_MapPlots/"+imgNameS,"map")''' foo=physIn3D(prodnMapA,imgNameS,unitS,0) imgNameS=prjS+"_Thkns_t"+str(tNow) unitS='m/yr' #print "thknsMapA=",thknsMapA,NPY.array(thknsMapA).shape foo=mapOutpParam(thknsMapA,imgNameS,unitS,0) '''foo=groundOverlay(LoLfHtWdA,\ "Thkns",\ "Carbonate thickness (m)",\ "_MapPlots/"+imgNameS,"map")''' foo=physIn3D(thknsMapA,imgNameS,unitS,0) imgNameS=prjS+"_Roughns_t"+str(tNow) unitS='RMS m @ <1m' #print "roughnsMapA=",roughnsMapA,NPY.array(roughnsMapA).shape foo=mapOutpParam(roughnsMapA,imgNameS,unitS,0) '''foo=groundOverlay(LoLfHtWdA,\ "Roughns",\ "Seabed roughness (RMS m @ <1m)",\ "_MapPlots/"+imgNameS,"map")''' foo=physIn3D(roughnsMapA,imgNameS,unitS,0) imgNameS=prjS+"_Depo_t"+str(tNow) unitS='m/yr' #print "depoMapA=",depoMapA,NPY.array(depoMapA).shape foo=mapOutpParam(depoMapA,imgNameS,unitS,0) '''foo=groundOverlay(LoLfHtWdA,\ "Depo",\ "Deposited thickness (m)",\ "_MapPlots/"+imgNameS,"map")''' foo=physIn3D(depoMapA,imgNameS,unitS,0) imgNameS=prjS+"_SedThkns_t"+str(tNow) unitS='m' #print "sedThknsMapA=",sedThknsMapA,NPY.array(sedThknsMapA).shape foo=mapOutpParam(sedThknsMapA,imgNameS,unitS,0) '''foo=groundOverlay(LoLfHtWdA,\ "SedThkns",\ "Sediment Thickness (m)",\ "_MapPlots/"+imgNameS,"map")''' foo=physIn3D(sedThknsMapA,imgNameS,unitS,0) #print "\tEnd of the reporting" if tNow>=tStop: print;print "Reached stoptime for the modelling " print "-"*60 break #Check some things before next tBench cycle... #print "Dims of crStocksMapA|etc:",NPY.array(crStocksMapA).shape print;print "End of the run "+"="*50 endClock=time.clock() #Make a thickness cross-section ----------------------------------- csY,csX=sectnA #From setup currently [4,-1] #(Y-section is used in the development case only; # later allow for X-sectn as well) thknsXSectA=[]; faciesXSectA=[] for csX in range(modlSZx): thknsTa=liftThknsData([csY,csX],"thkns") thknsXSectA.append(thknsTa) #Panel across the map at Y=csY faciesTa=liftDomntData([csY,csX],"live") faciesXSectA.append(faciesTa) #print "thknsXSectA,faciesXSectA=",thknsXSectA,faciesXSectA thknsXSectA=NPY.array(thknsXSectA).transpose() faciesXSectA=NPY.array(faciesXSectA).transpose() elevnSectA=elevnA[csY] #Draw the stratigraphic cross-section graphic ... foo=plotStratig(thknsXSectA,faciesXSectA,elevnSectA,prjS+"_strat_t"+\ str(tStop+1)+"_y"+str(csY),0) #Draw graphical plots of the time series on a section ----------- foo=graphPopulo(sectnA) #print;print ">>groundOverlay: Make the GoogleEarth display..." #foo=groundOverlay([],"","","","end") return # ========================================================================= def show2D(yxGridA,nameS,intQ): #This is just a gadget to display a 2D grid concisely onscreen #Always loud ! def strIntMult(elem): return str(int(elem*factor)) #print "yxGridA=",yxGridA if intQ: factor=1. else: factor=100. print "grid display:",nameS #print [repr(map(strIntMult,yxGA)) for yxGA in yxGridA] print "\n".join([repr(map(strIntMult,yxGA)) for yxGA in yxGridA]) return # ========================================================================= def bilearInterp(yxGridA,zGridA,yxA): #[1,1] is at bottom left #zGridA: [z11,z21,z22,z21] clockwise from lowerleft y,x=yxA y1,x1,y2,x2=yxGridA #LL then UR y,x z11,z21,z22,z12=zGridA Z=z11/((x2-x1)*(y2-y1))*((x2-x)*(y2-y))+\ z21/((x2-x1)*(y2-y1))*((x-x1)*(y2-y))+\ z12/((x2-x1)*(y2-y1))*((x2-x)*(y-y1))+\ z22/((x2-x1)*(y2-y1))*((x-x1)*(y-y1)) print "<>calcRoughness: liveA,deadA,vacance=",\ liveA,deadA,vacance relief1=0.; relief2=0.; area=0. for pCr in range(numOrgnsms): #Live... hght=kbCoverageA[pCr][0]/3. relief1+=liveA[pCr]*hght relief2+=liveA[pCr]*hght**2 #Note: /3 for diam to radius if liveA[pCr]<=0.0001: liveA[pCr]=0. #Effectively zero area+=liveA[pCr] if loudQ: print "Live: hght,relief1,relief2,area=",hght,relief1,relief2,area #Dead... hght=kbSkeletalA[pCr][1] #Collapsed clastsize relief1+=deadA[pCr]*hght relief2+=deadA[pCr]*hght**2 if deadA[pCr]<=0.0001: liveA[pCr]=0. #Effectively zero area+=deadA[pCr] if loudQ: print "Dead: hght,relief1,relief2,area=",hght,relief1,relief2,area #relief also includes vacancy (assume z=0) #relief1+=0.; relief2+=0. area+=vacance if loudQ: print "Final: relief1,relief2,area=",relief1,relief2,area roughns=NPY.sqrt(relief2/area) if loudQ: print "rms roughness in m: roughns=",roughns return roughns # ========================================================================= def calcVacancy(popsA): #Very simple calculation - but complicated issue # especially when ODE is not giving perfect answers... loudQ=False vacancy=1.-sum(popsA) if vacancy<0.: if loudQ: print "Warning: sum of stocks >1.0 ! sum(popsA)=",sum(popsA) vacancy=stockNull return vacancy # ========================================================================= def distribSED(sedsA,roughnsA): #After each t-step it redistributes the sediment... #( will need the sed,bathy,energy and slope maps.) loudQ=False if loudQ: print;print ">>distribSED:" nsewA=[0.,90.,180.,270.] nghbrA=[[+1,0],[0,+1],[-1,0],[0,-1]] #Clckw from cellabove depoA=[[0.]*modlSZx for i in xrange(modlSZy)] #t-step change in sediments #Eventually will ? # Continue the process till all sediment is home, # or for a set number of steps, or for a max distance ? depoSteps=3 #How many interations of the process for dS in range(depoSteps): for sY in range(modlSZy): for sX in range(modlSZx): minSedThk=roughnsA[sY][sX] #Usually keeps <10cm in place dWd=seaLEVL-elevnA[sY][sX] #Kaufman et al Diffusion in Strat models C0=50000.; C1=0.05 #Their standard values m2/yr diffsnC=C0*NPY.exp(-C1*dWd) if sedsA[sY][sX]>=minSedThk/10.: #Have enough sediment ? (tenth of roughness) gradsA=[0.]*4; moveA=[] slope,aspect=getSHORE(sY,sX) aspect=aspect % 360. if loudQ: print "minSedThk,dWd,slope,aspect=",\ minSedThk,dWd,slope,aspect #Account for direction ... #(will send it only through the 4 orthogonal boundaries) for cN in range(4): if abs(aspect-nsewA[cN])<=90.: #Ie downslope only gradsA[cN]=NPY.cos((aspect-nsewA[cN])*d2r) #print "aspect,gradsA=",aspect,gradsA for grad in gradsA: moveA.append(NPY.tan(slope*d2r)*grad*sedsA[sY][sX]) #ie gradient or diffusion can move it out of cell # diffusn send it equally in 4 directions; goes to 0 with depth if loudQ: print "NPY.tan(slope*d2r),diffsnC,grad,moveA=",\ NPY.tan(slope*d2r),diffsnC,grad,slope,moveA '''#Sediment pickup... #After: "http://journals.tdl.org/ICCE/article/viewFile/1263/pdf_79" # and "http://www.ecoh.co.jp/e/tech/pdf/APACE-sediment%283%29.pdf" betas=0.0045E0.64 rhoS=2.5; rhoW=1.0; dGrz=0.5 #mm K=0.8*exp(-2.5*d) if dWdX<0.: Qp=betas/((rhoS-rhoW)*gAccel*eWd)*dWdX else: Qp=0. #dWdX is Dissipation rate=dWaveEnergyFlux/ChangeInWaterLevel''' depoA[sY][sX]+=-sum(moveA) #Home cell #Deliver to neighbours... for nN in range(4): nA=nghbrA[nN] if sY+nA[0]>=0 and sY+nA[0]=0 and sX+nA[1]>setIniStock: sStart,nY,nX",sStart,nY,nX habiSuitA=[stockNull]*numOrgnsms addSome=float(numOrgnsms)/float(numOrgnsms+1) #eg 0.9 for 9 spp #Will allow for some vacant ground if sStart=="suitab": for iCr in range(numOrgnsms): habiSuitA[iCr],foo=calcMilieu(enviroA,iCr) #(Note: same habitat suitability funtion as used through model) #(No need here for trophicsInCellA, so use foo) #(Note: may be zero) #Normalize to <1.0 sumIniStocks=sum(habiSuitA) if loudQ: print "habiSuitA,sumIniStocks=",habiSuitA,sumIniStocks if sumIniStocks>0.: iniStocksA=[hSA/sumIniStocks*addSome for hSA in habiSuitA] else: iniStocksA=[0. for hSA in habiSuitA] else: #(Default) iniStocksA=[1./float(numOrgnsms+1)]*numOrgnsms #(Note: the extra one is for barren seabed) #Apply the frame fraction... for isaCr in range(numOrgnsms): iniStocksA[isaCr]=iniStocksA[isaCr]*kbCoverageA[isaCr][1] #Note: kbCoverageA[isaCr][1] is the living fraction #print "iniStocksA[isaCr],kbCoverageA[isaCr][1]=",\ # iniStocksA[isaCr],kbCoverageA[isaCr][1] #Limit lower bounds to stockNull iniStocksA=[max(iSA,stockNull) for iSA in iniStocksA] if loudQ: print "<>plotStratig:" #Cut the stratal matrices ------------------------------------------ #(stratThknsYXa,stratFaciesYXa come in via global) #faciesProfileA=stratFaciesYXa[:,sectnCut,:] xsectnShapeA=thknsProfileA.shape #[T,x] if loudQ: print "xsectnShapeA=",xsectnShapeA print " thknsProfileA=",thknsProfileA #print " faciesProfileA=",faciesProfileA print " elevnSectA=",elevnSectA #------------------------------------------------------------------- #prepare the plot space... fig2 = PLT.figure() #First draw the bathymetric profile (at page bottom) ---- ax22=PLT.subplot(212) mapSectnA=[xsS+0.5 for xsS in range(xsectnShapeA[1])] #(This moves scale along 0.5 units) #A choice introduced by CJJ 01Jun2011 ... plotBathyTypeS="log" #"log" or "lin" if plotBathyTypeS=="log": mapValueA=[NPY.log10(max(-eSA,0.1)) for eSA in elevnSectA] yPlotLabelS='Log10 WD (m)' elYlim=[3.5,-1.] else: mapValueA=elevnSectA[:] yPlotLabelS='Elevations (m)' elYlim=[10.,-3000.] print "bathy plot: mapValueA=",mapValueA ax22.plot(mapSectnA,mapValueA,'k--') #ax22.grid(True) ax22.set_xlim(-1,xsectnShapeA[1]) ax22.set_ylim(elYlim[0],elYlim[1]) #title('Profile elevations') ax22.set_ylabel(yPlotLabelS) ax22.set_xlabel("Map Section Cell (#)") # Shrink current axis's height by 10% on the bottom... box = ax22.get_position() ax22.set_position([box.x0, box.y0 + box.height * 0.33, box.width, box.height * 0.67]) #Now do the stratig panel -------------------------- ax21=PLT.subplot(211) maxT=0. faciesLeg=[[],[]] legendPatchA=[] tChron=tRept #Will draw time line each tChron time steps (usu =tRept) #chronVerts=[[] for i in range(tStop/tChron)] #Note: chron=0 is the first (non-zero) one #Change the matrix into a structure of the box shapes... for spX in range(xsectnShapeA[1]): #for x-width (ie work on each column) fspX=float(spX) #(Note: cells have origin at their count - like on the maps) tBot=0. tTop=0. #thknsProfileA[-1,spX] #Bottom cell in xth column if loudQ: print "spX,tBot,tTop=",spX,fspX,tBot,tTop print " maxT,nCr,colorsA[nCr]=" #To be continued for spT in range(xsectnShapeA[0]): #for y-height (and start from bottom) fspT=float(spT) if loudQ: print "fspT,fspX",fspT,fspX tBot=tTop #Previous value tTop+=thknsProfileA[spT,spX] maxT=max(maxT,tTop) nCr=faciesProfileA[spT,spX] if tTop==tBot: #No thickness alteration pass else: if loudQ: print " ",spT,spT,tBot,tTop,maxT,\ nCr,colorsA[nCr] #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, ] path = PATH.Path(verts, codes) #The drawing... patch=PATCHES.PathPatch(path, \ facecolor=colorsA[nCr], \ edgecolor=colorsA[nCr], lw=0) ax21.add_patch(patch) #Add to the legend lists if necessary... if not nCr in faciesLeg[0]: legendPatchA.append(patch) #(Note: need to use append because only some of # total organism collection will be plotted.) faciesLeg[0].append(nCr) faciesLeg[1].append(kbOrgnsmA[nCr]) #Every so often add a white time line... #(Here just compile the lines...) '''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=chronVerts for chronV in chronVerts: print "chronV=",chronV PLT.plot(chronV[:][0], chronV[:][1], lw=1, color='w')''' ax21.set_ylabel("Accumulation (m)") #Following the style of ... # "http://stackoverflow.com/questions/4700614/how-to-put-the-legend-out-of-the-plot" ax21.set_xlim(-1,xsectnShapeA[1]) ax21.set_ylim(0,maxT*1.5) # Shrink current axis's height by 10% on the bottom... box = ax21.get_position() ax21.set_position([box.x0, box.y0 + box.height * 0.33, box.width, box.height * 0.67]) #Legend ------------------------- leg=ax21.legend(legendPatchA,faciesLeg[1],loc='upper center',labelsep=0,\ bbox_to_anchor=(0.5, -0.2),ncol=3) for t in leg.get_texts(): t.set_fontsize(6) # the legend text fontsize''' print "nameS,runYMDHMs=",nameS,runYMDHMs PLT.title(nameS+"_"+runYMDHMs) fig2.savefig(dirS+runDirS+"_GraphPlots/"+nameS+'.png') if drawQ==1: fig2.show() if loudQ: print "<>liftThknsData: gYXa,charS,tStop=",gYXa,charS,tStop thknsTa=[] #Takes outputs for tThkns in range(tStop): #print "Reading:",dirS+runDirS+"_StratDumps/_"+\ # charS+"_"+str(tThkns)+".txt" thknsInF=open(dirS+runDirS+"_StratDumps/_"+\ charS+"_"+str(tThkns)+".txt","r") doneQ=False gY=-1 #Row counter #Header... thknsInS=thknsInF.readline() while thknsInF: thknsInS=thknsInF.readline() if thknsInS=="": break gY+=1 #print "thknsInS=",thknsInS if gY==gYXa[0]: #Decode the line and get the data... thknsInA=thknsInS.strip().strip(";").split(",") thknsVal=thknsInA[gYXa[1]] #print "thknsInA,tThkns,gY=",\ # thknsInA,tThkns,gY,thknsVal thknsTa.append(float(thknsVal)) break thknsInF.close() if loudQ: print "<>liftDomntData: charS,gYXs=",charS,gYXs #Get the dominant species ------------------------------------- faciesTa=[0]*(tStop+1) #Takes the outputs #The files here are per grid-cell... #print "Reading:",dirS+runDirS+"_StockDumps/_"+\ # charS+"_"+str(tFacies)+".txt" try: faciesInF=open(dirS+runDirS+"_StockDumps/_live_"+gYXs+".txt","r") except: #Must be on land (return Null) ... return faciesTa while faciesInF: #(No header in this case) faciesInS=faciesInF.readline() if faciesInS=="": break #print "faciesInS=",faciesInS #Decode the line and get the data... faciesInA=faciesInS.strip().split(";") tFacies=int(faciesInA[1].replace("T=","")) faciesA=faciesInA[2].replace("Stocks=","").split(",") if tFacies<=tStop: #print "tFacies,faciesA=",tFacies,faciesA #(faciesA is the list of organism stocks at time tFacies; # note that because of Burst, tFacies may occur several # times in the file; they will overwrite to keep the last) faciesTa[tFacies]=NPY.argmax(map(float,faciesA)) #(ie Index of dominant organism - or first on if a tie) #(Note that because of Burst the file may contain tFacies>tStop) faciesInF.close() #print "<>getEvents:" #(Note: This module lifted from CarboCELL 20Apr2011 and modified) eventsA=[] #7 elements evTimesA=[] '''Example 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(dirS+"_"+prjS+"_EVENTS.txt","r") eventsS=eventsF.readline() #Header evOutS=eventsS #Prepares for dumping to log file while eventsF: eventsS=eventsF.readline() if eventsS=="": break #EOF evOutS+=eventsS #For the log file output eventsS=eventsS.strip() #print "eventsS=",eventsS #Deal with commented lines... if eventsS[0]=="#": pass else: evA=eventsS.split("\t") #print "evA=",evA try: mortTargetsA=map(int,evA[3]\ .strip("[").strip("]").split(",")) except: mortTargetsA=[] try: immiTargetsA=map(int,evA[5]\ .strip("[").strip("]").split(",")) except: immiTargetsA=[] eventsA.append([evA[0],int(evA[1]),\ float(evA[2]),mortTargetsA,\ float(evA[4]),immiTargetsA]) nReptd=1 while int(evA[1])*nReptd>putPops: tPut=",tPut yxAddrS=str(nY)+"_"+str(nX) popsOutF=open(dirS+runDirS+"_StockDumps/_live_"+yxAddrS+".txt","a") popsOutF.write("YX="+yxAddrS+"; T="+str(tPut)+\ "; Stocks="+",".join(map(str,stocksNowInCellA))+"\n") popsOutF.close() return # ============================================================================ def getStocksAtT(pY,pX,tGet): tStockA=[]; tLiveA=[]; tDeadA=[] #Prepare for append loudQ=False if [pY,pX]==showA and loudQ: print ">>getStocksAtT:, pY,pX,tGet",pY,pX,tGet #Obtain the pop stocks data for time t for the cell from # the dumped file ... yxAddrS=str(pY)+"_"+str(pX) recrdsA=["live","dead"] #Address these two different records (file collections) for recrdsS in recrdsA: try: stockGetF=open(dirS+runDirS+"_StockDumps/_"+recrdsS+"_"+yxAddrS+".txt","r") except: #On land ? No graph possible. if loudQ: print "Cell on land" tLiveA=[[stockNull]*numOrgnsms]*(tStop+1) if [pY,pX]==showA and loudQ: print "len(allPopuloA),tStop=",len(allPopuloA),tStop return allPopuloA #Null result #Can go on to make a graph... dumpTs="." #An array of the times in the dump as they are read nStock=0 while stockGetF: stockGetS=stockGetF.readline() if stockGetS=="": break #EOF stockGetS=stockGetS.strip() #print "stockGetS=",stockGetS if stockGetS=="": pass #Empty else: stockGetA=stockGetS.replace(";","|").replace("=","|").split("|") yxRead,tRead=[stockGetA[1],float(stockGetA[3])] #YX,T #print "tRead,stockGetA=",tRead,stockGetA if tRead==tGet: nStock+=1 tStockA=map(float,stockGetA[-1].split(",")) #This will keep only the last one of tGet if [pY,pX]==showA and loudQ: print "tLiveA,nPop=",tStockA,nStock stockGetF.close() if recrdsS=="live": tLiveA=tStockA[:] elif recrdsS=="dead": tDeadA=tStockA[:] if [pY,pX]==showA and loudQ: print "<>putDead: tPut=",tPut yxAddrS=str(nY)+"_"+str(nX) deadOutF=open(dirS+runDirS+"_StockDumps/_dead_"+yxAddrS+".txt","a") deadOutF.write("YX="+yxAddrS+"; T="+str(tPut)+\ "; Dead="+",".join(map(str,deadNowInCellA))+"\n") deadOutF.close() return # ============================================================================ def getSetUp(): #Read the author-edited setup file ... # and assign the values to variables in the program setUpF=open(dirS+"_"+prjS+"_SETUP.txt","r") setUpS=setUpF.read() #All in one gulp setUpS=setUpS.strip()+"\n" setUpF.close() #Dump this to the run log file ... logF.write("\n\n"+""*60+"\n") logF.write(setUpS) logF.write("\n"+""*60+"\n") #Echo the file to the run folder... echoF=open(dirS+runDirS+"_"+prjS+"_SETUP.echo","w") echoF.write(setUpS) echoF.close() setUpDataA=[] paramsA=['prjA','lolaA','modlSZa','modlRESN','showA','sectnA',\ 'tStart','tStop','tStep','tRept','tCompute','tBurstLen',\ 'sStart'] formatsA=['textA','floatA','intA','float','intA','intA',\ 'int','int','int','int','int','int',\ 'text'] valuesA=[] setUpLinesA=setUpS.split("\n") for setUpLine in setUpLinesA: if setUpLine.strip()<>"": setUpDataS=setUpLine.split("\t")[0] if setUpDataS.strip()[0]<>"#": setUpDataS=setUpDataS.strip()\ .replace("[","").replace("]","")\ .replace('"',"") setUpValuA=setUpDataS.split("=") print "setUpValuA=|",setUpValuA,"|",setUpValuA[:][0] try: datIdx=paramsA.index(setUpValuA[0]) except: print "Revise the startup - parameter ",setUpValuA[0] raw_input (" not expected !") if formatsA[datIdx]=="float": setUpValuA[1]=float(setUpValuA[1]) elif formatsA[datIdx]=="int": setUpValuA[1]=int(setUpValuA[1]) elif formatsA[datIdx]=="intA": setUpValuA[1]=map(int,setUpValuA[1].split(",")) elif formatsA[datIdx]=="floatA": setUpValuA[1]=map(float,setUpValuA[1].split(",")) elif formatsA[datIdx]=="textA": setUpValuA[1]=setUpValuA[1].split(",") elif formatsA[datIdx]=="text": pass #Unchanged valuesA.append(setUpValuA[1]) #print "datIdx,paramsA[],formatsA[],valuesA=",\ # datIdx,paramsA[datIdx],formatsA[datIdx],valuesA[datIdx] print "paramsA,formatsA,valuesA=",\ paramsA,formatsA,valuesA print "<>getPopsSeries:" #Obtain the pop stocks data for the cell from # the dumped file in the form of a time series for each organism ... #Note: But the data is held in the files as time slices (sets) yxAddrS=str(pY)+"_"+str(pX) try: popDumpF=open(dirS+runDirS+"_StockDumps/_live_"+yxAddrS+".txt","r") except: #On land ? No graph possible. if loudQ: print "Cell on land" allPopuloA=[[stockNull]*numOrgnsms]*(tStop+1) if loudQ: print "len(allPopuloA),tStop=",len(allPopuloA),tStop return allPopuloA #Null result #Can go on to make a graph... dumpTs="." #An array of the times in the dump as they are read nPop=0 while popDumpF: popDumpS=popDumpF.readline() if popDumpS=="": break #EOF popDumpS=popDumpS.strip() #print "popDumpS=",popDumpS if popDumpS=="": pass #Empty else: popDumpA=popDumpS.replace(";","|").replace("=","|").split("|") dumpYX,dumpT=[popDumpA[1],float(popDumpA[3])] #YX,T #print "dumpT,popDumpA=",dumpT,popDumpA if ( not "."+str(int(dumpT))+"." in dumpTs ): if dumpT<=int(tStop): nPop+=1 populoA=map(float,popDumpA[-1].split(",")) if [pY,pX]==showA and loudQ: print "populoA,nPop=",populoA,nPop allPopuloA.append(populoA) #The stocks dumpTs+=str(int(dumpT))+"." popDumpF.close() if [pY,pX]==showA and loudQ: print "<>graphPopulo: " if loudQ: print " reptCells,reportCellsA=",reptCells,reportCellsA csY,csX=sectnA #From setup for pY in [csY]: for pX in range(modlSZx): allPopuloA=getPopsSeries(pY,pX) allPopuloA=NPY.array(allPopuloA).transpose() #print "allPopuloA(transposed)=",\ # NPY.array(allPopuloA).shape,allPopuloA # ---------------------------------------------- wd=seaLEVL-elevnA[pY][pX] print "Graphing cell: pY,pX,wd,aPA.shape {eg (n,n+1)}=",\ pY,pX,wd,NPY.array(allPopuloA).shape #(Note: Originally this wd test was on final wd, so # if the cell grew above SL, not plotted) pltTitleS=prjS+"_stocks_"+str(pY)+"yx"+str(pX)+\ "_"+str(int(wd))+"wd" #print "pltTitleS=",pltTitleS #Display in log form ? (CJJ Oct 2010) if logPlotQ: for cB in range(len(allPopuloA)): allPopuloA[cB]=map(NPY.log10,NPY.array(allPopuloA[cB])) #print "allPopuloA[cB]=",allPopuloA[cB] #Otherwise unchanged ... if showCellQ and loudQ: print "All log10 stocks: NPY.array(allPopuloA).shape=",\ NPY.array(allPopuloA).shape #print "\n".join(map(str,allPopuloA)) yrsA=range(tStop+1.) #+1 for the ending year foo=graphStockSeqnce(yrsA,allPopuloA,pltTitleS,0) #else: on land return #======================================================================= def getTEMP(lolaA): loudQ=False if loudQ: print;print ">>getTEMP: lolaA,modlRESN=",lolaA,modlRESN #(elevnA comes in through global) #Will report the seasonal range of the 1deg square - at this wd packTYPs="-1: print "Sought data at: seekStep,seekLoLaA=",\ seekStep,seekLoLaA''' #Now get the real values ------------------------------------ seasnlA=[] #holds the 4 seasonal values for zN in zNa: cellPOSa=[-88]*4 for iSSN in range(4): cellPOSa[iSSN]=iSSN*(180*360*33)+zN*(180*360)+grdPOS foo=tempGRDf.seek((cellPOSa[iSSN])*recordSZE) #(... Finds the one before) gridVAL=tempGRDf.read(recordSZE) #Reads actual (dataVAL,)=STR.unpack(packTYPs,gridVAL) seasnlA.append(int(dataVAL)) #print "cellPOSa,seasnlA,wd=",cellPOSa,seasnlA,wd # ----------------------------------------------- while -99 in seasnlA: seasnlA.remove(-99) if len(seasnlA)==0: seasnlA=[-99] tmAx.append(NPY.median(seasnlA)) ''' if NPY.median(seasnlA)<-10: print "nY,nX,cellPOSa,seasnlA,wd,zNa=",\ nY,nX,cellPOSa,seasnlA,wd,zNa print " seekLoLaA,grdPOS=",seekLoLaA,grdPOS #raw_input("!")''' tempA.append(tmAx) #print "tempA,lolaA=","\n".join(map(str,tempA)),lolaA tempGRDf.close() #Also make a plot in matplotlib... unitS="deg C" foo=mapOutpParam(tempA,"tempA"+"_input",unitS,0) if loudQ: print "<>getELEVN:" #definition and dimensions of the grid -------------------- print " GEBCO 1-minute Bathymetry: ",\ btyDIRs+"GEBCO_1mn/"+"GridOne.grd" fINgrd=open(btyDIRs+"GEBCO_1mn/"+"GridOne.grd","rb") packTYPs=">getHyperElevn ..................................." if loudQ: print;print "lolaSEEKa=",lolaSEEKa #Bilinear (ie hyperbolic) interpolation '''fBiLin=open(prjS+"_haveAlook.txt","a") nP=0''' grdSZEa=[21601, 10801] grdRESa=[1./60., 1./60.] grdBOXa=[-180., -90., 180., 90.] #print "grdSZEa,grdRESa,grdBOXa=",grdSZEa,grdRESa,grdBOXa dataVAL=NPY.array([0.]*4) #print "1_min_GEBCO(point): grdPOS=",grdPOS #print btyDIRs+"GEBCO_1mn/"+"GridOne.grd" fINgrd=open(btyDIRs+"GEBCO_1mn/"+"GridOne.grd","rb") packTYPs="i"; recordSZE=STR.calcsize(packTYPs) #print "unpack size: packTYPs,",packTYPs,STR.calcsize(packTYPs) seekPOSa=[int((lolaSEEKa[0]-grdBOXa[0])/grdRESa[0]),\ grdSZEa[1]-int((lolaSEEKa[1]-grdBOXa[1])/grdRESa[1])] #X,Y (Y is top down !) Note: DO want in not nint #print "lolaSEEKa,seekPOSa=",lolaSEEKa,seekPOSa #Find the discrepancy for first cell (LL of 4) dX,dY=[lolaSEEKa[0]-(seekPOSa[0]*grdRESa[0]+grdBOXa[0]),\ lolaSEEKa[1]-((grdSZEa[1]-seekPOSa[1])*grdRESa[1]+grdBOXa[1])] if dX>grdRESa[0] or dY>grdRESa[1]: print "binomial seprns > cellsize ! dX,dY=",dX,dY raw_input("Stopped !") #Get the 4point up-Z around point... dPOSa=[[0,0],[+1,0],[0,+1],[+1,+1]] #[X,Y] #See: "http://www.geocomputation.org/1999/082/gc_082.htm" dPi=-1 for dPOS in dPOSa: dPi+=1 seekA=[seekPOSa[0]+dPOS[0],seekPOSa[1]-dPOS[1]] #Note: - because Y grid pos is from top down) if (seekA[1]<0 or seekA[1]>grdSZEa[1]) or \ (seekA[0]<0 or seekA[0]>grdSZEa[0]): print "getASSOCgrid: Out of grid" dataVAL[dPi]=iNUL else: #OK... grdPOS=seekA[1]*grdSZEa[0]+seekA[0] #X,Y #print "Point in grid: dPOS,grdPOS=",dPOS,grdPOS foo=fINgrd.seek(grdPOS*recordSZE) #Finds one before gridVAL=fINgrd.read(recordSZE) #Reads actual dataVALa=STR.unpack(packTYPs,gridVAL) dataVAL[dPi]=dataVALa[0] #print "dataVALa,dataVAL=",dataVALa,dataVAL if loudQ: print "4 Points LL>>aclckw: seekA,dataVAL[dPi]=",\ seekA,dataVAL[dPi] #Where does the Station lie in the 4 nodes surrounding ? '''dX,dY=[(lolaSEEKa[0]-grdBOXa[0])/grdRESa[0] % 1.,\ (lolaSEEKa[1]-grdBOXa[1])/grdRESa[1] % 1.]''' #ie: Cell offset=real modulus of 1. in cellsize units #(These are the ratio distances from the bottom left hand corner) #Do the Bilinear=Hyperbolic interpolation: #See: "http://www.geocomputation.org/1999/082/gc_082.htm" a00=dataVAL[0] a10=(dataVAL[1]-dataVAL[0]) a01=(dataVAL[2]-dataVAL[0]) a11=(dataVAL[0]-dataVAL[1]-dataVAL[2]+dataVAL[3]) if loudQ: print "dataVAL,a00,a10,a01,a11=",dataVAL,a00,a10,a01,a11 hX=dX/grdRESa[0] hY=dY/grdRESa[1] rHYPO=a00+a10*hX+a01*hY+a11*hX*hY '''#Monitor... for dPi in range(4): nP+=1 dPOS=dPOSa[dPi] seekA=[seekPOSa[0]+dPOS[0],seekPOSa[1]+dPOS[1]] fBiLin.write(str(nP)+","+str(dPOS[0])+","+str(dPOS[1])+","+\ str(seekA[0])+","+str(seekA[1])+","+\ str(dataVAL[dPi])+","+\ str(lolaSEEKa[0]+grdRESa[0]*dPOS[0])+","+\ str(lolaSEEKa[1]+grdRESa[1]*dPOS[1])+\ "\n") fBiLin.close()''' #(This was tested by CJJ 16Aug04 in ‘tst_bathy_hyper.xls’ # and verified ) if loudQ: print "<>get3x3elevn ..................................." #print "lolaSEEKa=",lolaSEEKa grdSZEa=[21601, 10801] grdRESa=[1./60., 1./60.] grdBOXa=[-180., -90., 180., 90.] #print "grdSZEa,grdRESa,grdBOXa=",grdSZEa,grdRESa,grdBOXa seekPOSa=[int((float(lolaSEEKa[0])-grdBOXa[0])/grdRESa[0]-0.5)+1,\ grdSZEa[1]-int((float(lolaSEEKa[1])-grdBOXa[1])/grdRESa[1]+0.5)-1] #print "Point: lolaSEEKa,seekPOSa=",lolaSEEKa,seekPOSa #*POSa are X,Y #Get the 3x3 block around point... dPOSa=[[-1,+1],[+0,+1],[+1,+1],\ [-1,+0],[+0,+0],[+1,+0],\ [-1,-1],[+0,-1],[+1,-1]] dataVAL=NPY.array([0]*9) #print "1_min_GEBCO(point): grdPOS=",grdPOS #print btyDIRs+"GEBCO_1mn/"+"GridOne.grd" fINgrd=open(btyDIRs+"GEBCO_1mn/"+"GridOne.grd","rb") packTYPs="i"; recordSZE=STR.calcsize(packTYPs) #print "unpack size: packTYPs,",packTYPs,STR.calcsize(packTYPs) dPi=-1 for dPOS in dPOSa: #Working through the 3x3 local matrix dPi+=1 gdPOSa=[seekPOSa[0]+dPOS[0],seekPOSa[1]+dPOS[1]] #print "Point: dPi,gdPOSa=",dPi,gdPOSa if (gdPOSa[1]<0 or gdPOSa[1]>grdSZEa[1]) or \ (gdPOSa[0]<0 or gdPOSa[0]>grdSZEa[0]): print "getASSOCgrid: Out of grid" dataVAL[dPi]=iNUL else: #OK... grdPOS=gdPOSa[1]*grdSZEa[0]+gdPOSa[0] #print "Point in grid: dPOS,grdPOS=",dPOS,grdPOS foo=fINgrd.seek(grdPOS*recordSZE) #Finds one before gridVAL=fINgrd.read(recordSZE) #Reads actual dataVALa=STR.unpack(packTYPs,gridVAL) dataVAL[dPi]=dataVALa[0] #print "dataVALa,dataVAL=",dataVALa,dataVAL dataVAL.resize(3,3) #print "dataVAL=",dataVAL fINgrd.close() #Compute some statistics... dV=NPY.array(dataVAL) #print "dV=",dV #Water depth (-ve grid elevation): gdELVN=float(dV[0,0]) #Usual slope index (degrees): gdSLP1=float((dV[0,0]+dV[0,1]*2.+dV[0,2])-(dV[2,0]+dV[2,1]*2.+dV[2,2]))\ /(8.*grdSZEa[0]*(60.*1832.)) #Y gdSLP2=float((dV[0,0]+dV[0,1]*2.+dV[0,2])-(dV[2,0]+dV[2,1]*2.+dV[2,2]))\ /(8.*grdSZEa[0]*(60.*1832.*NPY.cos(float(lolaSEEKa[1])*d2r))) #X gdSLP=NPY.sqrt(gdSLP1**2+gdSLP2**2) #Radian gradient (corrected for cellsize, m heights) #Dartnell roughness (metres): gdRGH=max(dV.ravel())-min(dV.ravel()) #Terrain Roughness Index (metres) (Valentine et al.,2004) gdTRI1=abs(dV[0,2]-dV[1,1])+abs(dV[1,2]-dV[1,1])+\ abs(dV[2,2]-dV[1,1])+abs(dV[0,1]-dV[1,1]) gdTRI2=abs(dV[2,1]-dV[1,1])+abs(dV[2,0]-dV[1,1])+\ abs(dV[1,0]-dV[1,1])+abs(dV[0,0]-dV[1,1]) #|z(-1,1) - z(0,0)| + |z(0,1) - z(0,0)| + #|z(1,1) - z(0,0)| + |z(-1,0) - z(0,0)| + #|z(1,0) - z(0,0)| + |z(1,-1) - z(0,0)| + #|z(0,-1) - z(0,0)| + |z(-1,-1) - z(0,0)| gdTRI=(gdTRI1+gdTRI2)/8. #Rugosity (Lundblad, 2004; Jenness, 2003) ---------- #circumCLa=[[-1,+1],[+0,+1],[+1,+1],\ # [-1,+0],[+1,+0],\ # [-1,-1],[+0,-1],[+1,-1]] #for circN in range(len(circumCLa)): # hypot=NPY.sqrt(dV[circumCLa(circN)]**2+dV[0,0]**2) #--------------------------------------------------- #print "gdELVN,gdSLP,NPY.arctan(gdSLP),gdRGH,gdTRI=",\ # gdELVN,gdSLP,NPY.arctan(gdSLP),gdRGH,gdTRI gdSLPpct=gdSLP*100. return gdELVN,gdSLPpct,gdRGH #Slope is pcnt gradient #======================================================================= def getWAVES(lolaA): loudQ=False #Global Wave climate --------------------------------------------- #print "Wavewatch III Global Wave Climate Grid" #Note: Hs is input as integer values from _Enviro/wwIII up to 10 m wavesDIRs="../../dbSEABED/_db9/_Enviro/wwIII/" print;print ">>getWAVES (Wavewatch NWW3):" pNUL=-999. grdSZEa=[1440, 628] #[nCOLS x,nROWS y] grdRESa=[1./4., 1./4.] grdBOXa=[-0.125,-78.125, pNUL, pNUL] #llX,llY,urX,urY #(Last 2 unneeded) (Note: starts from Greenwich !) #print "grdSZEa,grdRESa,grdBOXa=",grdSZEa,grdRESa,grdBOXa datTYPa=["hsNWW3","tpNWW3","dpNWW3"] unitsA=["m","sec","deg"] #Height,period,direction dt=-1 for datTYPs in datTYPa: wavesGridNameS=wavesDIRs+datTYPs+".grd" print "wavesGridNameS=",wavesGridNameS #print "hs WWIII: grdPOS,seekPOSa=",grdPOS, seekPOSa wavesInF=open(wavesGridNameS,"rb") packTYPs="i"; recordSZE=STR.calcsize(packTYPs) dt+=1 unitS=unitsA[dt] dPi=0; wavesA=[] #Initializations #(Move over the modlSZx*modlSZy grid) for nY in range(modlSZy): wavesYa=[] #Initialize the y-th row for nX in range(modlSZy): #At 100x100 this equates to PeterB's grid seekLoLaA=[lolaA[0]+modlRESN*float(nX-(modlSZx/2)),\ lolaA[1]+modlRESN*float(nY-(modlSZy/2))] #print "Point raw: seekLoLaA=",seekLoLaA seekLoLaA=[((seekLoLaA[0]+360.) % 360.),\ ((seekLoLaA[1]+90) % 180.)-90.] # (rounds to 180deg EW, 90deg NS) #print "Point rounded: seekLoLaA=",seekLoLaA seekPOSa=[\ int((float(seekLoLaA[0])-grdBOXa[0])/grdRESa[0]-0.5)+1,\ grdSZEa[1]-int((float(seekLoLaA[1])-grdBOXa[1])/grdRESa[1]+0.5)-1\ ] #print "Point: seekPOSa=",seekPOSa #*POSa are X,Y if (seekPOSa[1]<0 or seekPOSa[1]>grdSZEa[1]) or \ (seekPOSa[0]<0 or seekPOSa[0]>grdSZEa[0]): print "Location out of **WWIII Grid: seekPOSa,seekLoLaA=",\ seekPOSa,seekLoLaA #siteCHARe[1] unchanged else: #OK... #print "Point in grid" grdPOS=seekPOSa[1]*grdSZEa[0]+seekPOSa[0]+1 foo=wavesInF.seek(grdPOS*recordSZE) #Finds one before gridVAL=wavesInF.read(recordSZE) #Reads actual wavesValA=STR.unpack(packTYPs,gridVAL) wavesVal=wavesValA[0] #Convert to shear stress at seabed based on local bathymetry # but regional (0.25dg) wave climate if wavesVal<0: wavesVal=float('nan') #Set nulls to NaN wavesYa.append(wavesVal) #(wavesYa is the y-collection at the x-level) #print "wavesYa=",wavesYa wavesA.append(wavesYa) wavesInF.close() #Also make a plot in matplotlib... foo=mapOutpParam(wavesA,datTYPs+"_input",unitS,0) #lc units ! #Note: '_input' distinguishes these from the stocks,etc maps #Assemble the output grids ... #print "datTYPs=",datTYPs if datTYPs=="hsNWW3": hsA=wavesA[:] if datTYPs=="tpNWW3": tpA=wavesA[:] if datTYPs=="dpNWW3": dpA=wavesA[:] #Note: these just give the basis for calculating energies for each cell # and for each time step, as the topography changes #print "<>calcENERGY:" waveTypeA=["deep","shallow","transitional"] wRho=1024.; gAccel=9.81 #kg/m3, m/s2 zEnrgyDissiA=[]; slopeA=[]; aspectA=[] for eY in range(modlSZy): zEnrgyDissiAx=[]; slopeAx=[]; aspectAx=[] for eX in range(modlSZx): eWd=seaLEVL-elevnA[eY][eX] #Intercept the on-land cells... if eWd<-1.0: zEnrgyDissiAx.append(float('nan')) slopeAx.append(float('nan')) aspectAx.append(float('nan')) else: #Submergent... #Calc Energy first --------------------------------- wvEnrgyDens=wRho*gAccel*(hsA[eY][eX]**2)/2. #Actually energy density (E/m2) wvOmeg=2.*pi/tpA[eY][eX] waveVa=[gAccel/wvOmeg,-99.,NPY.sqrt(gAccel*eWd)] #Ie: deep,transitional,shallow #Linear interpolation for transitional cases - interim CJJ waveVa[1]=(waveVa[0]+waveVa[2])/2. waveLintrim=waveVa[0]*tpA[eY][eX] if waveLintrim/2.eWd: waveType=2 #Shallow else: waveType=1 #Transitional wvLen=waveVa[waveType]*tpA[eY][eX] wvVel=waveVa[waveType] wvEnrgySrfc=wRho*gAccel*(hsA[eY][eX])**2/16. if loudQ: print print "eWd,wvEnrgySrfc,wvOmeg,wvVel,wvLen,hsA[eY][eX]=",\ eWd,wvEnrgySrfc,wvOmeg,wvVel,wvLen,hsA[eY][eX] #Then Dissipation ---------------------------------- #(After: CERC - Nearshore Wave Breaking and Decay # Technical Report CERC-93-11 July 1993 # by Jane M. Smith Coastal Engineering Research Center # "http://www.dtic.mil/cgi-bin/GetTRDoc?Location=U2&doc=GetTRDoc.pdf&AD=ADA268810") wvCf=0.15 #Friction coefficient (value approx 0.05-0.25) Cgrp=wvVel/2. Kappa=0.15 #Empirical decay constant CERC wvRmsHt=hsA[eY][eX]/2.82 #Offshore value (frmla after NOAA) wvPer=tpA[eY][eX] wvNumKp=2.*pi/wvLen if loudQ: print "Cgrp,wvRmsHt,wvPer,wvNumKp=",Cgrp,wvRmsHt,wvPer,wvNumKp slope,aspect=getSHORE(eY,eX) aspect=aspect % 360. #Angle wave making to shore; zero if wave going downslope wvTheta=abs((dpA[eY][eX]-180.)-aspect) if wvTheta>90.: wvTheta=0. #Waves are going offshore if loudQ: print "dpA[eY][eX],wvTheta,aspect=",dpA[eY][eX],wvTheta,aspect slopeAx.append(slope) aspectAx.append(aspect) if wvTheta==0.: #No slope effect (just bottom friction)... d_ECGcT_dX=(wRho*wvCf)/(6.*pi)*((2.*pi*wvRmsHt)\ /(wvPer*NPY.sinh(wvNumKp*eWd)))**3 else: #Bottom slope is involved... ECg=wvEnrgySrfc*Cgrp #Energy flux assoc with stable wv height #ECgS is the stable energy flux associated with the stable wave height Gam=0.4 #Approximate value of Gam Hstabl=Gam*eWd wvEnrgyStabl=wRho*gAccel*(Hstabl**2)/16. ECgS=wvEnrgyStabl*Cgrp #Flux; CgrpS=Cgrp is CJJ approximation if loudQ: print "wvEnrgySrfc,ECg,wvEnrgyStabl,ECgS=",\ wvEnrgySrfc,ECg,wvEnrgyStabl,ECgS d_ECGcT_dX=Kappa/eWd*(ECg-ECgS)+(wRho*wvCf)/(6.*pi)*\ ((2.*pi*wvRmsHt)/(wvPer*NPY.sinh(wvNumKp*eWd*d2r)))**3 #where ECGcT=ECg*NPY.cos(wvTheta*d2r), so... d_ECg_dX=d_ECGcT_dX/NPY.cos(wvTheta*d2r) if loudQ: print "NPY.sinh(wvNumKp*eWd*d2r),wvNumKp*eWd*d2r=",\ NPY.sinh(wvNumKp*eWd*d2r),wvNumKp*eWd*d2r wvEnrgyDissi=d_ECg_dX #Dissipation rate in J/m2/m #(change in energy flux/m) if loudQ: print "d_ECg_dX,d_ECGcT_dX,wvEnrgyDissi=",\ d_ECg_dX,d_ECGcT_dX,wvEnrgyDissi #What is the energy at depth z in the water ? zEnrgyDissi=wvEnrgyDissi*NPY.exp(-2.*pi*eWd/wvLen) if loudQ: print "waveTypeA,zEnrgyDissi",\ waveTypeA[waveType],zEnrgyDissi #Seem to be some faults in CERC method... if abs(zEnrgyDissi)<1.: zEnrgyDissi=1. #Joules (not much) if zEnrgyDissi<0.: if loudQ: print "-ve Dissi ! zEnrgyDissi=",\ zEnrgyDissi,"(set to 0.)" zEnrgyDissi=0. zEnrgyDissiAx.append(zEnrgyDissi) zEnrgyDissiA.append(zEnrgyDissiAx) slopeA.append(slopeAx) aspectA.append(aspectAx) if loudQ: print "zEnrgyDissiA=",zEnrgyDissiA,len(zEnrgyDissiA),\ len(zEnrgyDissiA[0]) #print "slopeA=",slopeA,len(slopeA),len(slopeA[0]) #print "aspectA=",aspectA,len(aspectA),len(aspectA[0]) #Also make a plot in matplotlib... #Note: '_input' distinguishes these from the stocks,etc maps unitS="Joules/m2/m" foo=mapOutpParam(zEnrgyDissiA,"WavDissip"+"_input",unitS,0) #lc units ! foo=physIn3D(zEnrgyDissiA,"WavDissip"+"_input",unitS,0) unitS="degrees" foo=mapOutpParam(slopeA,"Slope"+"_input",unitS,0) #lc units ! foo=mapOutpParam(aspectA,"Aspect"+"_input",unitS,0) #lc units ! foo=physIn3D(slopeA,"Slope"+"_input",unitS,0) foo=physIn3D(aspectA,"Aspect"+"_input",unitS,0) #print "<>getSHORE: bufferA=",bufferA loudQ=False #Finds the direction and sense of the shoreline # using the GIS 3x3 approach (as a starter) bY=sY+1;bX=sX+1 #print ">>getSHORE: elevnA,sY,sX=",sY,sX,elevnA[sY,sX] #print " bY,bX=",bY,bX #print "bufferA[]=",bufferA[bY-1:bY+2,bX-1:bX+2] zA=[] for ibY in range(-1,2): #print bufferA[bY+ibY][bX-1:bX+2] zA.append(bufferA[bY+ibY][bX-1:bX+2]) zA=NPY.ravel(zA) #(Buffer provides cells outside elevnA) #print "zA=",zA #Slope & aspect... dA=[modlRESN*60.*1836.,modlRESN*60.*1836.*NPY.cos(lolaA[1]*d2r)] b = (zA[2] + 2.*zA[5] + zA[8] - zA[0] - 2.*zA[3] - zA[6])\ / (8.*dA[1]) #dZdX c = (zA[0] + 2.*zA[1] + zA[2] - zA[6] - 2.*zA[7] - zA[8])\ / (8.*dA[0]) #dZdY if loudQ: print "b,c,dA=",b,c,dA slope = NPY.arctan(NPY.sqrt(b**2 + c**2))/d2r #Now in degrees aspect = 270.+NPY.arctan(b/c)/d2r-90.*(c/abs(c)) if loudQ: print "<>getAQUA (MODIS Aqua):" #definition and dimensions of the grid -------------------- packTYPs="f"; recordSZE=STR.calcsize(packTYPs) #print "unpack size: packTYPs,",packTYPs,STR.calcsize(packTYPs) grdSZEa=[4320,2160] grdRESa=[1./12., 1./12.] grdBOXa=[-180., -90., 180., 90.] #print "grdSZEa,grdRESa,grdBOXa=",grdSZEa,grdRESa,grdBOXa #---------------------------------------------------------- datTYPa=["chl","par","rrs","flh"] unitsA=["chlor-a mg/m^3","umol/m2/sec","sr^-1","(dimensionless)"] dt=-1 for datTYPs in datTYPa: aquaINf = open(chlDIRs+'SDS_l3m_'+datTYPs+'.dat', 'rb') dt+=1 unitS=unitsA[dt] dPi=0; datA=[] #Initializations #(Move over the modlSZx*modlSZy grid) for nY in range(modlSZy): datYa=[] #Initialize the y-th row for nX in range(modlSZx): #At 100x100 this equates to PeterB's grid seekLoLaA=[lolaA[0]+modlRESN*float(nX-modlSZx/2),\ lolaA[1]+modlRESN*float(nY-modlSZy/2)] #print "Point: seekLoLaA=",seekLoLaA seekPOSa=[\ int((float(seekLoLaA[0])-grdBOXa[0])/grdRESa[0]-0.5)+1,\ grdSZEa[1]-int((float(seekLoLaA[1])-grdBOXa[1])/grdRESa[1]+0.5)-1\ ] #print "Point: seekPOSa=",seekPOSa #*POSa are X,Y if (seekPOSa[1]<0 or seekPOSa[1]>grdSZEa[1]) or \ (seekPOSa[0]<0 or seekPOSa[0]>grdSZEa[0]): #print "seekPOSa=",seekPOSa print "getASSOCgrid: Out of grid" datYa.append(fNUL) else: #OK... grdPOS=seekPOSa[1]*grdSZEa[0]+seekPOSa[0] #print "Point in grid: dPOS,grdPOS=",dPOS,grdPOS foo=aquaINf.seek(grdPOS*recordSZE) #Finds one before gridVAL=aquaINf.read(recordSZE) #Reads actual aquaVal=STR.unpack(packTYPs,gridVAL)[0] if aquaVal<0: aquaVal=float('nan') #Set nulls to NaN datYa.append(aquaVal) #(datYa is the y-collection at the x-level) #print "elvnYa=",elvnYa datA.append(datYa) aquaINf.close() #Also make a plot in matplotlib... foo=mapOutpParam(datA,datTYPs+"_input",unitS,0) #lc units ! #Note: '_input' distinguishes these from the stocks,etc maps #Assemble the output grids ... #print "datTYPs=",datTYPs if datTYPs=="chl": chlA=datA[:] if datTYPs=="par": parA=datA[:] if datTYPs=="rrs": rrsA=datA[:] if datTYPs=="flh": flhA=datA[:] #print "<=0.: kPAR=0.121*(chlA[iY][iX])**0.428 benIrradA[0]=par*NPY.exp(-kPAR*wdI) else: benIrradA[0]=-1. #Null case (eg land) '''#Secchi Disk ---------------------------------- zSD is the Secchi Disk observed depth # Holmes formula... if wdI>5: benIrradA[1]=par*NPY.exp(-(1.44)*wdI) else: benIrradA[1]=par*NPY.exp(-(1.7/zSD)*wdI) benIrradA[2]=par*NPY.exp(-(2.6/(2.5+zSD)-0.048)*zSD) #------------------------------------------------''' #print "benIrradA: wdI,par,benIrradA",wdI,par,benIrradA irradAx.append(benIrradA[0]) #Take the Morel value only (fits with SatelliteData) #Par from MODIS are in "Einsteins m-2 day-1" # and also convert it to molar units... #Gattuso++: "Irradiance data expressed in energy units # were converted to molar units using a conversion factor # of 2.5×10**18 quanta s−1 watt−1 or 4.2 μmol photons # m−2 s−1 watt−1 (Morel and Smith, 1974)." irradA.append(irradAx) #print "irradA=",irradA #Also make a plot in matplotlib... unitS="mol photons/m^2/day" foo=mapOutpParam(irradA,"irradA"+"_input",unitS,0) foo=physIn3D(irradA,"irradA"+"_input",unitS,0) #---------------------------------------------------------------- #print "<>calcTrophics:" #Use the tanh function as used by Kleypas (ReefHab), # Gattuso++, and (originally) Bosscher & Schlager 1992. '''Fig. 5. Arbitrary P −E curve (A), diel change of net primary production (B) and changes in daily Ec as a function of daily irradiance (C) for photosynthetic organisms and communities. The primary production of an organism was calculated using the hyperbolic tangent function npp=gppmax×tanh(Ez/Ek )+ra where: npp is the rate of net primary production, gppmax is the maximum rate of gross primary production (set at 100), Ez is the benthic irradiance, Ek is the irradiance at which the initial slope of the P −E curve intersects the horizontal asymptote (set at 50 μmol photons m−2 s−1) and ra is the rate of dark respiration of the autotrophs (set at −20). Ez and Ek are in μmol photons m−2 s−1. The diel change in irradiance was modeled using a sine curve and using a photoperiod of 12 h dark and 12 h light. The rate of net primary production was calculated assuming that the rates of dark respiration of the heterotrophs and autotrophs are equal. ra was therefore simply added to the npp of the organism. In this generic example, the instantaneous Ec phot. (i.e. for the organism) and Ec comm. (i.e. for the community) are, respectively, 101 and 212 μmol photons m−2 s−1 (panels A and B). In panel (C), daily Ec (continuous lines) is twice the instantaneous Ec (dashed line) and the shaded area indicates the range of daily Ec for communities. This area is enclosed within an upper line which assumes no photoacclimation and a lower line which assumes a photoacclimation parallel to the one reported for individual organisms. The thick blue line shows the range of daily compensation irradiance. Irradiance in μmol photons m−2 d−1 is calculated using the relationship 0.0432 × irradiance (in μmol photons m−2 s−1) assuming a photoperiod of 12 h dark and 12 h light. ''' dkRespirn=-20. #ie Heterotrophic activity (ie not photosynthesis) #crIrradK=250. #See discussion in Kleypas p537 #But Min Ez in table 7 in Gattuso+2006: take average crIrradK=1. #Note: genrlGrwth is often set to 100% (as in Kleypas) genrlGrwth=100. #??Bad -- Input is daily ??? #dlyIrrad=irrad/365. #(Because input irrad is annual) #print "dlyIrrad,NPY.tanh(dlyIrrad/crIrradK)=",dlyIrrad,NPY.tanh(dlyIrrad/crIrradK ) fxnlGrwth=max(genrlGrwth*NPY.tanh(irradValue/crIrradK)+dkRespirn,0.)/100. #(CJJ: 100. is to convert from %; # Result should not be less than 0.) if loudQ: print "<'+"\n") kmlOUTf.write(''+"\n") kmlOUTf.write(' '+"\n") kmlOUTf.write(' '+prjS+''+"\n") kmlOUTf.write(' '+'numericalCARBONATE generates'+\ ''+"\n") #Put these in at top to be pasted to bottom when # main program was terminated but kml needs finishing kmlOUTf.write(' \n") kmlOUTf.write(' \n") elif kmlStatus=="map": kmlOUTf.write(' '+"\n") kmlOUTf.write(' '+nameS+''+"\n") kmlOUTf.write(' '+descnS+''+"\n") kmlOUTf.write(' '+"\n") kmlOUTf.write(' '+imageS+'.png'+"\n") kmlOUTf.write(' '+"\n") kmlOUTf.write(' '+"\n") kmlOUTf.write(' '+str(imgLocnA[1]+imgLocnA[3])+''+"\n") kmlOUTf.write(' '+str(imgLocnA[1])+''+"\n") kmlOUTf.write(' '+str(imgLocnA[0]+imgLocnA[2])+''+"\n") kmlOUTf.write(' '+str(imgLocnA[0])+''+"\n") kmlOUTf.write(' 0.0'+"\n") kmlOUTf.write(' '+"\n") kmlOUTf.write(' '+"\n") elif kmlStatus=="end": kmlOUTf.write(' '+"\n") kmlOUTf.write(''+"\n") kmlOUTf.close() return #======================================================================= def mapOutpParam(nmA,nameS,unitS,showQ): #Plots any n,m array (nmA) on a map base of the indices #The other dimension o, serves to colour the points loudQ=False if loudQ: print;print ">>mapOutpParam: nameS=",nameS print "nmA=",nmA #Note special cases if "log_wd" in nameS: nmA=map(lambda a: NPY.log10(a), nmA) #Optional #Take out the nan's for nm in range(len(nmA)): nmEs="|".join(map(str,nmA[nm])) nmEs=nmEs.replace('nan','-1.') if loudQ: print "nmEs=",nmEs nmA[nm]=map(float,nmEs.split("|")) #if "crStockSum" in nameS: # print "nmA=","\n".join(map(str,nmA)) #First make the coord and flattened data arrays... rx=NPY.arange(0.,float(modlSZx)+0.01,1.) ry=NPY.arange(0.,float(modlSZx)+0.01,1.) X,Y=NPY.meshgrid(rx,ry) fig1=PLT.figure() ax1=PLT.subplot(111) PLT.pcolor(X,Y,nmA,cmap=PLT.cm.jet) ax1.set_aspect(1.) cb=PLT.colorbar() cb.set_label(unitS) PLT.title(nameS+":"+runYMDHMs) if showQ==1: fig1.show() #Then, draw land as patch -------- patches=drawLandCells() for patch in patches: ax1.add_patch(patch) if "_input" in nameS: #Is an environmental input fig1.savefig(dirS+"_ENVIRO/"+nameS+'.png') else: #Is model output fig1.savefig(dirS+runDirS+"_MapPlots/"+nameS+'.png') fig1.clf() PLT.close() return #======================================================================= def mapDomOrgnsms(nmA,nmB,nameS,showQ): #Plots the integer n,m array (nmA) as symbols # with size (from nmB) and colour (from nmA) #The legend shows the taxon names loudQ=False if loudQ: print;print ">>mapDomOrgnsms: nameS=",nameS print "nmA=",nmA #First make the coord and flattened data arrays... rx=NPY.arange(0.5,float(modlSZx)+0.5,1.) ry=NPY.arange(0.5,float(modlSZy)+0.5,1.) X,Y=NPY.meshgrid(rx,ry) x=NPY.ravel(X); y=NPY.ravel(Y) nma=NPY.ravel(nmA) nmb=NPY.ravel(nmB) #cz=[hexColorsA[na] for na in nma] #Previous version: try to get back to it ! cz=[] for na in nma: #print na,hexColorsA[na] if na>=0: cz.append(hexColorsA[na]) else: #Deal with null cells... cz.append('gray') #print cz #This coding allows the colour scheme to adapt auuto to # number of organisms as it increases sz=NPY.array([max(0.0,float(nb))*100. for nb in nmb]) '''if tNow==60: showCoord=showA[1]*modlSZx+showA[0] print "t=60, nma,b,cz,sz=",zip(nma,nmb,cz,sz)[showCoord]''' #print "nma,b,cz,sz=",zip(nma,nmb,cz,sz) fig3=PLT.figure() ax3=PLT.subplot(111) #Main plot ----------------------------- #use nmA for size and color ax3.scatter(x,y,s=sz,c=cz,marker='o') ax3.set_aspect(1.) PLT.xlim(0,modlSZx) PLT.ylim(0,modlSZy) #print nameS+":"+runYMDHMs PLT.title(nameS+":"+runYMDHMs) #Then, draw land as patch -------- patches=drawLandCells() for patch in patches: ax3.add_patch(patch) #Move the whole lot over to left by 0.1... bbox=ax3.get_position() ax3.set_position([bbox.x0-0.1,bbox.y0,bbox.width,bbox.height]) for nC in range(numOrgnsms): aR=float(nC+1)/float(numOrgnsms) leg=PLT.text(float(modlSZx)+0.02, float(modlSZy)-2.*aR,\ str(nC)+" "+kbOrgnsmA[nC], \ color=colorsA[nC],fontsize=8) #fontname='Courier',\ #horizontalalignment='left',\ #verticalalignment='center',) #transform = ax2.transAxes,\ #) fig3.savefig(dirS+runDirS+"_MapPlots/"+nameS+'.png') fig3.clf() PLT.close() return #======================================================================= def drawLandCells(): #Uses maskYXa from global #print "maskYXa=",maskYXa patches=[] for yxA in maskYXa: fspY,fspX=yxA #print fspY,fspX #Draw the main patchwork... verts = [ (fspX, fspY), # left, bottom (fspX, fspY+1.), # left, top (fspX+1., fspY+1.), # right, top (fspX+1., fspY), # right, bottom (fspX, fspY), # (Closure) ignored ] codes = [PATH.Path.MOVETO, PATH.Path.LINETO, PATH.Path.LINETO, PATH.Path.LINETO, PATH.Path.CLOSEPOLY, ] path = PATH.Path(verts, codes) #The drawing... patches.append(PATCHES.PathPatch(path, \ facecolor="0.9", \ edgecolor="0.9", lw=1)) return patches #A drawing structure #======================================================================= def graphStockSeqnce(a,bA,nameS,showQ): #Used to plot all the separate taxapop histories on one graph loudQ=False if loudQ: print ">>graphStockSeqnce: nameS=",nameS print "len(a),NPY.array(bA).shape=",len(a),NPY.array(bA).shape #a=range(-1,tStop+1) #A fix bA=NPY.array([b[1:] for b in bA]) #A fix, also if len(a)<>len(bA[0]): print "len(a),NPY.array(bA).shape,len(nameS)",\ len(a),NPY.array(bA).shape,len(nameS) print "Warning: array mismatch: lengths|a,bA[0] " return oA=kbOrgnsmA[:] fig4 = PLT.figure() ax4=PLT.subplot(111) labelA=['']*len(oA) lineA=['']*len(oA) for bL in range(len(oA)): #print "bN,bL,a,hexColorsA[bL]=",numOrgnsms,bL,a,hexColorsA[bL] #print "bA[bL]=",bA[bL],len(a),numOrgnsms lineA[bL]=PLT.plot(a, bA[bL], \ color=hexColorsA[bL], \ label=oA[bL]) #Add the vacant bed plot ... vacA=[] stockByTa=NPY.transpose(bA) for sBTa in stockByTa: vac=0. for sBT in sBTa: vac+=10.**sBT vacA.append(max(1.-vac,stockNull)) if logPlotQ: vacA=NPY.log10(vacA) #Else unchanged lineA.append(PLT.plot(a, vacA,'k:',label="bare_ground")) oA.append("bare_ground") if loudQ: print "vacA,oA=",vacA,oA PLT.title(nameS+"_"+runYMDHMs) unitsA=["Areal Stock (fractional)","Time (yrs)"] ax4.set_xlabel(unitsA[0]) ax4.set_xlabel(unitsA[1]) #Follow # "http://stackoverflow.com/questions/4700614/how-to-put-the-legend-out-of-the-plot" #Note: ylim not necessary if logPlotQ: ax4.set_ylim(-4.,1.) #log10 values else: ax4.set_ylim(0.,1.5) #Linear value plotting # Shrink current axis's height by 10% on the bottom... box = ax4.get_position() ax4.set_position([box.x0, box.y0 + box.height * 0.33, box.width, box.height * 0.67]) #Legend ------------------------- leg=ax4.legend(lineA,oA,loc='upper center',labelsep=0,\ bbox_to_anchor=(0.5, -0.2),ncol=3) for t in leg.get_texts(): t.set_fontsize(6) # the legend text fontsize fig4.savefig(dirS+runDirS+"_GraphPlots/"+nameS+'.png') if showQ==1: fig4.show() #print "<>calcMilieu: enviroA=",enviroA #(These come in by global) print "iCr,kbWdHabittA,kbTempHabittA,kbParHabittA=",\ iCr,kbWdHabittA[iCr],kbTempHabittA[iCr],kbParHabittA[iCr] suitabA=[]; suitabSa=[]; crTrophic=[-99.]*3 #Note: crTrophic is an array of [photosyntetic,autonutrient,particulate] # adjustments of growth potential (CJJ May 2011) for enviro in enviroA: #ie Each parameter... if loudQ: print "enviro=",enviro #enviro holds the name and the value #Will assign a suitability only if there is a # confluence of the enviro and organism inputs. # Therefore important to have organism tolerances in KB. if enviro[0]=="wd" and kbWdHabittA[iCr]<>0.: #Water depth... intHabiR=interpHabiRange(enviro,kbWdHabittA[iCr]) suitabA.append(intHabiR) suitabSa.append("wd") elif enviro[0]=="temp" and kbTempHabittA[iCr]<>0.: #Temperature... intHabiR=interpHabiRange(enviro,kbTempHabittA[iCr]) suitabA.append(intHabiR) suitabSa.append("temp") elif enviro[0]=="par" and kbParHabittA[iCr]<>0.: if loudQ: print "At par: kbParHabittA[iCr],enviro=",\ kbParHabittA[iCr],enviro #Light suitability... intHabiR=interpHabiRange(enviro,kbParHabittA[iCr]) suitabA.append(intHabiR) suitabSa.append("par") #Also, photosynthetically available radiation (light)... # will adjust the KB growth rate (potential) proportionally crTrophic[0]=calcTrophics(enviro,iCr) '''elif enviro[0]=="chl" and kbNutrHabittA[iCr]<>0.: #Auto-nutrient trophic acalar... pass elif enviro[0]=="nutr" and kbChlHabittA[iCr]<>0.: #Auto-nutrient trophic acalar... pass''' if loudQ: print "Suitability: iCr,suitabA,suitabSa",\ iCr,suitabA,suitabSa while -99. in suitabA: suitabA.pop(suitabA.index(-99.)) if suitabA==[]: suitabA=[0.] suitabA=[max(sA,0.) for sA in suitabA] #This is an experimental situation ... #(? Use logit instead ?) #1. crSuitab=NPY.median(suitabA) #2. #crSuitab=min(suitabA) #3. #crSuitab=NPY.average(suitabA) if loudQ: print "<>interpHabiRange: enviro,habitatS=",enviro,habitatS enviroParam,enviroPrpty=enviro #Name and value #Detect KB errors ------------------------------------------------ #(i) Intercept null cases (eg [-99,-99])... #(This applies where a parameter is irrelevant to an organism, # like light for Lophelia pertusa) if ",-99" in habitatS: #print "<2: #An entry does not consist of a param/abund pair... raw_input("Unpaired KB input data (fix KB):"+habitatS) sys.quit("Stopped") #print habitatA #Ok to transfer to an array now... habValuMnMx=[min([hab[0] for hab in habitatA]),\ max([hab[0] for hab in habitatA])] habPcntMnMx=[min([hab[1] for hab in habitatA]),\ max([hab[1] for hab in habitatA])] if loudQ: print "habValuMnMx,habPcntMnMx=",habValuMnMx,habPcntMnMx #Raise max% to 100%... if habPcntMnMx[1]<100.: #print "habPcntMnMx[1]=",habPcntMnMx[1],habPcntMnMx #print "Pre: habitatA=",habitatA habitatA=[[habA[0],habA[1]/habPcntMnMx[1]*100.] for habA in habitatA] if loudQ: print "Post normalizing: habitatA=",habitatA habitatA=fixHabitatArray(habitatA) #Repair of the ends of the arrays enviroPrpty=float(enviroPrpty) if loudQ: print "habitatA=",habitatA print "enviroPrpty=",enviroPrpty #Now do the interpolation ------------------------------------- #It finds the habitat suitability in terms of one enviro variable if enviroPrpty>=habitatA[-1][0]: #If greater than highest, adopt zero #habiSuitA=habitatA[-1][1] habiSuitA=0. elif enviroPrpty<=habitatA[0][0]: #If less than lowest, adopt zero habiSuitA=0. else: for kkPi in range(1,len(habitatA)): #print kkPi,enviroPrpty,habitatA[kkPi][0] if enviroPrpty<=habitatA[kkPi][0]: #Lies between pair... #print "Gotcha: enviroPrpty,habitatA[kkPi-1,kkPi]",\ # enviroPrpty,habitatA[kkPi-1],habitatA[kkPi] habiSuitA=habitatA[kkPi][1]-(habitatA[kkPi][0]-enviroPrpty)* \ (habitatA[kkPi][1]-habitatA[kkPi-1][1])/ \ (habitatA[kkPi][0]-habitatA[kkPi-1][0]) break habiSuitA=habiSuitA/100. #Change from percent if loudQ: print "<0.: #Right side of array has an unclosed abundance # (should really go to zero) lastHabi=[habitatA[-2][0],habitatA[-1][0]] newValue=max(lastHabi[1]+habitatA[-1][1]/100.*\ (lastHabi[1]-lastHabi[0]),-1.) rgtHabitatA=[newValue,0.] #Note difference from [[]] lftHabitatA didLRq[1]=True if loudQ: print "Unclosed right in KB lineset (0-added):",\ lastHabi,rgtHabitatA #Pull it all together... if loudQ: print "Finishing fixing: ",list(habitatA),\ lftHabitatA,rgtHabitatA,didLRq if didLRq[0]: lftHabitatA.extend(list(habitatA)) habA=lftHabitatA[:] else: habA=habitatA[:] #print habA,didLRq if didLRq[1]: habA.append(rgtHabitatA) #(The loop deeply converts the array to list) if loudQ: print "After fixing: habA=",habA,didLRq habitatA=NPY.array(habA) if loudQ: print "End of error fixing: habitatA=",habitatA return habitatA # ============================================================================ def themAndUs(iCr,stocksNowInCellA,suitabsInCellA,trophicsInCellA): #This organizes all the iCr parameters into a Them/Us system... #The inputs come through global #These are the steady parameters (not competn, initial) # but also include eventsA #Growth is anything that increases the number of fecund individuals. # So increasing colony size, number of colonies both count. #Units and definitions... # growth is in fxn of colony size/yr # asperity is in fxn of colony area # suitab is in fxn of 100% areal coverage # recruit is in fxn of pop # competn is fxn via the interaction matrix loudQ=False if showCellQ and loudQ: print;print ">>themAndUs: kbOrgnsmA[iCr]=",kbOrgnsmA[iCr] print " iCr,stocksNowInCellA=",iCr,stocksNowInCellA print " suitabsInCellA,trophicsInCellA=",\ suitabsInCellA,trophicsInCellA #print "kbGrowthA=",kbGrowthA growthA=[0.,0.]; suitabA=[0.,0.]; competnA=[0.,0.] immigrA=[0.,0.]; mortalA=[0.,0.]; lvefxnA=[0.,0.] xA=[0.,0.] popGrA=[]; popSuA=[]; popCoA=[]; popMoA=[]; popFfA=[]; popImA=[] for iCrKB in range(numOrgnsms): #print "stocksNowInCellA=",stocksNowInCellA crStockNowInCellA=stocksNowInCellA[iCrKB] #Assign KB values to the equation variables ================== if iCrKB==iCr: #Us ! #Calculate basic environmental effects ------------------- #Growth: #i. Annulus is increase in area # (linear_growth+colony radius) grwthArea,colnyArea=calcAnnulus(kbGrowthA[iCrKB],kbCoverageA[iCrKB][0]) growthA[0]=(grwthArea+\ kbReprodA[iCrKB][0][0]*kbReprodA[iCrKB][0][1])*\ kbCoverageA[iCrKB][1] #ie: (growth_annulus+stock*clone_numbr*clone_success)*lvefxn #print "iCrKB,grwthArea,colnyArea,kbReprodA[iCrKB],kbCoverageA[iCrKB]=",\ # iCrKB,grwthArea,colnyArea,kbReprodA[iCrKB],kbCoverageA[iCrKB] #(Area of growth incl clones) #growthA finishes as a fractional value #Suitabilities: suitabA[0]=suitabsInCellA[iCrKB] #Is fraction /1. #Frame fraction lvefxnA[0]=kbCoverageA[iCrKB][1] #Only used for immigrants #Mortalities: mortalA[0]=max(kbMortalityA[iCrKB][0],1./kbMortalityA[iCrKB][1]) #print "mortality-us: mortalA[0],iCr=",mortalA[0],kbMortalityA[iCrKB],iCr #(In terms of rate=fractional) # [absolute#/m2 per stock unit] colonization) xA[0]=crStockNowInCellA else: #Them ! #Calculate basic environmental effects ------------------- grwthArea,colnyArea=calcAnnulus(kbGrowthA[iCrKB],kbCoverageA[iCrKB][0]) grow=(grwthArea+\ kbReprodA[iCrKB][0][0]*kbReprodA[iCrKB][0][1])*\ kbCoverageA[iCrKB][1] #(Ie: % to real, then Linear to areal) suit=suitabsInCellA[iCrKB] mort=max(kbMortalityA[iCrKB][0],1./kbMortalityA[iCrKB][1]) #ie max of (Native %,Reciprocal of longevity) #print "mortality-other: mort,iCrKB=",mort,kbMortalityA[iCrKB],iCrKB #Frame fraction lvefxn=kbCoverageA[iCrKB][1] #(This weighted-sums, preparing for the averages) popGrA.append(crStockNowInCellA*grow) popSuA.append(crStockNowInCellA*suit) popMoA.append(crStockNowInCellA*mort) popFfA.append(crStockNowInCellA*lvefxn) xA[1]+=crStockNowInCellA #Use median for growth, suitabA and recruit... if xA[1]>0.: growthA[1]=NPY.average(popGrA)/xA[1] suitabA[1]=NPY.average(popSuA)/xA[1] mortalA[1]=NPY.average(popMoA)/xA[1] lvefxnA[1]=NPY.average(popFfA)/xA[1] #print "mortality-rest: mortalA[1]=",mortalA[1] else: growthA[1]=0.; suitabA[1]=0.; mortalA[1]=0.; lvefxnA[1]=0. #Update recruitments with local (added to external)... immigrA=calcImmi(xA,immigrA,suitabA) #Competition methods ... #Get the competition from populations and colony sizes --- #Assume everybody is a bad neighbour and that the # interaction is proporional to impingement. #Think of a 1m2 quadrat. #Method A: fixed compA=[0.0,0.05] #Paranoid model all rest are 5% against 'us' #Was set to 5% from 20% by CJJ Apr 2011 competnA=compA[:] '''Too fancy ! #Method B: geometric (by ratios of circumferences) # bigger is more competitive competnA=[0.]*2 #Initializing compA=[0.0,0.05] #Assumed competition biases for iCrKB in range(numOrgnsms): if iCrKB==iCr: #Us competnA[0]=compA[0]*kbCoverageA[iCr][0]/kbCoverageA[iCrKB][0] #ie Competition by us [0] against all the rest [1] by iCr else: #Them popCoA.append(compA[1]*kbCoverageA[iCrKB][0]/kbCoverageA[iCr][0]) #ie Competition against iCr [0] by all rest # (if iCr is small; this value goes up) competnA[1]=NPY.average(popCoA)''' if showCellQ and loudQ: print;print "growthA=",growthA,"(Area)" print "suitabA=",suitabA,"(Fractional)" print "lvefxnA=",lvefxnA,"(Fractional)" print "immigrA=",immigrA,"(Area)" print "mortalA=",mortalA,"(Fractional?)" print "competnA=",competnA,"(Fractional)" print "xA=",xA,"(Area)" #raw_input("show me") #foo=checkSumStocks(stocksNowInCellA,"endOfThemAndUs") #Each of these is a size two list... return growthA,suitabA,lvefxnA,immigrA,mortalA,competnA,xA # ============================================================================ def calcAnnulus(linrGrwth,diam): colnyArea=pi*(diam/2.)**2 #Note: kbCoverageA 0= diameter (m); # 1=live_area (fxn); 2=capacity (K? #/m2) grwthArea=colnyArea-pi*(diam/2.-linrGrwth)**2 # ie. area(radius-(radius-linear_growth) ); that is, # the growth always attains UP TO the final size) return grwthArea,colnyArea # ============================================================================ def applyEvent(iCrKB,stocksNowInCellA): loudQ=False #Gets values from the table... evMortA,evImmiA=[[0.,0.],[0.,0.]] #Defaults for them&us evMort=0.; evImmi=0. eventQ=False for evA in eventsA: #Like: ["disease",5,0.2,[2,4],0,[],-,1] if tNow % evA[1]==0 and tNow<>0: eventQ=True #An event coincides ! for evCr in evA[3]: '''if showCellQ and loudQ: print "Applied -ve event: tNow,iCrKB,evCr,evA=",\ tNow,iCrKB,evCr,evA''' if evCr==iCrKB: '''if showCellQ and loudQ: print "evMortA[0] b4 after:",\ evMortA[0],evMortA[0]+evA[2]''' evMortA[0]+=evA[2] #This allows for more than 1 event being possible # (will take average) else: evMort+=evA[2]*stocksNowInCellA[evCr] for evCr in evA[5]: '''if showCellQ and loudQ: print "Applied +ve event: tNow,iCrKB,evCr,evA=",\ tNow,iCrKB,evCr,evA''' if evCr==iCrKB: '''if showCellQ and loudQ: print "evImmiA[0] b4 after:",\ evImmiA[0],evImmiA[0]+evA[4]''' evImmiA[0]+=evA[4] #This allows for more than 1 event being possible # (will take average) else: evImmi+=evA[4] if eventQ: #Get (un/weighted) averages ... #print "stocksNowInCellA & pop=",\ # stocksNowInCellA,stocksNowInCellA.pop(iCrKB) themStocks=sum(stocksNowInCellA)-stocksNowInCellA[iCrKB] if themStocks>0.: evMortA[1]=evMort/themStocks else: evMortA[1]=0. #ie averaged over the sum of 'them' stocks #if loudQ: print "Mort Event:",evCr,evMortA,evImmiA evImmiA[1]=evImmi/float(numOrgnsms-1) #Averaged over count of 'them' organisms #if loudQ: print "Immi Event:",evCr,evMortA,evImmiA #Creates arrays [evMortUs,evMortThem],[evImmiUs,evImmiThem] if showCellQ and loudQ: #print " evMort,evImmi=",evMort,evImmi print " evMortA,evImmiA",evMortA,evImmiA return evMortA,evImmiA # ============================================================================ def calcImmi(popA,immigrA,suitabA): #For time being assume only in-cell recruitment (CJJ Oct10) popA=[max(pop,stockNull) for pop in popA] #None allowed to be zero vacancy=calcVacancy(popA) #Fractional area immigrA=[immigrA[0]+vacancy*suitabA[0]*popA[0],\ immigrA[1]+vacancy*suitabA[1]*popA[1]] #Is by us/them #print "<>calcLotStocks:" print "iCr,kbOrgnsmA,stocksNowInCellA,suitabsInCellA,trophicsInCellA=" for iCr in range(numOrgnsms): print iCr,kbOrgnsmA[iCr],stocksNowInCellA[iCr],\ suitabsInCellA[iCr],trophicsInCellA[iCr] #By limiting tA, we limit the lenth of the burst... tA=map(float,range(tBurstRun)) stockSeriesInCellA=[]; deadSeriesInCellA=[] #List of stock through time for each organism dim=tBurstLen*numOrgnsms #Note: crStockSeqA=[Stock of organism,Stock of the rest]*tBurstLen # stockSeriesInCellA time-arrays of stocks per organism # stocksNowInCellA=Array of Stock of each organism per cell for iCr in range(numOrgnsms): foo=checkSumStocks(stocksNowInCellA,">>calcLotStocks") #Get the coefficients ----------------- #Unchanging in the burst-run growthA,suitabA,lvefxnA,immigrA,mortalA,competnA,xA\ =themAndUs(iCr,stocksNowInCellA,\ suitabsInCellA,trophicsInCellA) #xA is the starting stock set #(Note: the order of these if's is important) if suitabA[1]>> 'Integration successful.' if infodict['message']<>"Integration successful.": print;print infodict['message']+"!"*10 print "iCr,nY,nX,tNow=",iCr,nY,nX,tNow print "xA,growthA,suitabA,competnA,mortalA,immigrA=",\ xA,growthA,suitabA,competnA,mortalA,immigrA raw_input("Was gibt ?") #Send to output array (and put floor value in to help graphing) crStockSeqA=[max(nS,stockNull) for nS in newStockA[:,0]] #print "crStockSeqA",crStockSeqA #(Note: above op changes crStockSeqA to a list.) else: crStockSeqA=[stockNull]*tBurstRun #print "Just set: crStockSeqA",crStockSeqA if loudQ: print "iCr,crStockSeqA.shape,crStockSeqA=",\ iCr,NPY.array(crStockSeqA).shape,crStockSeqA #stockSeriesInCellA accumulates each organism's time series in a row... stockSeriesInCellA.append(crStockSeqA) if loudQ: print "[mortalA[0]]*tBurstLen=",[mortalA[0]]*tBurstLen print "addMortA[0],stockSeriesInCellA[iCr]",\ addMortA[0],stockSeriesInCellA[iCr] deadSeriesInCellA.append(\ [sSIC*mortalA[0] for sSIC in stockSeriesInCellA[iCr]]) #[Note: deadSeriesInCellA keeps the tally for deposits...] #The first entry holds the event mort (if any) ... #print "deadSeriesInCellA[iCr]=",deadSeriesInCellA[iCr] deadSeriesInCellA[iCr][0]+=addMortA[0]*stockSeriesInCellA[iCr][0] if loudQ: print "deadSeriesInCellA[iCr]=",deadSeriesInCellA[iCr] print "stockSeriesInCellA[iCr]=",stockSeriesInCellA[iCr] stockSetsInCellA=NPY.array(stockSeriesInCellA).transpose() deadSetsInCellA=NPY.array(deadSeriesInCellA).transpose() #print "stockSetsInCellA.shape=",stockSetsInCellA.shape #print "deadSetsInCellA.shape=",deadSetsInCellA.shape #If the max sum in the series is >1.0 # can do some normalization... #stockSeriesInCellA=normaliseBurstPops(stockSeriesInCellA) #CJJ to check whether this is needed/beneficial... if loudQ: print "<=1.0: #On land ! maskYXa.append([nY,nX]) if loudQ: print "Land mask set: maskYXa=",maskYXa maskYXa=NPY.array(maskYXa) #Note: maskYXa goes out through global return # ============================================================================ def reportEnviro(tempA,chlA,parA,irradA,zEnrgyDissiA): loudQ=False if loudQ: print;print ">>reportEnviro: tNow=",tNow #All the KB biol data come in through global repOUTf=open(dirS+runDirS+"_StockDumps/"+\ "_enviro_"+str(tNow)+".rep","w") #Initial file (tStart=0) #Stage A: General repOUTa=["\n"] #Starts, and will also serve as a line space repOUTa.append("#File=Population state file for CarboLOT") repOUTa.append("#Project code|name="+":".join(prjA)) repOUTa.append("#Run date/person="+runYMDHMs+"LOCAL by CJJ") repOUTa.append("") repOUTa.append("#General "+"-"*60) repOUTa.append(str(tNow)+"\t"+"=Time Now (yrs)") repOUTa.append(",".join(map(str,lolaA))+"\t"+\ "=Model Centre (degr WGS84)") repOUTa.append(str(modlRESN)+"\t"+\ "=Model Resolution (degr WGS84)") repOUTa.append(",".join(map(str,LoLfHtWdA))+"\t"+\ "=Model LL LoLa (degr WGS84)") repOUTf.write("\n".join(map(str,repOUTa))+"\n") #Stage B: EnviroGrids repOUTa=["\n"] #Will serve as a line space repOUTa.append("#Enviro Grids "+"-"*60) #Elevation ----------------- repOUTa.append(",".join(map(str,modlSZa))+"\t"+\ "=Elevations (m)") for elevn in elevnA: repOUTa.append(",".join(map(str,elevn))+";") repOUTa.append("#"+"-"*60) repOUTa.append("\n") #Bottom Temp ----------------- repOUTa.append(",".join(map(str,modlSZa))+"\t"+\ "=Bottom temperature (degC)") for temp in tempA: repOUTa.append(",".join(map(str,temp))+";") repOUTa.append("#"+"-"*60) repOUTa.append("\n") #Surface PAR ----------------- repOUTa.append(",".join(map(str,modlSZa))+"\t"+\ "=Surface PAR Irradiance (moles/m2/day)") for par in parA: repOUTa.append(",".join(map(str,par))+";") repOUTa.append("#"+"-"*60) repOUTa.append("\n") #Chl ----------------- repOUTa.append(",".join(map(str,modlSZa))+"\t"+\ "=Surface Chlorophyll (mg/ltr)") for chl in chlA: repOUTa.append(",".join(map(str,chl))+";") repOUTa.append("#"+"-"*60) #Benthic PAR ----------------- repOUTa.append(",".join(map(str,modlSZa))+"\t"+\ "=Benthic PAR Irradiance (moles/m2/day)") for irrad in irradA: repOUTa.append(",".join(map(str,irrad))+";") repOUTa.append("#"+"-"*60) repOUTa.append("\n") #Energy ----------------- repOUTa.append(",".join(map(str,modlSZa))+"\t"+\ "=Wave Energy Dissip at Z (J/m2)") for zEnrgyDissi in zEnrgyDissiA: repOUTa.append(",".join(map(str,zEnrgyDissi))+";") repOUTa.append("#"+"-"*60) #Put the output ... repOUTf.write("\n".join(map(str,repOUTa))+"\n") repOUTf.close() #print "<>calcSedProduction: " #Sum the linear_growth_rates*calcification for iCr in range(numOrgnsms): #prodConst = 0.01 #From PeterB # This is a factor to transform a population number into a # volume of produced carbonate sediment #prodYXa[nY][nX]=creatrA[nY][nX] * prodConst #if showCellQ and loudQ: print "stocksA=",stocksA grwthArea,colnyArea=calcAnnulus(kbGrowthA[iCr],kbCoverageA[iCr][0]) crIncrse=stocksA[iCr]*grwthArea/colnyArea*kbCoverageA[iCr][1] #fxnl area increase #Multiply by crunched vol (kbSkltlA[:][0]) so reduced # to solid level and correct for porosity and density crProdn=crIncrse*0.005*kbSkeletalA[iCr][0]*(skltlPorosity*2870.) #0.005 represents the 5mm live 'skin' # ie. stock*area_growth*frame_fraction*porosity*mineral_density if showCellQ and loudQ: if crProdn>1.0: print "prodn too big: iCr,grwthArea,colnyArea=",iCr,grwthArea,colnyArea print " stocksA[iCr]=",stocksA[iCr] print " kbCoverageA[iCr](sz,lvefxn,capac)=",kbCoverageA[iCr] print " kbGrowthA[iCr](grw)=",kbGrowthA[iCr] print " kbSkeletalA[iCr][0](disintVol)=",kbSkeletalA[iCr][0] print "crIncrse,crProdn,kbSkeletalA[iCr][0],kbBinPresncA=",\ crIncrse,crProdn,kbSkeletalA[iCr][0],kbBinPresncA #Distribute the production to the zBins for iBin in range(4): prodnA[iBin]+=crProdn*kbBinPresncA[iBin] # ie organism growth in each z-bin (growth is per year) # in kg/m2/yr (hence 2870. density) if showCellQ and loudQ: print "<>calcSedBuildup: deadsA",deadsA #Sum the linear_growth_rates*calcification for iCr in range(numOrgnsms): if showCellQ and loudQ: print "iCr,deadsA[iCr]=",iCr,deadsA[iCr] print " kbSkeletalA[iCr](disintVol,disintGrz,-,frFxn)=",kbSkeletalA[iCr] #The thickness is the 'stock' of deads (area) # * framefraction * crunched sizes # to give the thickness (metres) crThkns=deadsA[iCr]*kbSkeletalA[iCr][3]*\ kbSkeletalA[iCr][0]*sum(kbBinPresncA[1:]) #(Note: Use only the above-surface bins) #Detritus is made at and above bed by bioerosion # at a rate of (say) 20% of the thickness per year # (see "http://www.jstor.org/stable/pdfplus/2460500.pdf?acceptTC=true") thknsA[0]+=crThkns*0.80; thknsA[1]+=crThkns*0.20 #(Note later this can be adjusted by the populations present of # bioeroders based on their hab suitab) #Will use dbSEABED syntax for lithology PET descns... '''Ultimately will produce strings like these ... insituPetA+="cr_"+str(iCr)+" -"+str(kbSkeletalA[iCr][1])+"m_szd"+\ " /"+str(thknsA[0])+"m_thk"+" ; " detriPetA+="cr_"+str(iCr)+"_sed"+" -"+str(detriGrz)+"m_szd"+\ " /"+str(thknsA[1])+"m_thk"+" ; "''' insituPetA.append([iCr,kbSkeletalA[iCr][1],thknsA[0]]) #dbSBD LTH detriGrz=0.5/1000. #(Note: interim grainsize m, =0.5mm =1 phi) detriPetA.append([iCr,detriGrz,thknsA[1]]) #dbSBD CMP if sum(thknsA)>1.0: #if showCellQ and loudQ: print "thkns too thick: thknsA,insituPetA,detriPetA=",\ thknsA,insituPetA,detriPetA print " deadsA=",deadsA return thknsA,insituPetA,detriPetA # ============================================================================ # ============================================================================ # Read in Knowledge Base def getKBdata(): global kbCoverageA,kbGrowthA,kbSkeletalA,\ kbZonesA,kbMortalityA,kbReprodA,\ kbWdHabittA,kbTempHabittA,kbParHabittA global kbCoverageU,kbGrowthU,kbSkeletalU,\ kbZonesU,kbMortalityU,kbReprodU,\ kbWdHabittU,kbTempHabittU,kbParHabittU global numOrgnsms,kbOrgnsmA kbInF=open("../_carboKB/"+kbFileName) #Presently kbFileName="_carboKB_7.okb" loudQ=False #Initializations... kbCoverageA=[]; kbGrowthA=[]; kbSkeletalA=[];\ kbZonesA=[]; kbMortalityA=[]; kbReprodA=[];\ kbWdHabittA=[]; kbTempHabittA=[];kbParHabittA=[] kbCoverageU=[]; kbGrowthU=[]; kbSkeletalU=[];\ kbZonesU=[]; kbMortalityU=[]; kbReprodU=[];\ kbWdHabittU=[]; kbTempHabittU=[];kbParHabittU=[] numOrgnsms=0; kbOrgnsmA=[] while kbInF: kbInS=kbInF.readline() if kbInS=="": break #EOF kbInS=kbInS.strip() if kbInS=="" or kbInS[0]=="#": pass else: kbInA=kbInS.replace("[","").replace("]","").split("\t") kbInA=[kIa.strip() for kIa in kbInA] #print "kbInA=",kbInA #Here we use the same variable names as # in the write statements in _carboKB) '''"Organism\t"+str(iCr)+"|"+kbOrgnsmA[iCr] "\n Covrg\t"+"|".join(map(str,kbCoverageA[iCr]))+\ "\t"+"|".join(map(str,kbCoverageU[iCr]))+\ "\n Grwth\t"+"["+str(kbGrowthA[iCr])+"]"+\ "\t"+"["+str(kbGrowthU[iCr])+"]"+\ "\n Skltl\t"+"|".join(map(str,kbSkeletalA[iCr]))+\ "\t"+"|".join(map(str,kbSkeletalU[iCr]))+\ "\n Zones\t"+"|".join(map(str,kbZonesA[iCr]))+\ "\t"+"|".join(map(str,kbZonesU[iCr]))+\ "\n HabWd\t"+str(kbWdHabittA[iCr])+\ "\t"+str(kbWdHabittU[iCr])+\ "\n HabTemp\t"+str(kbTempHabittA[iCr])+\ "\t"+str(kbTempHabittU[iCr])+\ "\n HabPar\t"+str(kbParHabittA[iCr])+\ "\t"+str(kbParHabittU[iCr])+\ "\n MortR\t"+"|".join(map(str,kbMortalityA[iCr]))+\ "\t"+"|".join(map(str,kbMortalityU[iCr]))+\ "\n ReprR\t"+"|".join(map(str,kbReprodA[iCr]))+\ "\t"+"|".join(map(str,kbReprodU[iCr]))+"|"+"\n"''' #Note ! All parameters must have something appended # for all organisms; else arrays go out of sync. if kbInA[0]=="Organism": kbOrgnsmA.append(kbInA[1].split("|")[1]) iCr=int(kbInA[1].split("|")[0]) numOrgnsms=max(numOrgnsms,iCr+1) if kbInA[0]=="Covrg": kbCoverageA.append(map(float,kbInA[1].split("|"))) kbCoverageU.append(kbInA[2].split("|")) if kbInA[0]=="Grwth": #(Only a single value here) kbGrowthA.append(float(kbInA[1])) kbGrowthU.append(kbInA[2]) if kbInA[0]=="Skltl": #(Here the last item is string) kbSkelA=kbInA[1].split("|") kbSkeletalA.append(\ [float(kbSkelA[0]),float(kbSkelA[1]),\ kbSkelA[2],float(kbSkelA[3])] ) kbSkeletalU.append(kbInA[2].split("|")) if kbInA[0]=="Zones": kbZonesA.append(map(float,kbInA[1].split("|"))) kbZonesU.append(kbInA[2].split("|")) if kbInA[0]=="HabWd": kbWdHabittA.extend(kbInA[1].split("|")) #String kbWdHabittU.extend(kbInA[2].split("|")) if kbInA[0]=="HabTemp": kbTempHabittA.extend(kbInA[1].split("|"))#String kbTempHabittU.extend(kbInA[2].split("|")) if kbInA[0]=="HabPar": kbParHabittA.extend(kbInA[1].split("|"))#String kbParHabittU.extend(kbInA[2].split("|")) if kbInA[0]=="MortR": kbMortalityA.append(map(float,kbInA[1].split("|"))) kbMortalityU.append(kbInA[2].split("|")) if kbInA[0]=="ReprR": #Note: kbReprod are double-level lists kbRpA=kbInA[1].split("|") kbRpA=[map(float,kbRa.split(",")) for kbRa in kbRpA] print "kbRpA=",kbRpA kbRpU=kbInA[2].split("|") kbRpU=[kbRu.split(",") for kbRu in kbRpU] print "kbRpU=",kbRpU #CJJ: Need to intervene here to tame the recruit numbers #We will assume for most species that 3 recruits/m2 is the norm kbRpA[1]=[100.,0.03] #100 individuals produced, 3 succeed kbRpU[1]=["#/yr", "fxn/#/m2"] #works out to num successes per m2 kbReprodA.append(kbRpA) kbReprodU.append(kbRpU) #(All entries remain string format) #Echo (in exactly same format # & with same code as in carboKB_*.py) ... for iCr in range(numOrgnsms): outA=["\n"] outA.append("#"+"-"*60) outA.append("\nOrganism\t"+str(iCr)+"|"+kbOrgnsmA[iCr]) outA.append("\n Covrg\t"+repr(kbCoverageA[iCr])) outA.append("\t"+repr(kbCoverageU[iCr])) outA.append("\n Grwth\t"+"["+repr(kbGrowthA[iCr])+"]") outA.append("\t"+"["+repr(kbGrowthU[iCr])+"]") outA.append("\n Skltl\t"+repr(kbSkeletalA[iCr])) outA.append("\t"+repr(kbSkeletalU[iCr])) outA.append("\n Zones\t"+repr(kbZonesA[iCr])) outA.append("\t"+repr(kbZonesU[iCr])) outA.append("\n HabWd\t"+repr(kbWdHabittA[iCr])) outA.append("\t"+repr(kbWdHabittU[iCr])) outA.append("\n HabTemp\t"+repr(kbTempHabittA[iCr])) outA.append("\t"+repr(kbTempHabittU[iCr])) outA.append("\n HabPar\t"+repr(kbParHabittA[iCr])) outA.append("\t"+repr(kbParHabittU[iCr])) outA.append("\n MortR\t"+repr(kbMortalityA[iCr])) outA.append("\t"+repr(kbMortalityU[iCr])) outA.append("\n ReprR\t"+repr(kbReprodA[iCr])) outA.append("\t"+repr(kbReprodU[iCr])+"|"+"\n") print "".join(outA) logF.write("".join(outA)) print "Finished input of OKB !" kbInF.close() return # ============================================================================ # ============================================================================ foo=myMAIN() finiClock=time.clock() setupClock=(startClock-initClock)/60. #Minutes runClock=(endClock-startClock)/60. finishClock=(finiClock-endClock)/60. tAdvantage=tStop/(runClock/60./24/365.25) #ModelYears/RunYears outS="Elapsed times: setupClock,runClock,finishClock,tAdvantage="+\ repr(setupClock)+"|"+\ repr(runClock)+"|"+\ repr(finishClock)+"| "+\ "mins ( x"+repr(int(tAdvantage))+" advantage)" print outS logF.write(outS+"\n") #Finish ----------------------- print "End of Simulation "+"="*50 logF.close()