]> SALOME platform Git repositories - modules/hydro.git/blob - src/HYDROTools/interpolZ.py
Salome HOME
Merge branch 'BR_H2018_3' into BR_2018_V8_5
[modules/hydro.git] / src / HYDROTools / interpolZ.py
1 # -*- coding: utf-8 -*-
2 """
3 Example of use case of interpolZ with the default values:
4
5 # --- case name in HYDRO
6 nomCas = 'inondation1'
7
8 # --- med file 2D(x,y) of the case produced by SMESH
9 fichierMaillage = '/home/B27118/projets/LNHE/garonne/inondation1.med'
10
11 # --- dictionary [med group name] = region name
12 dicoGroupeRegion= dict(litMineur          = 'inondation1_litMineur',
13                        litMajeurDroite    = 'inondation1_riveDroite',
14                        litMajeurGauche    = 'inondation1_riveGauche',
15
16 # --- Z interpolation on the bathymety/altimetry on the mesh nodes
17 interpolZ(nomCas, fichierMaillage, dicoGroupeRegion)
18 """
19 __revision__ = "V3.01"
20 # -----------------------------------------------------------------------------
21
22 # -----------------------------------------------------------------------------
23
24 import salome
25
26 salome.salome_init()
27 theStudy = salome.myStudy
28 theStudyId = salome.myStudyId
29
30 import numpy as np
31 import MEDLoader as ml
32 import MEDCoupling as mc
33
34 # -----------------------------------------------------------------------------
35
36 from med import medfile
37 from med import medmesh
38 from med import medfield
39 from med import medenum
40
41 # -----------------------------------------------------------------------------
42
43 import HYDROPy
44
45 class MyInterpolator( HYDROPy.HYDROData_IInterpolator ):
46   """
47   Class MyInterpolator
48   """
49   def __init__ (self) :
50     """
51 Constructor
52     """
53   def GetAltitudeForPoint( self, theCoordX, theCoordY ):
54     """
55     Function
56     """
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
59     #z2 = (z - 74.0)*10
60     z2 = z
61     return z2
62
63 # -----------------------------------------------------------------------------
64
65
66 def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., regions_interp_method=None, m3d=True, xyzFile=False, verbose=False):
67   """
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.
70
71   In:
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.
75                       Key: face group name
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).
78                      default value is 90.
79     interpolMethod: integer value or dict
80                     if integer : 
81                     0 = nearest point on bathymetry
82                     1 = linear interpolation
83                     if dict : key is a region name, value is a type of interpolation (0 or 1)
84                     if None : default type used (0)
85     m3d: True/False to produce a 3D mesh. Default is False.
86     xyzFile: True/False to write an ascii file with xyz for every node. Default is False.
87   Out:
88     statz: statistique for z
89                       Key: face group name
90                       Value: (minz, maxz, meanz, stdz, v05z, v95z)
91   Out:
92     return <fichierMaillage>F.med : med file with Z value in a field "BOTTOM"
93                                     Option: Z value also in Z coordinate if m3D is true
94     return <fichierMaillage>.xyz : text file with X, Y, Z values
95   """
96   statz = dict()
97   erreur = 0
98   message = ""
99
100   while not erreur:
101
102     if verbose:
103       ligne = "nomCas: %s" % nomCas
104       ligne += "\ninterpolMethods: " % regions_interp_method
105       if (zUndef != None ):
106         ligne += "\nzUndef: %f" % zUndef
107       ligne += "\nm3d: %d" % m3d
108       print (ligne)
109
110     doc = HYDROPy.HYDROData_Document.Document(theStudyId)
111     cas = doc.FindObjectByName(nomCas)
112     print ( "cas : ", cas)
113     custom_inter = MyInterpolator()
114
115     basename = fichierMaillage[:-4]
116     fichierFMaillage = basename + 'F.med'
117
118     print ("dicoGroupeRegion = ", dicoGroupeRegion)
119     ligne = "fichierMaillage  = %s" % fichierMaillage
120     ligne += "\nfichierFMaillage = %s" % fichierFMaillage
121     if xyzFile:
122       fichierFonds = basename + '.xyz'
123       ligne += "\nfichierFonds     = %s" % fichierFonds
124     print (ligne)
125 #
126 # 1. Reads the mesh
127 #
128     meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
129     if verbose:
130       print (meshMEDFileRead)
131 #
132 # 2. Checks the names of the groups of faces
133 #
134     t_group_n = meshMEDFileRead.getGroupsNames()
135     dicoGroupeRegion_0 = dict()
136     nb_pb = 0
137     for gr_face_name in dicoGroupeRegion:
138       saux = gr_face_name.strip()
139       if saux not in t_group_n:
140         message += "Group: '" + gr_face_name + "'\n"
141         nb_pb += 1
142       else :
143         dicoGroupeRegion_0[saux] =  dicoGroupeRegion[gr_face_name]
144     if verbose:
145       ligne = "Number of problems: %d" % nb_pb
146       print (ligne)
147 #
148     if nb_pb > 0:
149       if nb_pb == 1:
150         message += "This group does"
151       else:
152         message += "These %d groups do" % nb_pb
153       message += " not belong to the mesh.\n"
154       message += "Please check the names of the group(s) of faces corresponding to each region of the HYDRO case.\n"
155       message += "Groups for this file:\n"
156       for group_n in t_group_n :
157         message += "'%s'\n" % group_n
158       erreur = 2
159       break
160 #
161 # 3. Gets the information about the nodes
162 #
163     nbnodes = meshMEDFileRead.getNumberOfNodes()
164     if verbose:
165       ligne = "Number of nodes: %d" % nbnodes
166       print (ligne)
167 #
168     coords = meshMEDFileRead.getCoords()
169     #print coords
170     if verbose:
171       nb_comp = coords.getNumberOfComponents()
172       l_info = coords.getInfoOnComponents()
173       ligne = ""
174       l_info_0=["X", "Y", "Z"]
175       for id_node in (0, 1, nbnodes-1):
176         ligne += "\nNode #%6d:" % id_node
177         for iaux in range(nb_comp):
178           if l_info[iaux]:
179             saux = l_info[iaux]
180           else:
181             saux = l_info_0[iaux]
182           ligne += " %s" % saux
183           ligne += "=%f" % coords[id_node, iaux]
184       print (ligne)
185 #
186 # 4. Exploration of every group of faces
187 #
188     tb_aux = np.zeros(nbnodes, dtype=np.bool)
189 #
190     bathy = np.zeros(nbnodes, dtype=np.float)
191     bathy.fill(zUndef)
192 #
193     for gr_face_name in dicoGroupeRegion_0:
194 #
195 #     4.1. Region connected to the group
196 #
197       nomreg = dicoGroupeRegion_0[gr_face_name]
198       ligne = "------- Region: '%s'" % nomreg
199       ligne += ", connected to group '%s'" % gr_face_name
200       print (ligne)
201       region = doc.FindObjectByName(nomreg)
202 #
203 #     4.2. Mesh of the group
204 #
205       mesh_of_the_group = meshMEDFileRead.getGroup(0, gr_face_name, False)
206       nbr_cells = mesh_of_the_group.getNumberOfCells()
207       if verbose:
208         ligne = "\t. Number of cells: %d" % nbr_cells
209         print (ligne)
210 #
211 #     4.3. Nodes of the meshes of the group
212 #          Every node is flagged in tb_aux
213 #
214       tb_aux.fill(False)
215       for id_elem in range(nbr_cells):
216         l_nodes = mesh_of_the_group.getNodeIdsOfCell(id_elem)
217         #print l_nodes
218         for id_node in l_nodes:
219           tb_aux[id_node] = True
220       np_aux = tb_aux.nonzero()
221       if len(np_aux[0]):
222         if verbose:
223           ligne = "\t. Number of nodes for this group: %d" % len(np_aux[0])
224           print (ligne)
225       #print ("np_aux:", np_aux)
226 #
227 #     4.4. Interpolation over the nodes of the meshes of the group
228 #
229       vx = list()
230       vy = list()
231       for nodeId in np_aux[0]:
232         vx.append(coords[nodeId, 0])
233         vy.append(coords[nodeId, 1])
234       #print ("vx:\n", vx)
235       #print ("vy:\n", vy)
236 #
237       interpolMethod = 0
238       if regions_interp_method is not None:
239         if isinstance(regions_interp_method, dict) and nomreg in regions_interp_method.keys():
240           interpolMethod = int(regions_interp_method[nomreg])
241         elif isinstance(regions_interp_method, int):
242           interpolMethod = regions_interp_method
243
244       #print ('interp', interpolMethod )
245       vz = cas.GetAltitudesForPoints(vx, vy, region, interpolMethod)
246 #
247       #print ("vz:\n", vz)
248       minz = np.amin(vz)
249       maxz = np.amax(vz)
250       meanz = np.mean(vz)
251       stdz = np.std(vz)
252       v05z = np.percentile(vz, 05)
253       v95z = np.percentile(vz, 95)
254 #
255       if verbose:
256         ligne = ".. Minimum: %f" % minz
257         ligne += ", maximum: %f" % maxz
258         ligne += ", mean: %f\n" % meanz
259         ligne += ".. stdeviation: %f" % stdz
260         ligne += ", v05z: %f" % v05z
261         ligne += ", v95z: %f" % v95z
262         print (ligne)
263 #
264 #     4.5. Storage of the z and of the statistics for this region
265 #
266       statz[gr_face_name] = (minz, maxz, meanz, stdz, v05z, v95z)
267 #
268       for iaux, nodeId in enumerate(np_aux[0]):
269         bathy[nodeId] = vz[iaux]
270 #
271 # 5. Minimum:
272 #    During the interpolation, if no value is available over a node, a default value
273 #    is set: -9999. It has no importance for the final computation, but if the field
274 #    or the mesh is displayed, it makes huge gap. To prevent this artefact, a more
275 #    convenient "undefined" value is set. This new undefined value is given by the user.
276 #
277 #    zUndefThreshold: the default value is -9000. It is tied with the value -9999. given
278 #                     by the interpolation when no value is defined.
279 #
280     zUndefThreshold = -9000.
281     if verbose:
282       ligne = "zUndefThreshold: %f" % zUndefThreshold
283       print (ligne)
284 #
285     #print "bathy :\n", bathy
286     np_aux_z = (bathy < zUndefThreshold).nonzero()
287     if verbose:
288       ligne = ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
289       print (ligne)
290     if len(np_aux_z[0]):
291       for iaux in np_aux_z[0]:
292         bathy[iaux] = zUndef
293 #
294 # 6. Option : xyz file
295 #
296     if xyzFile:
297 #
298       if verbose:
299         ligne = ".. Ecriture du champ de bathymétrie sur le fichier :\n%s" % fichierFonds
300         print (ligne)
301 #
302       with open(fichierFonds, "w") as fo :
303         for nodeId in range(nbnodes):
304           ligne = "%11.3f %11.3f %11.3f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
305           fo.write(ligne)
306 #
307 # 7. Final MED file
308 # 7.1. Transformation of the bathymetry as a double array
309 #
310     bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
311     bathy_dd.setInfoOnComponents(["Z [m]"])
312 #
313 # 7.2. If requested, modification of the z coordinate
314 #
315     if m3d:
316       coords3D = ml.DataArrayDouble.Meld([coords[:,0:2], bathy_dd])
317       coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
318       #print "coords3D =\n", coords3D
319       meshMEDFileRead.setCoords(coords3D)
320 #
321 # 7.3. Writes the mesh
322 #
323     if verbose:
324       if m3d:
325         saux = " 3D"
326       else:
327         saux = ""
328       ligne = ".. Ecriture du maillage" + saux
329       ligne += " sur le fichier :\n%s" % fichierFMaillage
330       print (ligne)
331 #
332     meshMEDFileRead.write(fichierFMaillage, 2)
333 #
334 # 7.4. Writes the field
335 #
336     med_field_name = "BOTTOM"
337     if verbose:
338       ligne = ".. Ecriture du champ '%s'" % med_field_name
339       print (ligne)
340 #
341     #print "bathy_dd =\n", bathy_dd
342     fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES)
343     fieldOnNodes.setName(med_field_name)
344     fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0))
345     fieldOnNodes.setArray(bathy_dd)
346 #   Ces valeurs d'instants sont mises pour assurer la lecture par TELEMAC
347 #   instant = 0.0
348 #   numero d'itération : 0
349 #   pas de numero d'ordre (-1)
350     fieldOnNodes.setTime(0.0, 0, -1)
351 #
352     fMEDFile_ch_d = ml.MEDFileField1TS()
353     fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
354     fMEDFile_ch_d.write(fichierFMaillage, 0)
355 #
356     break
357 #
358   if erreur:
359     print message
360 #
361   return statz
362
363
364 def interpolZ_B(bathyName, fichierMaillage, gr_face_name, zUndef=90., interp_method=0, m3d=True, xyzFile=False, verbose=False):
365   """
366   interpolZ_B takes a 2D (x,y) mesh and calls the active instance of module HYDRO
367   to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node.
368
369   In:
370     bathyName: Bathymetry Name in module HYDRO
371     fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case
372     gr_face_name:  face group name
373     zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
374                      default value is 90.
375     interp_method: interpolation method
376                     0 = nearest point on bathymetry
377                     1 = linear interpolation
378     m3d: True/False to produce a 3D mesh. Default is True.
379     xyzFile: True/False to write an ascii file with xyz for every node. Default is False.
380   Out:
381     if OK : statz_gr_face_name: statistic for z: (minz, maxz, meanz, stdz, v05z, v95z)
382         else FALSE
383   Out:
384     return <fichierMaillage>F.med : med file with Z value in a field "BOTTOM"
385                                     Option: Z value also in Z coordinate if m3D is true
386     return <fichierMaillage>.xyz : text file with X, Y, Z values
387   """
388   doc = HYDROPy.HYDROData_Document.Document(theStudyId)
389   bathy_obj = doc.FindObjectByName(bathyName)
390   print ( "bathy : ", bathy_obj)
391   if bathy_obj is None:
392     print ( "bathy is None")
393     return False
394
395   basename = fichierMaillage[:-4]
396   fichierFMaillage = basename + 'F.med'
397
398   ligne = "fichierMaillage  = %s" % fichierMaillage
399   ligne += "\nfichierFMaillage = %s" % fichierFMaillage
400   
401   if xyzFile:
402     fichierFonds = basename + '.xyz'
403     ligne += "\nfichierFonds     = %s" % fichierFonds
404   print (ligne)
405 #
406 # 1. Reads the mesh
407 #
408   meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
409   if verbose:
410     print (meshMEDFileRead)
411 #
412 # 2. Checks the names of the groups of faces
413 #
414   t_group_n = meshMEDFileRead.getGroupsNames()
415   gr_face_name_tr = gr_face_name.strip()
416   if gr_face_name_tr not in t_group_n:
417     print "Group not found"
418     return False          
419 #
420 # 3. Gets the information about the nodes
421 #
422   nbnodes = meshMEDFileRead.getNumberOfNodes()
423   if verbose:
424     ligne = "Number of nodes: %d" % nbnodes
425     print (ligne)
426 #
427   coords = meshMEDFileRead.getCoords()
428   #print coords
429   if verbose:
430     nb_comp = coords.getNumberOfComponents()
431     l_info = coords.getInfoOnComponents()
432     ligne = ""
433     l_info_0=["X", "Y", "Z"]
434     for id_node in (0, 1, nbnodes-1):
435       ligne += "\nNode #%6d:" % id_node
436       for iaux in range(nb_comp):
437         if l_info[iaux]:
438           saux = l_info[iaux]
439         else:
440           saux = l_info_0[iaux]
441         ligne += " %s" % saux
442         ligne += "=%f" % coords[id_node, iaux]
443     print (ligne)
444 #
445 # 4. Exploration of every group of faces
446 #
447   tb_aux = np.zeros(nbnodes, dtype=np.bool)
448 #
449   bathy = np.zeros(nbnodes, dtype=np.float)
450   bathy.fill(zUndef)
451
452 #
453 # 4.1. Mesh of the group
454 #
455   mesh_of_the_group = meshMEDFileRead.getGroup(0, gr_face_name_tr, False)
456   nbr_cells = mesh_of_the_group.getNumberOfCells()
457   if verbose:
458     ligne = "\t. Number of cells: %d" % nbr_cells
459     print (ligne)
460 #
461 # 4.2. Nodes of the meshes of the group
462 #      Every node is flagged in tb_aux
463 #
464   tb_aux.fill(False)
465   for id_elem in range(nbr_cells):
466     l_nodes = mesh_of_the_group.getNodeIdsOfCell(id_elem)
467     #print l_nodes
468     for id_node in l_nodes:
469       tb_aux[id_node] = True
470   np_aux = tb_aux.nonzero()
471   if len(np_aux[0]):
472     if verbose:
473       ligne = "\t. Number of nodes for this group: %d" % len(np_aux[0])
474       print (ligne)
475   #print ("np_aux:", np_aux)
476 #
477 # 4.3. Interpolation over the nodes of the meshes of the group
478 #
479   vx = list()
480   vy = list()
481   for nodeId in np_aux[0]:
482     vx.append(coords[nodeId, 0])
483     vy.append(coords[nodeId, 1])
484   #print ("vx:\n", vx)
485   #print ("vy:\n", vy)
486 #
487   vz = bathy_obj.GetAltitudesForPoints(vx, vy, interp_method)
488 #
489   #print ("vz:\n", vz)
490   minz = np.amin(vz)
491   maxz = np.amax(vz)
492   meanz = np.mean(vz)
493   stdz = np.std(vz)
494   v05z = np.percentile(vz, 05)
495   v95z = np.percentile(vz, 95)
496 #
497   if verbose:
498     ligne = ".. Minimum: %f" % minz
499     ligne += ", maximum: %f" % maxz
500     ligne += ", mean: %f\n" % meanz
501     ligne += ".. stdeviation: %f" % stdz
502     ligne += ", v05z: %f" % v05z
503     ligne += ", v95z: %f" % v95z
504     print (ligne)
505 #
506 # 4.4. Storage of the z and of the statistics for this region
507 #
508   statz_gr_face_name = (minz, maxz, meanz, stdz, v05z, v95z)
509 #
510   for iaux, nodeId in enumerate(np_aux[0]):
511     bathy[nodeId] = vz[iaux]
512 #
513 # 5. Minimum:
514 #    During the interpolation, if no value is available over a node, a default value
515 #    is set: -9999. It has no importance for the final computation, but if the field
516 #    or the mesh is displayed, it makes huge gap. To prevent this artefact, a more
517 #    convenient "undefined" value is set. This new undefined value is given by the user.
518 #
519 #    zUndefThreshold: the default value is -9000. It is tied with the value -9999. given
520 #                     by the interpolation when no value is defined.
521 #
522   zUndefThreshold = -9000.
523   if verbose:
524     ligne = "zUndefThreshold: %f" % zUndefThreshold
525     print (ligne)
526 #
527   #print "bathy :\n", bathy
528   np_aux_z = (bathy < zUndefThreshold).nonzero()
529   if verbose:
530     ligne = ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
531     print (ligne)
532   if len(np_aux_z[0]):
533     for iaux in np_aux_z[0]:
534       bathy[iaux] = zUndef
535 #
536 # 6. Option : xyz file
537 #
538   if xyzFile:
539     if verbose:
540       ligne = ".. Ecriture du champ de bathymétrie sur le fichier :\n%s" % fichierFonds
541       print (ligne)
542 #
543     with open(fichierFonds, "w") as fo :
544       for nodeId in range(nbnodes):
545         ligne = "%11.3f %11.3f %11.3f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
546         fo.write(ligne)
547 #
548 # 7. Final MED file
549 # 7.1. Transformation of the bathymetry as a double array
550 #
551   bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
552   bathy_dd.setInfoOnComponents(["Z [m]"])
553 #
554 # 7.2. If requested, modification of the z coordinate
555 #
556   if m3d:
557     coords3D = ml.DataArrayDouble.Meld([coords[:,0:2], bathy_dd])
558     coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
559     #print "coords3D =\n", coords3D
560     meshMEDFileRead.setCoords(coords3D)
561 #
562 # 7.3. Writes the mesh
563 #
564   if verbose:
565     if m3d:
566       saux = " 3D"
567     else:
568       saux = ""
569     ligne = ".. Ecriture du maillage" + saux
570     ligne += " sur le fichier :\n%s" % fichierFMaillage
571     print (ligne)
572 #
573   meshMEDFileRead.write(fichierFMaillage, 2)
574 #
575 # 7.4. Writes the field
576 #
577   med_field_name = "BOTTOM"
578   if verbose:
579     ligne = ".. Ecriture du champ '%s'" % med_field_name
580     print (ligne)
581 #
582   #print "bathy_dd =\n", bathy_dd
583   fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES)
584   fieldOnNodes.setName(med_field_name)
585   fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0))
586   fieldOnNodes.setArray(bathy_dd)
587 # Ces valeurs d'instants sont mises pour assurer la lecture par TELEMAC
588 # instant = 0.0
589 # numero d'itération : 0
590 # pas de numero d'ordre (-1)
591   fieldOnNodes.setTime(0.0, 0, -1)
592 #
593   fMEDFile_ch_d = ml.MEDFileField1TS()
594   fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
595   fMEDFile_ch_d.write(fichierFMaillage, 0)
596 #
597   return statz_gr_face_name