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 # -----------------------------------------------------------------------------
29 import MEDLoader as ml
30 import medcoupling as mc
32 # -----------------------------------------------------------------------------
34 from med import medfile
35 from med import medmesh
36 from med import medfield
37 from med import medenum
39 # -----------------------------------------------------------------------------
43 class MyInterpolator( HYDROPy.HYDROData_IInterpolator ):
51 def GetAltitudeForPoint( self, theCoordX, theCoordY ):
55 alt_obj = HYDROPy.HYDROData_IInterpolator.GetAltitudeObject( self )
56 z = alt_obj.GetAltitudeForPoint( theCoordX, theCoordY ) # Custom calculation takes the base value and changes it to test
61 # -----------------------------------------------------------------------------
64 def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., regions_interp_method=None, m3d=True, xyzFile=False, verbose=False):
66 interpolZ takes a 2D (x,y) mesh and calls the active instance of module HYDRO
67 to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node.
70 nomCas: Calculation Case Name in module HYDRO
71 fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case
72 dicoGroupeRegion: python dictionary giving the correspondance of mesh groups to HYDRO regions.
74 Value: region name in the HYDRO Case
75 zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
77 interpolMethod: integer value or dict
79 0 = nearest point on bathymetry
80 1 = linear interpolation
81 if dict : key is a region name, value is a type of interpolation (0 or 1)
82 if None : default type used (0)
83 m3d: True/False to produce a 3D mesh. Default is False.
84 xyzFile: True/False to write an ascii file with xyz for every node. Default is False.
86 statz: statistique for z
88 Value: (minz, maxz, meanz, stdz, v05z, v95z)
90 return <fichierMaillage>F.med : med file with Z value in a field "BOTTOM"
91 Option: Z value also in Z coordinate if m3D is true
92 return <fichierMaillage>.xyz : text file with X, Y, Z values
101 ligne = "nomCas: %s" % nomCas
102 ligne += "\ninterpolMethods: %s" % regions_interp_method
103 if (zUndef != None ):
104 ligne += "\nzUndef: %f" % zUndef
105 ligne += "\nm3d: %d" % m3d
108 doc = HYDROPy.HYDROData_Document.Document()
109 cas = doc.FindObjectByName(nomCas)
110 print ( "cas : ", cas)
111 custom_inter = MyInterpolator()
113 basename = fichierMaillage[:-4]
114 fichierFMaillage = basename + 'F.med'
116 print ("dicoGroupeRegion = ", dicoGroupeRegion)
117 ligne = "fichierMaillage = %s" % fichierMaillage
118 ligne += "\nfichierFMaillage = %s" % fichierFMaillage
120 fichierFonds = basename + '.xyz'
121 ligne += "\nfichierFonds = %s" % fichierFonds
126 meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
128 print (meshMEDFileRead)
130 # 2. Checks the names of the groups of faces
132 t_group_n = meshMEDFileRead.getGroupsNames()
133 dicoGroupeRegion_0 = dict()
135 for gr_face_name in dicoGroupeRegion:
136 saux = gr_face_name.strip()
137 if saux not in t_group_n:
138 message += "Group: '" + gr_face_name + "'\n"
141 dicoGroupeRegion_0[saux] = dicoGroupeRegion[gr_face_name]
143 ligne = "Number of problems: %d" % nb_pb
148 message += "This group does"
150 message += "These %d groups do" % nb_pb
151 message += " not belong to the mesh.\n"
152 message += "Please check the names of the group(s) of faces corresponding to each region of the HYDRO case.\n"
153 message += "Groups for this file:\n"
154 for group_n in t_group_n :
155 message += "'%s'\n" % group_n
159 # 3. Gets the information about the nodes
161 nbnodes = meshMEDFileRead.getNumberOfNodes()
163 ligne = "Number of nodes: %d" % nbnodes
166 coords = meshMEDFileRead.getCoords()
169 nb_comp = coords.getNumberOfComponents()
170 l_info = coords.getInfoOnComponents()
172 l_info_0=["X", "Y", "Z"]
173 for id_node in (0, 1, nbnodes-1):
174 ligne += "\nNode #%6d:" % id_node
175 for iaux in range(nb_comp):
179 saux = l_info_0[iaux]
180 ligne += " %s" % saux
181 ligne += "=%f" % coords[id_node, iaux]
184 # 4. Exploration of every group of faces
186 tb_aux = np.zeros(nbnodes, dtype=np.bool)
188 bathy = np.zeros(nbnodes, dtype=np.float)
191 for gr_face_name in dicoGroupeRegion_0:
193 # 4.1. Region connected to the group
195 nomreg = dicoGroupeRegion_0[gr_face_name]
196 ligne = "------- Region: '%s'" % nomreg
197 ligne += ", connected to group '%s'" % gr_face_name
199 region = doc.FindObjectByName(nomreg)
201 # 4.2. Mesh of the group
203 mesh_of_the_group = meshMEDFileRead.getGroup(0, gr_face_name, False)
204 nbr_cells = mesh_of_the_group.getNumberOfCells()
206 ligne = "\t. Number of cells: %d" % nbr_cells
209 # 4.3. Nodes of the meshes of the group
210 # Every node is flagged in tb_aux
213 for id_elem in range(nbr_cells):
214 l_nodes = mesh_of_the_group.getNodeIdsOfCell(id_elem)
216 for id_node in l_nodes:
217 tb_aux[id_node] = True
218 np_aux = tb_aux.nonzero()
221 ligne = "\t. Number of nodes for this group: %d" % len(np_aux[0])
223 #print ("np_aux:", np_aux)
225 # 4.4. Interpolation over the nodes of the meshes of the group
229 for nid in np_aux[0]:
231 vx.append(coords[nodeId, 0])
232 vy.append(coords[nodeId, 1])
237 if regions_interp_method is not None:
238 if isinstance(regions_interp_method, dict) and nomreg in list(regions_interp_method.keys()):
239 interpolMethod = int(regions_interp_method[nomreg])
240 elif isinstance(regions_interp_method, int):
241 interpolMethod = regions_interp_method
243 #print ('interp', interpolMethod )
244 vz = cas.GetAltitudesForPoints(vx, vy, region, interpolMethod)
251 v05z = np.percentile(vz, 0o5)
252 v95z = np.percentile(vz, 95)
255 ligne = ".. Minimum: %f" % minz
256 ligne += ", maximum: %f" % maxz
257 ligne += ", mean: %f\n" % meanz
258 ligne += ".. stdeviation: %f" % stdz
259 ligne += ", v05z: %f" % v05z
260 ligne += ", v95z: %f" % v95z
263 # 4.5. Storage of the z and of the statistics for this region
265 statz[gr_face_name] = (minz, maxz, meanz, stdz, v05z, v95z)
267 for iaux, nodeId in enumerate(np_aux[0]):
268 bathy[nodeId] = vz[iaux]
271 # During the interpolation, if no value is available over a node, a default value
272 # is set: -9999. It has no importance for the final computation, but if the field
273 # or the mesh is displayed, it makes huge gap. To prevent this artefact, a more
274 # convenient "undefined" value is set. This new undefined value is given by the user.
276 # zUndefThreshold: the default value is -9000. It is tied with the value -9999. given
277 # by the interpolation when no value is defined.
279 zUndefThreshold = -9000.
281 ligne = "zUndefThreshold: %f" % zUndefThreshold
284 #print "bathy :\n", bathy
285 np_aux_z = (bathy < zUndefThreshold).nonzero()
287 ligne = ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
290 for iaux in np_aux_z[0]:
293 # 6. Option : xyz file
298 ligne = ".. Ecriture du champ de bathymétrie sur le fichier :\n%s" % fichierFonds
301 with open(fichierFonds, "w") as fo :
302 for nodeId in range(nbnodes):
303 ligne = "%11.3f %11.3f %11.3f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
307 # 7.1. Transformation of the bathymetry as a double array
309 bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
310 bathy_dd.setInfoOnComponents(["Z [m]"])
312 # 7.2. If requested, modification of the z coordinate
315 coords3D = ml.DataArrayDouble.Meld([coords[:,0:2], bathy_dd])
316 coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
317 #print "coords3D =\n", coords3D
318 meshMEDFileRead.setCoords(coords3D)
320 # 7.3. Writes the mesh
327 ligne = ".. Ecriture du maillage" + saux
328 ligne += " sur le fichier :\n%s" % fichierFMaillage
331 meshMEDFileRead.write(fichierFMaillage, 2)
333 # 7.4. Writes the field
335 med_field_name = "BOTTOM"
337 ligne = ".. Ecriture du champ '%s'" % med_field_name
340 #print "bathy_dd =\n", bathy_dd
341 fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES)
342 fieldOnNodes.setName(med_field_name)
343 fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0))
344 fieldOnNodes.setArray(bathy_dd)
345 # Ces valeurs d'instants sont mises pour assurer la lecture par TELEMAC
347 # numero d'itération : 0
348 # pas de numero d'ordre (-1)
349 fieldOnNodes.setTime(0.0, 0, -1)
351 fMEDFile_ch_d = ml.MEDFileField1TS()
352 fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
353 fMEDFile_ch_d.write(fichierFMaillage, 0)
363 def interpolZ_B(bathyName, fichierMaillage, gr_face_name, zUndef=90., interp_method=0, m3d=True, xyzFile=False, verbose=False):
365 interpolZ_B takes a 2D (x,y) mesh and calls the active instance of module HYDRO
366 to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node.
369 bathyName: Bathymetry Name in module HYDRO
370 fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case
371 gr_face_name: face group name
372 zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
374 interp_method: interpolation method
375 0 = nearest point on bathymetry
376 1 = linear interpolation
377 m3d: True/False to produce a 3D mesh. Default is True.
378 xyzFile: True/False to write an ascii file with xyz for every node. Default is False.
380 if OK : statz_gr_face_name: statistic for z: (minz, maxz, meanz, stdz, v05z, v95z)
383 return <fichierMaillage>F.med : med file with Z value in a field "BOTTOM"
384 Option: Z value also in Z coordinate if m3D is true
385 return <fichierMaillage>.xyz : text file with X, Y, Z values
387 doc = HYDROPy.HYDROData_Document.Document()
388 bathy_obj = doc.FindObjectByName(bathyName)
389 print( "bathy : ", bathy_obj)
390 if bathy_obj is None:
391 print ( "bathy is None")
394 basename = fichierMaillage[:-4]
395 fichierFMaillage = basename + 'F.med'
397 ligne = "fichierMaillage = %s" % fichierMaillage
398 ligne += "\nfichierFMaillage = %s" % fichierFMaillage
401 fichierFonds = basename + '.xyz'
402 ligne += "\nfichierFonds = %s" % fichierFonds
407 meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
409 print (meshMEDFileRead)
411 # 2. Checks the names of the groups of faces
413 t_group_n = meshMEDFileRead.getGroupsNames()
414 gr_face_name_tr = gr_face_name.strip()
415 if gr_face_name_tr not in t_group_n:
416 print("Group not found")
419 # 3. Gets the information about the nodes
421 nbnodes = meshMEDFileRead.getNumberOfNodes()
423 ligne = "Number of nodes: %d" % nbnodes
426 coords = meshMEDFileRead.getCoords()
429 nb_comp = coords.getNumberOfComponents()
430 l_info = coords.getInfoOnComponents()
432 l_info_0=["X", "Y", "Z"]
433 for id_node in (0, 1, nbnodes-1):
434 ligne += "\nNode #%6d:" % id_node
435 for iaux in range(nb_comp):
439 saux = l_info_0[iaux]
440 ligne += " %s" % saux
441 ligne += "=%f" % coords[id_node, iaux]
444 # 4. Exploration of every group of faces
446 tb_aux = np.zeros(nbnodes, dtype=np.bool)
448 #bathy = np.zeros(nbnodes, dtype=np.float)
450 bathy = coords[:,2].toNumPyArray()
453 # 4.1. Mesh of the group
455 mesh_of_the_group = meshMEDFileRead.getGroup(0, gr_face_name_tr, False)
456 nbr_cells = mesh_of_the_group.getNumberOfCells()
458 ligne = "\t. Number of cells: %d" % nbr_cells
461 # 4.2. Nodes of the meshes of the group
462 # Every node is flagged in tb_aux
465 for id_elem in range(nbr_cells):
466 l_nodes = mesh_of_the_group.getNodeIdsOfCell(id_elem)
468 for id_node in l_nodes:
469 tb_aux[id_node] = True
470 np_aux = tb_aux.nonzero()
473 ligne = "\t. Number of nodes for this group: %d" % len(np_aux[0])
475 #print ("np_aux:", np_aux)
477 # 4.3. Interpolation over the nodes of the meshes of the group
481 for nid in np_aux[0]:
483 vx.append(coords[nodeId, 0])
484 vy.append(coords[nodeId, 1])
488 vz = bathy_obj.GetAltitudesForPoints(vx, vy, interp_method)
495 v05z = np.percentile(vz, 0o5)
496 v95z = np.percentile(vz, 95)
499 ligne = ".. Minimum: %f" % minz
500 ligne += ", maximum: %f" % maxz
501 ligne += ", mean: %f\n" % meanz
502 ligne += ".. stdeviation: %f" % stdz
503 ligne += ", v05z: %f" % v05z
504 ligne += ", v95z: %f" % v95z
507 # 4.4. Storage of the z and of the statistics for this region
509 statz_gr_face_name = (minz, maxz, meanz, stdz, v05z, v95z)
511 for iaux, nodeId in enumerate(np_aux[0]):
512 bathy[nodeId] = vz[iaux]
515 # During the interpolation, if no value is available over a node, a default value
516 # is set: -9999. It has no importance for the final computation, but if the field
517 # or the mesh is displayed, it makes huge gap. To prevent this artefact, a more
518 # convenient "undefined" value is set. This new undefined value is given by the user.
520 # zUndefThreshold: the default value is -9000. It is tied with the value -9999. given
521 # by the interpolation when no value is defined.
523 zUndefThreshold = -9000.
525 ligne = "zUndefThreshold: %f" % zUndefThreshold
528 #print "bathy :\n", bathy
529 np_aux_z = (bathy < zUndefThreshold).nonzero()
531 ligne = ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
534 for iaux in np_aux_z[0]:
537 # 6. Option : xyz file
541 ligne = ".. Ecriture du champ de bathymétrie sur le fichier :\n%s" % fichierFonds
544 with open(fichierFonds, "w") as fo :
545 for nodeId in range(nbnodes):
546 ligne = "%11.3f %11.3f %11.3f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
550 # 7.1. Transformation of the bathymetry as a double array
552 bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
553 bathy_dd.setInfoOnComponents(["Z [m]"])
555 # 7.2. If requested, modification of the z coordinate
558 coords3D = ml.DataArrayDouble.Meld([coords[:,0:2], bathy_dd])
559 coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
560 #print "coords3D =\n", coords3D
561 meshMEDFileRead.setCoords(coords3D)
563 # 7.3. Writes the mesh
570 ligne = ".. Ecriture du maillage" + saux
571 ligne += " sur le fichier :\n%s" % fichierFMaillage
574 meshMEDFileRead.write(fichierFMaillage, 2)
576 # 7.4. Writes the field
578 med_field_name = "BOTTOM"
580 ligne = ".. Ecriture du champ '%s'" % med_field_name
583 #print "bathy_dd =\n", bathy_dd
584 fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES)
585 fieldOnNodes.setName(med_field_name)
586 fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0))
587 fieldOnNodes.setArray(bathy_dd)
588 # Ces valeurs d'instants sont mises pour assurer la lecture par TELEMAC
590 # numero d'itération : 0
591 # pas de numero d'ordre (-1)
592 fieldOnNodes.setTime(0.0, 0, -1)
594 fMEDFile_ch_d = ml.MEDFileField1TS()
595 fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
596 fMEDFile_ch_d.write(fichierFMaillage, 0)
598 return statz_gr_face_name