diff --git a/RMA2SERAPHIN20191118.py b/RMA2SERAPHIN20191118.py new file mode 100644 index 0000000..15f9c8e --- /dev/null +++ b/RMA2SERAPHIN20191118.py @@ -0,0 +1,373 @@ + +# coding: utf-8 + +# In[64]: + +import struct +import matplotlib.pyplot as plt +import math +from py_rmatools import rma + +plt.rcParams.update({'figure.max_open_warning': 0}) + + +# In[65]: + +meshFilename = 'HWQ010.rm1' +channelWidth = 100 +RMAfilename = 'HWQ277_2010_WQ' + +#If RMA11 +constNum = [1,2,3,5,6,7,8,10,11,13,14,15,18] + + +# In[66]: + +def isElementOneD(nodelist): + if len(nodelist) == 2: + return True + return False + +def isElementSquare(nodelist): + if len(nodelist) == 4: + return True + return False + +def square2Triangle(ElementNum): + nodelist = ElementDict[ElementNum] + if isElementSquare(nodelist): + ElementDict[ElementNum] = [nodelist[0], nodelist[1], nodelist[2]] + ElementList.append(max(ElementList) + 1) + ElementDict[ElementList[-1]]= [nodelist[0], nodelist[2], nodelist[3]] + +def oneD2triangle(ElementNum): + if isElementOneD(ElementDict[ElementNum]): + nAe = ElementDict[ElementNum][0] #nAe Node A existing + nBe = ElementDict[ElementNum][1] + + if not nAe in node1Dduplicate: node1Dduplicate[nAe] = [] + if not nBe in node1Dduplicate: node1Dduplicate[nBe] = [] + + xA = nodeDict[nAe][0] + xB = nodeDict[nBe][0] + yA = nodeDict[nAe][1] + yB = nodeDict[nBe][1] + + normalVec = [-(yB - yA),(xB - xA)] + dist = math.sqrt(normalVec[0]**2 + normalVec[1]**2) + normalVec[0] = normalVec[0] / dist + normalVec[1] = normalVec[1] / dist + xA2 = xA + channelWidth * normalVec[0] + xB2 = xB + channelWidth * normalVec[0] + yA2 = yA + channelWidth * normalVec[1] + yB2 = yB + channelWidth * normalVec[1] + + + nA = max(NodeList) + 1 + nB = max(NodeList) + 2 + + node1Dduplicate[nAe].append(nA) + node1Dduplicate[nBe].append(nB) + + node2nodevalue[nA] = nAe + node2nodevalue[nB] = nBe + + + NodeList.append(nA) + NodeList.append(nB) + nodeDict[nA] = [xA2, yA2, -1.01] + nodeDict[nB] = [xB2, yB2, -1.01] + + newEle = max(ElementList) + 1 + ElementList .append(newEle) + ElementDict[ElementNum] = [nAe, nA, nBe] + ElementDict[newEle] = [nA, nB, nBe] + +def RMA11toSerafin(): + f = open('{}.slf'.format(RMAfilename), 'wb') + + f.write(struct.pack(">l",80)) + str='{0: >80}'.format('SERAFIN ') + f.write(str.encode('ascii')) + f.write(struct.pack(">l",80)) + + f.write(struct.pack(">l",8)) + f.write(struct.pack(">l",len(constNum))) + f.write(struct.pack(">l",0)) + f.write(struct.pack(">l",8)) + + for c in constName: + f.write(struct.pack(">l",32)) + str='{0: <32}'.format(c) + f.write(str.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",R.year)) + f.write(struct.pack(">l",1)) + 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",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)) + + + + while R.next(): + #for i in range(3): + f.write(struct.pack(">l",4)) + f.write(struct.pack(">f",R.time * 3600)) + f.write(struct.pack(">l",4)) + + for c in constNum: + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + if key in node2nodevalue.keys(): + f.write(struct.pack(">f",R.constit[c][node2nodevalue[key]])) + else: + f.write(struct.pack(">f",R.constit[c][key])) + + f.write(struct.pack(">l",len(NodeList)*4)) + + #R.next() + + f.close() + + + +def RMA2toSerafin(): + f = open('{}.slf'.format(RMAfilename), 'wb') + + f.write(struct.pack(">l",80)) + str='{0: >80}'.format('SERAFIN ') + f.write(str.encode('ascii')) + f.write(struct.pack(">l",80)) + + f.write(struct.pack(">l",8)) + f.write(struct.pack(">l",len(constName))) + f.write(struct.pack(">l",0)) + f.write(struct.pack(">l",8)) + + for c in constName: + f.write(struct.pack(">l",32)) + str='{0: <32}'.format(c) + f.write(str.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",R.year)) + f.write(struct.pack(">l",1)) + 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",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)) + + + + while R.next(): + + f.write(struct.pack(">l",4)) + f.write(struct.pack(">f",R.time * 3600)) + f.write(struct.pack(">l",4)) + + + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + if key in node2nodevalue.keys(): + f.write(struct.pack(">f",R.xvel[node2nodevalue[key]])) + else: + f.write(struct.pack(">f",R.xvel[key])) + f.write(struct.pack(">l",len(NodeList)*4)) + + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + if key in node2nodevalue.keys(): + f.write(struct.pack(">f",R.yvel[node2nodevalue[key]])) + else: + f.write(struct.pack(">f",R.yvel[key])) + f.write(struct.pack(">l",len(NodeList)*4)) + + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + if key in node2nodevalue.keys(): + f.write(struct.pack(">f",R.depth[node2nodevalue[key]])) + else: + f.write(struct.pack(">f",R.depth[key])) + f.write(struct.pack(">l",len(NodeList)*4)) + + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + if key in node2nodevalue.keys(): + f.write(struct.pack(">f",R.elevation[node2nodevalue[key]])) + else: + f.write(struct.pack(">f",R.elevation[key])) + f.write(struct.pack(">l",len(NodeList)*4)) + + + R.next() + + f.close() + + +# In[67]: + +#Read mesh file and extract node (except mid node) and elements - plus convert 1D element to 2D for vizualisation +NodeList = [] +ElementList = [] +ElementDict = {} +nodeDict = {} +node1Dduplicate = {} #Original Number: List of Duplicates +node2nodevalue = {} #link between the node number and the node value to use + #(e.g. if node 10 is a 1D node: 10 is not duplicate so {1:1}, + #but node 2050 (duplicate of 10) (1D to 2D) the value of the duplicated + #node will be the same as the original so we might have {2050: 10}) + +with open(meshFilename) as f: + line = f.readline() + line = f.readline() + line = f.readline() + line = f.readline() + + cpt = 1 + while line and line != ' 9999\n': + temp = line.split() + ElementDict[int(temp[0])] = [int(temp[i]) for i in range(1,9,2) if temp[i] != '0' and int(temp[9]) < 100] + ElementList.append(int(temp[0])) + line = f.readline() + + for key, value in ElementDict.items(): + NodeList.extend(value) + + NodeList = list(set(NodeList)) + + line = f.readline() + while line and line != ' 9999\n': + temp = line.split() + if int(temp[0]) in NodeList: + nodeDict[int(temp[0])] = [float(temp[1]),float(temp[2]),float(temp[3])] + line = f.readline() + + +for e in ElementList: + oneD2triangle(e) + square2Triangle(e) + +for key in list(ElementDict): #Remove Special Element 902..... + if len(ElementDict[key]) != 3: + print(key, ElementDict[key]) + ElementDict.pop(key) + ElementList.remove(key) + +nodeOrdered = {} +cpt = 1 +for key, value in nodeDict.items(): + nodeOrdered[key] = cpt + cpt +=1 + + + +# # Open and Read First Step of the RMA File and Save a Serafin + +# In[72]: + +R=rma() +R.open(RMAfilename) +R.next() + +if R.type==b'RMA11 ': + constName = [] + for c in constNum: + constName.append(R.constit_name[c].decode("utf-8")) + RMA11toSerafin() + +if R.type==b'RMA2 ': + constName = ['X-VEL','Y-VEL','DEPTH','FREE SURFACE'] + RMA2toSerafin() + diff --git a/RMA2SERAPHIN_3DRMA11_20191118.py b/RMA2SERAPHIN_3DRMA11_20191118.py new file mode 100644 index 0000000..36d61be --- /dev/null +++ b/RMA2SERAPHIN_3DRMA11_20191118.py @@ -0,0 +1,738 @@ + +# coding: utf-8 + +# In[64]: + +import struct +import matplotlib.pyplot as plt +import math +from py_rmatools_v03 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}) + + +# In[65]: + +meshFilename = 'pok005.rm1' +channelWidth = 100 +RMAfilenames = ['POA034'] + +#If RMA11 +constNum = [1] + +#If RMA11 3D + +#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' (percentile only work for RMA11 FOR now) +Options = [8] +percentiles = [0.5,0.8,0.9,0.95,1] +NodeLayersFile = 'NodeLayers005.txt' #the node layer file was created by getting the nodes number from the RMA Outfile (just run RMA11 for a few steps) +Nnodes = 6674 #Number of nodes (including mid side node) - thisn umber also correspond to the last node number in the nodelayers file +# In[66]: + +def isElementOneD(nodelist1): + if len(nodelist1) == 3: + return True + return False + +def isElementSquare(nodelist1): + if len(nodelist1) == 8: + return True + return False + +def isElementTriangle(nodelist1): + if len(nodelist1) == 6: + return True + return False + +def square2Triangle(ElementNum): + nodelist1 = ElementDict[ElementNum] + if isElementSquare(nodelist1): + ElementDict[ElementNum] = [nodelist1[0], nodelist1[1], nodelist1[7]] + ElementList.append(max(ElementList) + 1) + ElementDict[ElementList[-1]]= [nodelist1[1], nodelist1[2], nodelist1[3]] + ElementList.append(max(ElementList) + 1) + ElementDict[ElementList[-1]]= [nodelist1[1], nodelist1[3], nodelist1[7]] + ElementList.append(max(ElementList) + 1) + ElementDict[ElementList[-1]]= [nodelist1[3], nodelist1[4], nodelist1[5]] + ElementList.append(max(ElementList) + 1) + ElementDict[ElementList[-1]]= [nodelist1[3], nodelist1[5], nodelist1[7]] + ElementList.append(max(ElementList) + 1) + ElementDict[ElementList[-1]]= [nodelist1[5], nodelist1[6], nodelist1[7]] + +def triangle2Triangles(ElementNum): + nodelist1 = ElementDict[ElementNum] + if isElementTriangle(nodelist1): + ElementDict[ElementNum] = [nodelist1[0], nodelist1[1], nodelist1[5]] + ElementList.append(max(ElementList) + 1) + ElementDict[ElementList[-1]]= [nodelist1[1], nodelist1[2], nodelist1[3]] + ElementList.append(max(ElementList) + 1) + ElementDict[ElementList[-1]]= [nodelist1[1], nodelist1[3], nodelist1[5]] + ElementList.append(max(ElementList) + 1) + ElementDict[ElementList[-1]]= [nodelist1[3], nodelist1[4], nodelist1[5]] + +def oneD2triangle(ElementNum): + if isElementOneD(ElementDict[ElementNum]): + nAe = ElementDict[ElementNum][0] #nAe Node A existing + nBe = ElementDict[ElementNum][1] + + if not nAe in node1Dduplicate: node1Dduplicate[nAe] = [] + if not nBe in node1Dduplicate: node1Dduplicate[nBe] = [] + + xA = nodeDict[nAe][0] + xB = nodeDict[nBe][0] + yA = nodeDict[nAe][1] + yB = nodeDict[nBe][1] + + normalVec = [-(yB - yA),(xB - xA)] + dist = math.sqrt(normalVec[0]**2 + normalVec[1]**2) + normalVec[0] = normalVec[0] / dist + normalVec[1] = normalVec[1] / dist + xA2 = xA + channelWidth * normalVec[0] + xB2 = xB + channelWidth * normalVec[0] + yA2 = yA + channelWidth * normalVec[1] + yB2 = yB + channelWidth * normalVec[1] + + + nA = max(NodeList) + 1 + nB = max(NodeList) + 2 + + node1Dduplicate[nAe].append(nA) + node1Dduplicate[nBe].append(nB) + + node2nodevalue[nA] = nAe + node2nodevalue[nB] = nBe + + + NodeList.append(nA) + NodeList.append(nB) + nodeDict[nA] = [xA2, yA2, -1.01] + nodeDict[nB] = [xB2, yB2, -1.01] + + newEle = max(ElementList) + 1 + ElementList .append(newEle) + ElementDict[ElementNum] = [nAe, nA, nBe] + ElementDict[newEle] = [nA, nB, nBe] + +def RMA2toSerafin(): + f = open('{}.slf'.format(RMAfilename), 'wb') + + f.write(struct.pack(">l",80)) + str='{0: >80}'.format('SERAFIN ') + f.write(str.encode('ascii')) + f.write(struct.pack(">l",80)) + + f.write(struct.pack(">l",8)) + f.write(struct.pack(">l",len(constName))) + f.write(struct.pack(">l",0)) + f.write(struct.pack(">l",8)) + + for c in constName: + f.write(struct.pack(">l",32)) + str='{0: <32}'.format(c) + f.write(str.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",R.year)) + f.write(struct.pack(">l",1)) + 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",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)) + + + + while R.next(): + + f.write(struct.pack(">l",4)) + f.write(struct.pack(">f",R.time * 3600)) + f.write(struct.pack(">l",4)) + + + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + if key in node2nodevalue.keys(): + f.write(struct.pack(">f",R.xvel[node2nodevalue[key]])) + else: + f.write(struct.pack(">f",R.xvel[key])) + f.write(struct.pack(">l",len(NodeList)*4)) + + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + if key in node2nodevalue.keys(): + f.write(struct.pack(">f",R.yvel[node2nodevalue[key]])) + else: + f.write(struct.pack(">f",R.yvel[key])) + f.write(struct.pack(">l",len(NodeList)*4)) + + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + if key in node2nodevalue.keys(): + f.write(struct.pack(">f",R.depth[node2nodevalue[key]])) + else: + f.write(struct.pack(">f",R.depth[key])) + f.write(struct.pack(">l",len(NodeList)*4)) + + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + if key in node2nodevalue.keys(): + f.write(struct.pack(">f",R.elevation[node2nodevalue[key]])) + else: + f.write(struct.pack(">f",R.elevation[key])) + f.write(struct.pack(">l",len(NodeList)*4)) + + + f.close() + +def RMA11toSerafin(option=1): + f = open('{}_{}.slf'.format(run,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(constName))) + f.write(struct.pack(">l",0)) + f.write(struct.pack(">l",8)) + + for idx,c in enumerate(constName): + f.write(struct.pack(">l",32)) + strtemp='{0: <32}'.format(str(idx)+c) + 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)) + + + + while R.next(): + #for i in range(3): + 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)) + + + + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + writeConst('X-VEL',key,f) + + f.write(struct.pack(">l",len(NodeList)*4)) + + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + writeConst('Y-VEL',key,f) + f.write(struct.pack(">l",len(NodeList)*4)) + + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + writeConst('DEPTH',key,f) + f.write(struct.pack(">l",len(NodeList)*4)) + + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + writeConst('FREE SURFACE',key,f) + f.write(struct.pack(">l",len(NodeList)*4)) + + for c in constNum: + f.write(struct.pack(">l",len(NodeList)*4)) + for key, value in nodeDict.items(): + writeConst(c,key,f) + f.write(struct.pack(">l",len(NodeList)*4)) + f.close() + + +def RMA10toSerafin(option=1): + f = open('{}_{}.slf'.format(run,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(constName))) + f.write(struct.pack(">l",0)) + f.write(struct.pack(">l",8)) + + for idx,c in enumerate(constName): + f.write(struct.pack(">l",32)) + strtemp='{0: <32}'.format(str(idx)+c) + 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)) + + + + while R.next(): + #for i in range(3): + 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)) + + + +# f.write(struct.pack(">l",len(NodeList)*4)) +# for key, value in nodeDict.items(): +# writeConst('X-VEL',key,f) +# f.write(struct.pack(">l",len(NodeList)*4)) + + + + for c in constName: + f.write(struct.pack(">l",len(NodeList)*4)) + data = [None] * len(nodeDict) + idx = 0 + for key, value in nodeDict.items(): + data[idx] = writeConst(c,key) + idx += 1 + f.write(b''.join(data)) + f.write(struct.pack(">l",len(NodeList)*4)) + f.close() + +def writeConst(param,key): +#get the surface node number +# ['X-VEL','Y-VEL','DEPTH','FREE SURFACE','Constituent 1'......] + + if param == 'X-VEL': + tempR = R.xvel + elif param == 'Y-VEL': + tempR = R.yvel + elif param == 'Z-VEL': + tempR = R.zvel + elif param == 'DEPTH': + tempR = R.depth + elif param == 'FREE SURFACE': + tempR = R.elevation + elif param == 'TEMPERATURE': + tempR = R.temperature + elif param == 'SALINITY': + tempR = R.salinity + elif param == 'SSED': + tempR = R.sussed + else: + tempR = R.constit[c] + + if key in node2nodevalue.keys(): + key = node2nodevalue[key] + #based on the surface node number and option selected calculate the value + if option == 'Max': + tempArr = [] + for n in nodelayer[key]: + tempArr.append(tempR[n]) + tempVal1 = max(tempArr) + elif isinstance(option, int): + temppos = (option - 1 ) * 2 + 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)) + return 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 = {} + + for key, value in nodeDict.items(): + if key in node2nodevalue.keys(): + key = node2nodevalue[key] + + dfAll[key] = np.array([]) + + 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[key] = np.append(dfAll[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",np.quantile(dfAll[key],p))) + + f.write(struct.pack(">l",len(NodeList)*4)) + +def readNodeLayers(fname,Nnodes): + fn = open(fname) + lines = fn.readlines() + + surfaceNode = int(lines[0]) + nodelayer[surfaceNode] = [] + for l in lines: + l = int(l) + if (l != surfaceNode) & (l < Nnodes + 1): + surfaceNode = l + nodelayer[surfaceNode] = [] + nodelayer[surfaceNode].append(l) + +# In[67]: + +#Read mesh file and extract node (except mid node) and elements - plus convert 1D element to 2D for vizualisation +nodelayer = {} +NodeList = [] +ElementList = [] +ElementDict = {} +nodeDict = {} +node1Dduplicate = {} #Original Number: List of Duplicates +node2nodevalue = {} #link between the node number and the node value to use + #(e.g. if node 10 is a 1D node: 10 is not duplicate so {1:1}, + #but node 2050 (duplicate of 10) (1D to 2D) the value of the duplicated + #node will be the same as the original so we might have {2050: 10}) + +with open(meshFilename) as f: + line = f.readline() + line = f.readline() + line = f.readline() + line = f.readline() + + cpt = 1 + while line and line != ' 9999\n': + + #temp = line.split() + + temp = re.findall('.....',line) + ElementDict[int(temp[0])] = [int(temp[i]) for i in range(1,9) if int(temp[i]) != 0 and int(temp[9]) < 100] + ElementList.append(int(temp[0])) + line = f.readline() + + for key, value in ElementDict.items(): + NodeList.extend(value) + + NodeList = list(set(NodeList)) + + line = f.readline() + while line and line != ' 9999\n': + + formatFix = (10,16,20,14,10,10) + temp = [] + for form in formatFix: + temp.append(line[:form]) + line = line[form:] + + #line.split() + if int(temp[0]) in NodeList: + nodeDict[int(temp[0])] = [float(temp[1]),float(temp[2]),float(temp[3])] + line = f.readline() + +print('Processing triangles') +ElementList2 = ElementList[:] +for e in ElementList2: + print(e) + if isElementTriangle(ElementDict[e]): + triangle2Triangles(e) + elif isElementOneD(ElementDict[e]): + oneD2triangle(e) + else: + square2Triangle(e) + +for key in list(ElementDict): #Remove Special Element 902..... + if len(ElementDict[key]) != 3: + print(key, ElementDict[key]) + ElementDict.pop(key) + ElementList.remove(key) + +nodeOrdered = {} +cpt = 1 +for key, value in nodeDict.items(): + nodeOrdered[key] = cpt + cpt +=1 + +# # Open and Read First Step of the RMA File and Save a Serafin + +# In[72]: +for run in RMAfilenames: + for option in Options: + RMAfilename = '{}'.format(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': + for c in constNum: + writePCTL(c,percentiles) + else: + RMA11toSerafin(option) + + if R.type=='RMA10 ': + constName = [] + readNodeLayers(NodeLayersFile,Nnodes) + constName = ['X-VEL','Y-VEL','Z-VEL','FREE SURFACE','SALINITY','TEMPERATURE'] + + + if option == 'Percentile': + for c in constNum: + writePCTL(c,percentiles) + else: + RMA10toSerafin(option) + + if R.type==b'RMA2 ': + constName = ['X-VEL','Y-VEL','DEPTH','FREE SURFACE'] + RMA2toSerafin() + + + diff --git a/py_rmatools.py b/py_rmatools.py index 0a350fd..6fc255a 100644 --- a/py_rmatools.py +++ b/py_rmatools.py @@ -136,9 +136,9 @@ class rma: np=int(a[2]) self.year=a[3] if ((nqal-5)!=(self.num_constits)): - print ("Warning - NQAL-5 (%d) on this timestep does not match header (%d)" % (nqal-5,self.num_constits)) + print ("Warning - NQAL-5 {} on this timestep {} does not match header {}".format(nqal-5,self.time,self.num_constits)) if (np!=self.num_nodes): - print ("Warning - NP (%d) on this timestep does not match header (%d)" % (np,self.num_nodes)) + print ("Warning - NP {} on this timestep {} does not match header {}".format(np,self.time,self.num_nodes)) b=unpack('%df' % nqal*np, self.file.read(4*nqal*np)) i=1 while i<=np: diff --git a/py_rmatools_v03.py b/py_rmatools_v03.py new file mode 100644 index 0000000..7e2af6e --- /dev/null +++ b/py_rmatools_v03.py @@ -0,0 +1,212 @@ +from struct import unpack +#from pylab import * +import sys +from array import array + + +# Updated 10/11/15 BMM : Added mesh. +# Updated 10/11/15 MD : Added dictionnary initialisation in mesh. +# Update 8/11/2019 MD : Improve performance (used array instead of unpack) + +class rma: + def __init__(self): + self.file = 'not opened yet' + + # May want to do some other things here later + + def open(self, filename): + self.file = open(('%s.rma' % (filename)), 'rb') + self.header = self.file.read(1000).decode("utf-8") + self.type = self.header[0:10] + self.title = self.header[100:172] + self.geometry = self.header[200:300] + self.num_nodes = int(self.header[40:50]) + self.num_elements = int(self.header[50:60]) + if self.type == 'RMA11 ': + self.num_constits = int(self.header[60:70]) + if len(self.header[80:90].strip()) == 0: + self.num_sedlayers = 0 + else: + self.num_sedlayers = int(self.header[80:90]) + self.constit_name = [] + self.constit_name.append("NULL") + i = 1 + print(self.num_constits) + print(self.header[300:1000]) + while i <= self.num_constits: + # print self.header[300:400] + self.constit_name.append(self.header[300 + (i - 1) * 8:308 + (i - 1) * 8]) + if self.header[300 + (i - 1) * 8:308 + (i - 1) * 8] == " SSED ": + self.constit_name.append(" BSHEAR") + self.constit_name.append("BedThick") + j = 1 + print(self.num_sedlayers) + while j <= self.num_sedlayers: + self.constit_name.append(" L%dThick" % j) + j = j + 1 + i = i + 1 + # print i + + # print self.num_constits + # print i + # if self.num_sedlayers > 0: + # self.num_constits=self.num_constits+2+self.num_sedlayers + + # Header format in RMA2 + + # HEADER(1:10) ='RMA2 ' + # HEADER(11:20)=DATEC = Date of run + # HEADER(21:30)=TIMEC = Time of run + # HEADER(31:40)=ZONEC = Time zone + # HEADER(41:50)=NP = Number of nodes + # HEADER(51:60)=NE = Number of elements + # HEADER(101:172)=TITLE = Title line + # HEADER(201:300)=FNAME = File name of binary geometry + + # Header format in RMA11 + # HEADER(1:10) - Model + # HEADER(11:20) - DATEC of model runtime + # HEADER(21:30) - TIMEC of model runtime + # HEADER(31:40) - ZONEC of model runtime + # HEADER(41:90) - '(5I10)' NP,NE,NQAL,IGS,LGS or '(5I10)' NP,NE,NQAL,ISEDS,NLAYT + # HEADER(101:172) - TITLE + # HEADER(301:540) - '(30A8)' CLABL + # HEADER(201:248) - FNAMD directory of the geometry + # HEADER(251:298) - FNAM2 filename of the geometry + + def mesh(self, filename): + self.meshfile = open(('%s.rm1' % (filename)), 'r') + l = self.meshfile.readline() + + while l[0:5] != ' 9999' and len(l) != 5: + l = self.meshfile.readline() + l = self.meshfile.readline() + self.node_e = {} + self.node_n = {} + self.node_z = {} + while l[0:10] != ' 9999' and len(l) != 10: + ll = l.split() + self.node_e[ll[0]] = ll[1] + self.node_n[ll[0]] = ll[2] + self.node_z[ll[0]] = ll[3] + l = self.meshfile.readline + + def next(self): + # This reads the next timestep and populates the variables. + # Reset all of the variables + self.time = 0.0 + self.year = 0 + self.xvel = [] + self.yvel = [] + self.zvel = [] + self.depth = [] + self.elevation = [] + self.temperature = [] + self.salinity = [] + self.sussed = [] + + # Add in an entry to fill position 0. Important since RMA files start numbering nodes from 1. + self.xvel.append(-999) + self.yvel.append(-999) + self.zvel.append(-999) + self.depth.append(-999) + self.elevation.append(-999) + self.temperature.append(-999) + self.salinity.append(-999) + self.sussed.append(-999) + + if self.type == 'RMA2 ': + # WRITE(IRMAFM) TETT,NP,IYRR,((VEL(J,K),J=1,3),K=1,NP),(WSEL(J),J=1,NP),(VDOT(3,K),K=1,NP) + t = self.file.read(12) + if t: + a = unpack('fii', t) + self.time = a[0] + np = int(a[1]) + self.year = a[2] + if (np != self.num_nodes): + print("Warning - NP (%d) on this timestep does not match header (%d)" % (np, self.num_nodes)) + #b = unpack('%df' % 5 * np, self.file.read(20 * np)) + b = array('f') + b.fromfile(self.file, 5 * np) + +# b = np.fromfile(self.file, dtype=float,count = 5 * np) + + tempA = [(x - 1) * 3 for x in range(1, np + 1)] + tempB = [np * 3 + (x - 1) for x in range(1, np + 1)] + self.xvel.extend([b[i] for i in tempA]) + self.yvel.extend([b[i + 1] for i in tempA]) + self.depth.extend([b[i + 2] for i in tempA]) + self.elevation.extend([b[i] for i in tempB]) + + if self.type == 'RMA11 ': + self.constit = [] + c = 0 + # Start at zero to leave an empty array space. Constits numbered 1 .. num_constits + while c <= self.num_constits: + self.constit.append([]) + # Put a null into the 0 position so that RMA node numbers are in the same numbered array position. + self.constit[c].append(-999) + c = c + 1 + # READ(file1,END=100) TETT1,NQAL,NP,IYRR, ((VEL(K,J),J=1,NP),K=1,3), (wd(j),j=1,np), (wsel(j),j=1,np), ((TCON1(K,J),J=1,NP),K=1,NQAL-5) + t = self.file.read(16) + if t: + a = unpack('fiii', t) + self.time = a[0] + nqal = int(a[1]) + np = int(a[2]) + self.year = a[3] + if ((nqal - 5) != (self.num_constits)): + print("Warning - NQAL-5 (%d) on this timestep does not match header (%d)" % ( + nqal - 5, self.num_constits)) + if (np != self.num_nodes): + print("Warning - NP (%d) on this timestep does not match header (%d)" % (np, self.num_nodes)) + + b = array('f') + b.fromfile(self.file, nqal * np) + #b = unpack('%df' % nqal * np, self.file.read(4 * nqal * np)) #4 because single precision float type have 4 bytes + + tempA = [(x - 1) * 3 for x in range(1, np + 1)] + tempB = [np * 3 + (x - 1) for x in range(1, np + 1)] + tempC = [np * 4 + (x - 1) for x in range(1, np + 1)] + self.xvel.extend([b[i] for i in tempA]) + self.yvel.extend([b[i + 1] for i in tempA]) + self.zvel.extend([b[i + 2] for i in tempA]) + self.depth.extend([b[i] for i in tempB]) + self.elevation.extend([b[i] for i in tempC]) + for c in range(1, self.num_constits + 1): + (self.constit[c].extend([b[np * ((c - 1) + 5) + (x - 1)] for x in range(1, np + 1)])) + if self.type == 'RMA10 ': + # WRITE(IRMAFM) TETT,NP,NDF,NE,IYRR,((VSING(K,J),K=1,NDF),VVEL(J),WSLL(J),J=1,NP),(DFCT(J),J=1,NE),(VSING(7,J),J=1,NP) + # WRITE(IRMAFM) TETT,NP,IYRR,((VEL(J,K),J=1,3),K=1,NP),(WSEL(J),J=1,NP),(VDOT(3,K),K=1,NP) + t = self.file.read(20) + if t: + a = unpack('fiiii', t) + self.time = a[0] + np = a[1] + ndf = 6 + ne = a[3] + self.year = a[4] + if (np != self.num_nodes): + print("Warning - NP1 (%d) on this timestep does not match header (%d)" % (np, self.num_nodes)) + tempRead = np * (3 + ndf) + ne + b = array('f') + b.fromfile(self.file, tempRead) +# b = numpy.fromfile(self.file.read(4 * tempRead), dtype=float) + #b = unpack('%df' % tempRead, self.file.read(4 * tempRead)) + i = 1 + while i <= np: + self.xvel.append(b[0 + (i - 1) * 8]) + self.yvel.append(b[1 + (i - 1) * 8]) + self.depth.append(b[2 + (i - 1) * 8]) + self.salinity.append(b[3 + (i - 1) * 8]) + self.temperature.append(b[4 + (i - 1) * 8]) + self.sussed.append(b[5 + (i - 1) * 8]) + self.zvel.append(b[6 + (i - 1) * 8]) + self.elevation.append(b[7 + (i - 1) * 8]) + i = i + 1 + + if len(self.xvel) == 1: + # Note that this is a 1 now as we've filled the zero array position + return False + else: + return True