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=-100., 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).
76 default value is -100.
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 zUndef +10. It is tied with the value -100. given
277 # by the interpolation when no value is defined.
279 zUndefThreshold = zUndef +10.
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 # 6. Option : xyz file
295 ligne = ".. Ecriture du champ de bathymétrie sur le fichier :\n%s" % fichierFonds
298 with open(fichierFonds, "w") as fo :
299 for nodeId in range(nbnodes):
300 ligne = "%11.3f %11.3f %11.3f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
304 # 7.1. Transformation of the bathymetry as a double array
306 bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
307 bathy_dd.setInfoOnComponents(["Z [m]"])
309 # 7.2. If requested, modification of the z coordinate
312 coords3D = ml.DataArrayDouble.Meld([coords[:,0:2], bathy_dd])
313 coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
314 #print "coords3D =\n", coords3D
315 meshMEDFileRead.setCoords(coords3D)
317 # 7.3. Writes the mesh
324 ligne = ".. Ecriture du maillage" + saux
325 ligne += " sur le fichier :\n%s" % fichierFMaillage
328 meshMEDFileRead.write(fichierFMaillage, 2)
330 # 7.4. Writes the field
332 med_field_name = "BOTTOM"
334 ligne = ".. Ecriture du champ '%s'" % med_field_name
337 #print "bathy_dd =\n", bathy_dd
338 fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES)
339 fieldOnNodes.setName(med_field_name)
340 fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0))
341 fieldOnNodes.setArray(bathy_dd)
342 # Ces valeurs d'instants sont mises pour assurer la lecture par TELEMAC
344 # numero d'itération : 0
345 # pas de numero d'ordre (-1)
346 fieldOnNodes.setTime(0.0, 0, -1)
348 fMEDFile_ch_d = ml.MEDFileField1TS()
349 fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
350 fMEDFile_ch_d.write(fichierFMaillage, 0)
360 def interpolZ_B(bathyName, fichierMaillage, gr_face_name, zUndef=-100., interp_method=0, m3d=True, xyzFile=False, verbose=False):
362 interpolZ_B takes a 2D (x,y) mesh and calls the active instance of module HYDRO
363 to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node.
366 bathyName: Bathymetry Name in module HYDRO
367 fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case
368 gr_face_name: face group name
369 zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
370 default value is -100.
371 interp_method: interpolation method
372 0 = nearest point on bathymetry
373 1 = linear interpolation
374 m3d: True/False to produce a 3D mesh. Default is True.
375 xyzFile: True/False to write an ascii file with xyz for every node. Default is False.
377 if OK : statz_gr_face_name: statistic for z: (minz, maxz, meanz, stdz, v05z, v95z)
380 return <fichierMaillage>F.med : med file with Z value in a field "BOTTOM"
381 Option: Z value also in Z coordinate if m3D is true
382 return <fichierMaillage>.xyz : text file with X, Y, Z values
384 doc = HYDROPy.HYDROData_Document.Document()
385 bathy_obj = doc.FindObjectByName(bathyName)
386 print( "bathy : ", bathy_obj)
387 if bathy_obj is None:
388 print ( "bathy is None")
391 basename = fichierMaillage[:-4]
392 fichierFMaillage = basename + 'F.med'
394 ligne = "fichierMaillage = %s" % fichierMaillage
395 ligne += "\nfichierFMaillage = %s" % fichierFMaillage
398 fichierFonds = basename + '.xyz'
399 ligne += "\nfichierFonds = %s" % fichierFonds
404 meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
406 print (meshMEDFileRead)
408 # 2. Checks the names of the groups of faces
410 t_group_n = meshMEDFileRead.getGroupsNames()
411 gr_face_name_tr = gr_face_name.strip()
412 if gr_face_name_tr not in t_group_n:
413 print("Group not found")
416 # 3. Gets the information about the nodes
418 nbnodes = meshMEDFileRead.getNumberOfNodes()
420 ligne = "Number of nodes: %d" % nbnodes
423 coords = meshMEDFileRead.getCoords()
426 nb_comp = coords.getNumberOfComponents()
427 l_info = coords.getInfoOnComponents()
429 l_info_0=["X", "Y", "Z"]
430 for id_node in (0, 1, nbnodes-1):
431 ligne += "\nNode #%6d:" % id_node
432 for iaux in range(nb_comp):
436 saux = l_info_0[iaux]
437 ligne += " %s" % saux
438 ligne += "=%f" % coords[id_node, iaux]
441 # 4. Exploration of every group of faces
443 tb_aux = np.zeros(nbnodes, dtype=np.bool)
445 bathy = np.zeros(nbnodes, dtype=np.float)
447 if (coords.getNumberOfComponents() >2):
448 bathy = coords[:,2].toNumPyArray()
450 print("=== WARNING! === Mesh has no altitude component z, z will be filled with zUndef = %s outside the group!"%zUndef)
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 zUndef +10. It is tied with the value -100. given
521 # by the interpolation when no value is defined.
523 zUndefThreshold = zUndef + 10.
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 # 6. Option : xyz file
538 ligne = ".. Ecriture du champ de bathymétrie sur le fichier :\n%s" % fichierFonds
541 with open(fichierFonds, "w") as fo :
542 for nodeId in range(nbnodes):
543 ligne = "%11.3f %11.3f %11.3f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
547 # 7.1. Transformation of the bathymetry as a double array
549 bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
550 bathy_dd.setInfoOnComponents(["Z [m]"])
552 # 7.2. If requested, modification of the z coordinate
555 coords3D = ml.DataArrayDouble.Meld([coords[:,0:2], bathy_dd])
556 coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
557 #print "coords3D =\n", coords3D
558 meshMEDFileRead.setCoords(coords3D)
560 # 7.3. Writes the mesh
567 ligne = ".. Ecriture du maillage" + saux
568 ligne += " sur le fichier :\n%s" % fichierFMaillage
571 meshMEDFileRead.write(fichierFMaillage, 2)
573 # 7.4. Writes the field
575 med_field_name = "BOTTOM"
577 ligne = ".. Ecriture du champ '%s'" % med_field_name
580 #print "bathy_dd =\n", bathy_dd
581 fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES)
582 fieldOnNodes.setName(med_field_name)
583 fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0))
584 fieldOnNodes.setArray(bathy_dd)
585 # Ces valeurs d'instants sont mises pour assurer la lecture par TELEMAC
587 # numero d'itération : 0
588 # pas de numero d'ordre (-1)
589 fieldOnNodes.setTime(0.0, 0, -1)
591 fMEDFile_ch_d = ml.MEDFileField1TS()
592 fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
593 fMEDFile_ch_d.write(fichierFMaillage, 0)
595 return statz_gr_face_name