1 # -*- coding: utf-8 -*-
3 Example of use case of interpolZ with the default values:
5 # --- case name in HYDRO
8 # --- med file 2D(x,y) of the case produced by SMESH
9 fichierMaillage = '/home/B27118/projets/LNHE/garonne/inondation1.med'
11 # --- dictionary [med group name] = region name
12 dicoGroupeRegion= dict(litMineur = 'inondation1_litMineur',
13 litMajeurDroite = 'inondation1_riveDroite',
14 litMajeurGauche = 'inondation1_riveGauche',
16 # --- Z interpolation on the bathymety/altimetry on the mesh nodes
17 interpolZ(nomCas, fichierMaillage, dicoGroupeRegion)
19 __revision__ = "V3.01"
20 # -----------------------------------------------------------------------------
22 # -----------------------------------------------------------------------------
27 theStudy = salome.myStudy
28 theStudyId = salome.myStudyId
31 import MEDLoader as ml
32 import MEDCoupling as mc
34 # -----------------------------------------------------------------------------
36 from med import medfile
37 from med import medmesh
38 from med import medfield
39 from med import medenum
41 # -----------------------------------------------------------------------------
45 class MyInterpolator( HYDROPy.HYDROData_IInterpolator ):
53 def GetAltitudeForPoint( self, theCoordX, theCoordY ):
57 alt_obj = HYDROPy.HYDROData_IInterpolator.GetAltitudeObject( self )
58 z = alt_obj.GetAltitudeForPoint( theCoordX, theCoordY ) # Custom calculation takes the base value and changes it to test
63 # -----------------------------------------------------------------------------
66 def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMethod=0, m3d=False, xyzFile=False, verbose=False):
68 interpolZ takes a 2D (x,y) mesh and calls the active instance of module HYDRO
69 to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node.
72 nomCas: Calculation Case Name in module HYDRO
73 fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case
74 dicoGroupeRegion: python dictionary giving the correspondance of mesh groups to HYDRO regions.
76 Value: region name in the HYDRO Case
77 zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
79 interpolMethod: integer value
80 0 = nearest point on bathymetry (default)
81 1 = linear interpolation
82 m3d: True/False to produce a 3D mesh. Default is False.
83 xyzFile: True/False to write an ascii file with xyz for every node. Default is False.
85 statz: statistique for z
87 Value: (minz, maxz, meanz, stdz, v05z, v95z)
89 return <fichierMaillage>F.med : med file with Z value in a field "BOTTOM"
90 Option: Z value also in Z coordinate if m3D is true
91 return <fichierMaillage>.xyz : text file with X, Y, Z values
100 ligne = "nomCas: %s" % nomCas
101 ligne += "\ninterpolMethod: %d" % interpolMethod
102 if (zUndef != None ):
103 ligne += "\nzUndef: %f" % zUndef
104 ligne += "\nm3d: %d" % m3d
107 doc = HYDROPy.HYDROData_Document.Document(theStudyId)
108 cas = doc.FindObjectByName(nomCas)
109 print ( "cas : ", cas)
110 custom_inter = MyInterpolator()
112 basename = fichierMaillage[:-4]
113 fichierFMaillage = basename + 'F.med'
115 print ("dicoGroupeRegion = ", dicoGroupeRegion)
116 ligne = "fichierMaillage = %s" % fichierMaillage
117 ligne += "\nfichierFMaillage = %s" % fichierFMaillage
119 fichierFonds = basename + '.xyz'
120 ligne += "\nfichierFonds = %s" % fichierFonds
125 meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
127 print (meshMEDFileRead)
129 # 2. Checks the names of the groups of faces
131 t_group_n = meshMEDFileRead.getGroupsNames()
132 dicoGroupeRegion_0 = dict()
134 for gr_face_name in dicoGroupeRegion:
135 saux = gr_face_name.strip()
136 if saux not in t_group_n:
137 message += "Group: '" + gr_face_name + "'\n"
140 dicoGroupeRegion_0[saux] = dicoGroupeRegion[gr_face_name]
142 ligne = "Number of problems: %d" % nb_pb
147 message += "This group does"
149 message += "These %d groups do" % nb_pb
150 message += " not belong to the mesh.\n"
151 message += "Please check the names of the group(s) of faces corresponding to each region of the HYDRO case.\n"
152 message += "Groups for this file:\n"
153 for group_n in t_group_n :
154 message += "'%s'\n" % group_n
158 # 3. Gets the information about the nodes
160 nbnodes = meshMEDFileRead.getNumberOfNodes()
162 ligne = "Number of nodes: %d" % nbnodes
165 coords = meshMEDFileRead.getCoords()
168 nb_comp = coords.getNumberOfComponents()
169 l_info = coords.getInfoOnComponents()
171 l_info_0=["X", "Y", "Z"]
172 for id_node in (0, 1, nbnodes-1):
173 ligne += "\nNode #%6d:" % id_node
174 for iaux in range(nb_comp):
178 saux = l_info_0[iaux]
179 ligne += " %s" % saux
180 ligne += "=%f" % coords[id_node, iaux]
183 # 4. Exploration of every group of faces
185 tb_aux = np.zeros(nbnodes, dtype=np.bool)
187 bathy = np.zeros(nbnodes, dtype=np.float)
190 for gr_face_name in dicoGroupeRegion_0:
192 # 4.1. Region connected to the group
194 nomreg = dicoGroupeRegion_0[gr_face_name]
195 ligne = "------- Region: '%s'" % nomreg
196 ligne += ", connected to group '%s'" % gr_face_name
198 region = doc.FindObjectByName(nomreg)
200 # 4.2. Mesh of the group
202 mesh_of_the_group = meshMEDFileRead.getGroup(0, gr_face_name, False)
203 nbr_cells = mesh_of_the_group.getNumberOfCells()
205 ligne = "\t. Number of cells: %d" % nbr_cells
208 # 4.3. Nodes of the meshes of the group
209 # Every node is flagged in tb_aux
212 for id_elem in range(nbr_cells):
213 l_nodes = mesh_of_the_group.getNodeIdsOfCell(id_elem)
215 for id_node in l_nodes:
216 tb_aux[id_node] = True
217 np_aux = tb_aux.nonzero()
220 ligne = "\t. Number of nodes for this group: %d" % len(np_aux[0])
222 #print ("np_aux:", np_aux)
224 # 4.4. Interpolation over the nodes of the meshes of the group
228 for nodeId in np_aux[0]:
229 vx.append(coords[nodeId, 0])
230 vy.append(coords[nodeId, 1])
234 vz = cas.GetAltitudesForPoints(vx, vy, region, interpolMethod)
241 v05z = np.percentile(vz, 05)
242 v95z = np.percentile(vz, 95)
245 ligne = ".. Minimum: %f" % minz
246 ligne += ", maximum: %f" % maxz
247 ligne += ", mean: %f\n" % meanz
248 ligne += ".. stdeviation: %f" % stdz
249 ligne += ", v05z: %f" % v05z
250 ligne += ", v95z: %f" % v95z
253 # 4.5. Storage of the z and of the statistics for this region
255 statz[gr_face_name] = (minz, maxz, meanz, stdz, v05z, v95z)
257 for iaux, nodeId in enumerate(np_aux[0]):
258 bathy[nodeId] = vz[iaux]
261 # During the interpolation, if no value is available over a node, a default value
262 # is set: -9999. It has no importance for the final computation, but if the field
263 # or the mesh is displayed, it makes huge gap. To prevent this artefact, a more
264 # convenient "undefined" value is set. This new undefined value is given by the user.
266 # zUndefThreshold: the default value is -9000. It is tied with the value -9999. given
267 # by the interpolation when no value is defined.
269 zUndefThreshold = -9000.
271 ligne = "zUndefThreshold: %f" % zUndefThreshold
274 #print "bathy :\n", bathy
275 np_aux_z = (bathy < zUndefThreshold).nonzero()
277 ligne = ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
280 for iaux in np_aux_z[0]:
283 # 6. Option : xyz file
288 ligne = ".. Ecriture du champ de bathymétrie sur le fichier :\n%s" % fichierFonds
291 with open(fichierFonds, "w") as fo :
292 for nodeId in range(nbnodes):
293 ligne = "%11.3f %11.3f %11.3f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
297 # 7.1. Transformation of the bathymetry as a double array
299 bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
300 bathy_dd.setInfoOnComponents(["Z [m]"])
302 # 7.2. If requested, modification of the z coordinate
305 coords3D = ml.DataArrayDouble.Meld([coords[:,0:2], bathy_dd])
306 coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
307 #print "coords3D =\n", coords3D
308 meshMEDFileRead.setCoords(coords3D)
310 # 7.3. Writes the mesh
317 ligne = ".. Ecriture du maillage" + saux
318 ligne += " sur le fichier :\n%s" % fichierFMaillage
321 meshMEDFileRead.write(fichierFMaillage, 2)
323 # 7.4. Writes the field
325 med_field_name = "BOTTOM"
327 ligne = ".. Ecriture du champ '%s'" % med_field_name
330 #print "bathy_dd =\n", bathy_dd
331 fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES)
332 fieldOnNodes.setName(med_field_name)
333 fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0))
334 fieldOnNodes.setArray(bathy_dd)
335 # Ces valeurs d'instants sont mises pour assurer la lecture par TELEMAC
337 # numero d'itération : 0
338 # pas de numero d'ordre (-1)
339 fieldOnNodes.setTime(0.0, 0, -1)
341 fMEDFile_ch_d = ml.MEDFileField1TS()
342 fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
343 fMEDFile_ch_d.write(fichierFMaillage, 0)