diff --git a/RMA2SERAPHIN_3DRMA11.py b/RMA2SERAPHIN_3DRMA11.py index 92afcc4..09f1eb5 100644 --- a/RMA2SERAPHIN_3DRMA11.py +++ b/RMA2SERAPHIN_3DRMA11.py @@ -9,6 +9,8 @@ import math from py_rmatools import rma import re from datetime import datetime, timedelta +import pandas as pd +import numpy as np plt.rcParams.update({'figure.max_open_warning': 0}) @@ -16,14 +18,16 @@ plt.rcParams.update({'figure.max_open_warning': 0}) meshFilename = 'bub005c.rm1' channelWidth = 100 -RMAfilename = 'BUB023_WQ' +RMAfilenames = ['BUB028','BUB029','BUB030','BUB031','BUB032'] #If RMA11 constNum = [1] #If RMA11 3D -Options = ['Max',1,4,8] #Option available 'Max' , an integer correspond to the layer number to extract the results, if 2D use options = [1] +#Options = [4,8,1,'Max'] #Option available 'Max' , an integer correspond to the layer number to extract the results, if 2D use options = [1] - 'Percentile' +Options = ['Percentile'] +percentiles = [0.5,0.8,0.9,0.95,1] NodeLayersFile = 'NodeLayers.txt' #the node layer file was created by getting the nodes number from the RMA Outfile (just run RMA11 for a few steps) Nnodes = 5769 #Number of nodes (including mid side node) - thisn umber also correspond to the last node number in the nodelayers file # In[66]: @@ -89,7 +93,7 @@ def oneD2triangle(ElementNum): ElementDict[newEle] = [nA, nB, nBe] def RMA11toSerafin(option=1): - f = open('{}_{}.slf'.format(RMAfilename,option), 'wb') + f = open('{}_{}.slf'.format(run,option), 'wb') f.write(struct.pack(">l",80)) strtemp='{0: >80}'.format('SERAFIN ') @@ -222,18 +226,128 @@ def writeConst(param,key,f): if option == 'Max': tempArr = [] for n in nodelayer[key]: - tempArr.append(tempR[key]) + tempArr.append(tempR[n]) tempVal1 = max(tempArr) elif isinstance(option, int): temppos = (option - 1 ) * 2 - temppos2 = min([temppos,len(nodelayer[key])]) - + temppos2 = min([temppos,len(nodelayer[key])-1]) tempVal1 = tempR[nodelayer[key][temppos2]] else: tempVal1 = tempR[key] #save the value in the selaphin file f.write(struct.pack(">f",tempVal1)) + +def writePCTL(c,pctl): + + f = open('{}_{}_{}.slf'.format(run,c,option), 'wb') + + f.write(struct.pack(">l",80)) + strtemp='{0: >80}'.format('SERAFIN ') + f.write(strtemp.encode('ascii')) + f.write(struct.pack(">l",80)) + + f.write(struct.pack(">l",8)) + f.write(struct.pack(">l",len(pctl))) + f.write(struct.pack(">l",0)) + f.write(struct.pack(">l",8)) + + for idx,p in enumerate(pctl): + f.write(struct.pack(">l",32)) + strtemp='{0: <32}'.format(str(idx)+ 'pct ' + str(p)) + f.write(strtemp.encode('ascii')) + f.write(struct.pack(">l",32)) + + + f.write(struct.pack(">l",40)) + f.write(struct.pack(">l",1)) + f.write(struct.pack(">l",0)) + f.write(struct.pack(">l",0)) + f.write(struct.pack(">l",0)) + f.write(struct.pack(">l",0)) + f.write(struct.pack(">l",0)) + f.write(struct.pack(">l",0)) + f.write(struct.pack(">l",0)) + f.write(struct.pack(">l",0)) + f.write(struct.pack(">l",1)) + f.write(struct.pack(">l",40)) + + f.write(struct.pack(">l",24)) + f.write(struct.pack(">l",startDate.year)) + f.write(struct.pack(">l",startDate.month)) + f.write(struct.pack(">l",startDate.day)) + f.write(struct.pack(">l",startDate.hour)) + f.write(struct.pack(">l",startDate.minute)) + f.write(struct.pack(">l",startDate.second)) + f.write(struct.pack(">l",24)) + + f.write(struct.pack(">l",16)) + f.write(struct.pack(">l",len(ElementList))) + f.write(struct.pack(">l",len(NodeList))) + f.write(struct.pack(">l",3)) + f.write(struct.pack(">l",1)) + f.write(struct.pack(">l",16)) + + + f.write(struct.pack(">l",len(ElementList)*3*4)) + for el in ElementList: + for nd in ElementDict[el]: + + f.write(struct.pack(">l",nodeOrdered[nd])) + f.write(struct.pack(">l",len(ElementList)*3*4)) + + f.write(struct.pack(">l",len(NodeList))) + for i in range(0,len(NodeList)): + f.write(struct.pack(">l",0)) + f.write(struct.pack(">l",len(NodeList))) + + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + f.write(struct.pack(">f",value[0])) + f.write(struct.pack(">l",len(NodeList)*4)) + + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + f.write(struct.pack(">f",value[1])) + f.write(struct.pack(">l",len(NodeList)*4)) + + columnName = [node2nodevalue[key] if key in node2nodevalue.keys() else key for key, value in nodeDict.items()] + global dfAll + dfAll = pd.DataFrame(columns=columnName) + cpt = 0 + + + currentDate = datetime(R.year,1,1) + timedelta(hours = R.time) + timeCurrentStep = currentDate - startDate + + f.write(struct.pack(">l",4)) + f.write(struct.pack(">f",timeCurrentStep.total_seconds())) + f.write(struct.pack(">l",4)) + + while R.next(): + + tempR = R.constit[c] + for key, value in nodeDict.items(): + if key in node2nodevalue.keys(): + key = node2nodevalue[key] + + tempArr = [] + for n in nodelayer[key]: + tempArr.append(tempR[n]) + tempVal1 = max(tempArr) + + dfAll.loc[cpt,key] = tempVal1 + cpt = cpt + 1 + + for p in pctl: + f.write(struct.pack(">l",len(NodeList)*4)) + + for key, value in nodeDict.items(): + if key in node2nodevalue.keys(): + key = node2nodevalue[key] + f.write(struct.pack(">f",dfAll[key].quantile(p))) + + f.write(struct.pack(">l",len(NodeList)*4)) def readNodeLayers(fname,Nnodes): fn = open(fname) @@ -314,24 +428,30 @@ for key, value in nodeDict.items(): nodeOrdered[key] = cpt cpt +=1 -Nnodes = max(NodeList) # # Open and Read First Step of the RMA File and Save a Serafin # In[72]: -for option in Options: - R=rma() - R.open(RMAfilename) - R.next() - startDate = datetime(R.year,1,1) + timedelta(hours = R.time) - if R.type==b'RMA11 ': - constName = [] - readNodeLayers(NodeLayersFile,Nnodes) - constName = ['X-VEL','Y-VEL','DEPTH','FREE SURFACE'] - for c in constNum: - constName.append(R.constit_name[c].decode("utf-8")) - print(c) - - RMA11toSerafin(option) +for run in RMAfilenames: + for option in Options: + RMAfilename = '{}/{}_WQ'.format(run,run) + R=rma() + R.open(RMAfilename) + R.next() + startDate = datetime(R.year,1,1) + timedelta(hours = R.time) + if R.type==b'RMA11 ': + constName = [] + readNodeLayers(NodeLayersFile,Nnodes) + constName = ['X-VEL','Y-VEL','DEPTH','FREE SURFACE'] + for c in constNum: + constName.append(R.constit_name[c].decode("utf-8")) + print(c) + + if option == 'Percentile': + df2 = pd.DataFrame() + for c in constNum: + writePCTL(c,percentiles) + else: + RMA11toSerafin(option) diff --git a/Seraphin files.pdf b/Seraphin files.pdf index 3845c17..09d4505 100644 Binary files a/Seraphin files.pdf and b/Seraphin files.pdf differ