Upgrade 3DRMA11 - add percentile - fix bugs

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:
+for run in RMAfilenames:
-    R=rma()
+    for option in Options:
-    R.open(RMAfilename)
+        RMAfilename = '{}/{}_WQ'.format(run,run)
-    R.next()
+        R=rma()
-    startDate = datetime(R.year,1,1) + timedelta(hours = R.time)
+        R.open(RMAfilename)
-    if R.type==b'RMA11     ':
+        R.next()
-        constName = []
+        startDate = datetime(R.year,1,1) + timedelta(hours = R.time)
-        readNodeLayers(NodeLayersFile,Nnodes)
+        if R.type==b'RMA11     ':
-        constName = ['X-VEL','Y-VEL','DEPTH','FREE SURFACE'] 
+            constName = []
-        for c in constNum:
+            readNodeLayers(NodeLayersFile,Nnodes)
-            constName.append(R.constit_name[c].decode("utf-8"))
+            constName = ['X-VEL','Y-VEL','DEPTH','FREE SURFACE'] 
-            print(c)
+            for c in constNum:
-
+                constName.append(R.constit_name[c].decode("utf-8"))
-        RMA11toSerafin(option)
+                print(c)
+        
+            if option == 'Percentile':
+                df2 = pd.DataFrame()
+                for c in constNum:
+                    writePCTL(c,percentiles)
+            else:
+                RMA11toSerafin(option)