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
  28 import numpy as np
  29 import MEDLoader as ml
  30 import medcoupling as mc
  31
  32 # -----------------------------------------------------------------------------
  33
  34 from med import medfile
  35 from med import medmesh
  36 from med import medfield
  37 from med import medenum
  38
  39 # -----------------------------------------------------------------------------
  40
  41 import HYDROPy
  42
  43 class MyInterpolator( HYDROPy.HYDROData_IInterpolator ):
  44   """
  45   Class MyInterpolator
  46   """
  47   def __init__ (self) :
  48     """
  49 Constructor
  50     """
  51   def GetAltitudeForPoint( self, theCoordX, theCoordY ):
  52     """
  53     Function
  54     """
  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
  57     #z2 = (z - 74.0)*10
  58     z2 = z
  59     return z2
  60
  61 # -----------------------------------------------------------------------------
  62
  63
  64 def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=-100., regions_interp_method=None, m3d=True, xyzFile=False, verbose=False):
  65   """
  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.
  68
  69   In:
  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.
  73                       Key: face group name
  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
  78                     if integer :
  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.
  85   Out:
  86     statz: statistique for z
  87                       Key: face group name
  88                       Value: (minz, maxz, meanz, stdz, v05z, v95z)
  89   Out:
  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
  93   """
  94   statz = dict()
  95   erreur = 0
  96   message = ""
  97
  98   while not erreur:
  99
 100     if verbose:
 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
 106       print (ligne)
 107
 108     doc = HYDROPy.HYDROData_Document.Document()
 109     cas = doc.FindObjectByName(nomCas)
 110     print ( "cas : ", cas)
 111     custom_inter = MyInterpolator()
 112
 113     basename = fichierMaillage[:-4]
 114     fichierFMaillage = basename + 'F.med'
 115
 116     print ("dicoGroupeRegion = ", dicoGroupeRegion)
 117     ligne = "fichierMaillage  = %s" % fichierMaillage
 118     ligne += "\nfichierFMaillage = %s" % fichierFMaillage
 119     if xyzFile:
 120       fichierFonds = basename + '.xyz'
 121       ligne += "\nfichierFonds     = %s" % fichierFonds
 122     print (ligne)
 123 #
 124 # 1. Reads the mesh
 125 #
 126     meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
 127     if verbose:
 128       print (meshMEDFileRead)
 129 #
 130 # 2. Checks the names of the groups of faces
 131 #
 132     t_group_n = meshMEDFileRead.getGroupsNames()
 133     dicoGroupeRegion_0 = dict()
 134     nb_pb = 0
 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"
 139         nb_pb += 1
 140       else :
 141         dicoGroupeRegion_0[saux] =  dicoGroupeRegion[gr_face_name]
 142     if verbose:
 143       ligne = "Number of problems: %d" % nb_pb
 144       print (ligne)
 145 #
 146     if nb_pb > 0:
 147       if nb_pb == 1:
 148         message += "This group does"
 149       else:
 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
 156       erreur = 2
 157       break
 158 #
 159 # 3. Gets the information about the nodes
 160 #
 161     nbnodes = meshMEDFileRead.getNumberOfNodes()
 162     if verbose:
 163       ligne = "Number of nodes: %d" % nbnodes
 164       print (ligne)
 165 #
 166     coords = meshMEDFileRead.getCoords()
 167     #print (coords)
 168     if verbose:
 169       nb_comp = coords.getNumberOfComponents()
 170       l_info = coords.getInfoOnComponents()
 171       ligne = ""
 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):
 176           if l_info[iaux]:
 177             saux = l_info[iaux]
 178           else:
 179             saux = l_info_0[iaux]
 180           ligne += " %s" % saux
 181           ligne += "=%f" % coords[id_node, iaux]
 182       print (ligne)
 183 #
 184 # 4. Exploration of every group of faces
 185 #
 186     tb_aux = np.zeros(nbnodes, dtype=np.bool)
 187 #
 188     bathy = np.zeros(nbnodes, dtype=np.float)
 189     bathy.fill(zUndef)
 190 #
 191     for gr_face_name in dicoGroupeRegion_0:
 192 #
 193 #     4.1. Region connected to the group
 194 #
 195       nomreg = dicoGroupeRegion_0[gr_face_name]
 196       ligne = "------- Region: '%s'" % nomreg
 197       ligne += ", connected to group '%s'" % gr_face_name
 198       print (ligne)
 199       region = doc.FindObjectByName(nomreg)
 200 #
 201 #     4.2. Mesh of the group
 202 #
 203       mesh_of_the_group = meshMEDFileRead.getGroup(0, gr_face_name, False)
 204       nbr_cells = mesh_of_the_group.getNumberOfCells()
 205       if verbose:
 206         ligne = "\t. Number of cells: %d" % nbr_cells
 207         print (ligne)
 208 #
 209 #     4.3. Nodes of the meshes of the group
 210 #          Every node is flagged in tb_aux
 211 #
 212       tb_aux.fill(False)
 213       for id_elem in range(nbr_cells):
 214         l_nodes = mesh_of_the_group.getNodeIdsOfCell(id_elem)
 215         #print l_nodes
 216         for id_node in l_nodes:
 217           tb_aux[id_node] = True
 218       np_aux = tb_aux.nonzero()
 219       if len(np_aux[0]):
 220         if verbose:
 221           ligne = "\t. Number of nodes for this group: %d" % len(np_aux[0])
 222           print (ligne)
 223       #print ("np_aux:", np_aux)
 224 #
 225 #     4.4. Interpolation over the nodes of the meshes of the group
 226 #
 227       vx = list()
 228       vy = list()
 229       for nid in np_aux[0]:
 230         nodeId = nid.item()
 231         vx.append(coords[nodeId, 0])
 232         vy.append(coords[nodeId, 1])
 233       #print ("vx:\n", vx)
 234       #print ("vy:\n", vy)
 235 #
 236       interpolMethod = 0
 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
 242
 243       #print ('interp', interpolMethod )
 244       vz = cas.GetAltitudesForPoints(vx, vy, region, interpolMethod)
 245 #
 246       #print ("vz:\n", vz)
 247       minz = np.amin(vz)
 248       maxz = np.amax(vz)
 249       meanz = np.mean(vz)
 250       stdz = np.std(vz)
 251       v05z = np.percentile(vz, 0o5)
 252       v95z = np.percentile(vz, 95)
 253 #
 254       if verbose:
 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
 261         print (ligne)
 262 #
 263 #     4.5. Storage of the z and of the statistics for this region
 264 #
 265       statz[gr_face_name] = (minz, maxz, meanz, stdz, v05z, v95z)
 266 #
 267       for iaux, nodeId in enumerate(np_aux[0]):
 268         bathy[nodeId] = vz[iaux]
 269 #
 270 # 5. Minimum:
 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.
 275 #
 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.
 278 #
 279     zUndefThreshold = zUndef +10.
 280     if verbose:
 281       ligne = "zUndefThreshold: %f" % zUndefThreshold
 282       print (ligne)
 283 #
 284     #print "bathy :\n", bathy
 285     np_aux_z = (bathy < zUndefThreshold).nonzero()
 286     if verbose:
 287       ligne = ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
 288       print (ligne)
 289 #
 290 # 6. Option : xyz file
 291 #
 292     if xyzFile:
 293 #
 294       if verbose:
 295         ligne = ".. Ecriture du champ de bathymétrie sur le fichier :\n%s" % fichierFonds
 296         print (ligne)
 297 #
 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])
 301           fo.write(ligne)
 302 #
 303 # 7. Final MED file
 304 # 7.1. Transformation of the bathymetry as a double array
 305 #
 306     bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
 307     bathy_dd.setInfoOnComponents(["Z [m]"])
 308 #
 309 # 7.2. If requested, modification of the z coordinate
 310 #
 311     if m3d:
 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)
 316 #
 317 # 7.3. Writes the mesh
 318 #
 319     if verbose:
 320       if m3d:
 321         saux = " 3D"
 322       else:
 323         saux = ""
 324       ligne = ".. Ecriture du maillage" + saux
 325       ligne += " sur le fichier :\n%s" % fichierFMaillage
 326       print (ligne)
 327 #
 328     meshMEDFileRead.write(fichierFMaillage, 2)
 329 #
 330 # 7.4. Writes the field
 331 #
 332     med_field_name = "BOTTOM"
 333     if verbose:
 334       ligne = ".. Ecriture du champ '%s'" % med_field_name
 335       print (ligne)
 336 #
 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
 343 #   instant = 0.0
 344 #   numero d'itération : 0
 345 #   pas de numero d'ordre (-1)
 346     fieldOnNodes.setTime(0.0, 0, -1)
 347 #
 348     fMEDFile_ch_d = ml.MEDFileField1TS()
 349     fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
 350     fMEDFile_ch_d.write(fichierFMaillage, 0)
 351 #
 352     break
 353 #
 354   if erreur:
 355     print(message)
 356 #
 357   return statz
 358
 359
 360 def interpolZ_B(bathyName, fichierMaillage, gr_face_name, zUndef=-100., interp_method=0, m3d=True, xyzFile=False, verbose=False):
 361   """
 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.
 364
 365   In:
 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.
 376   Out:
 377     if OK : statz_gr_face_name: statistic for z: (minz, maxz, meanz, stdz, v05z, v95z)
 378         else FALSE
 379   Out:
 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
 383   """
 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")
 389     return False
 390
 391   basename = fichierMaillage[:-4]
 392   fichierFMaillage = basename + 'F.med'
 393
 394   ligne = "fichierMaillage  = %s" % fichierMaillage
 395   ligne += "\nfichierFMaillage = %s" % fichierFMaillage
 396
 397   if xyzFile:
 398     fichierFonds = basename + '.xyz'
 399     ligne += "\nfichierFonds     = %s" % fichierFonds
 400   print (ligne)
 401 #
 402 # 1. Reads the mesh
 403 #
 404   meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
 405   if verbose:
 406     print (meshMEDFileRead)
 407 #
 408 # 2. Checks the names of the groups of faces
 409 #
 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")
 414     return False
 415 #
 416 # 3. Gets the information about the nodes
 417 #
 418   nbnodes = meshMEDFileRead.getNumberOfNodes()
 419   if verbose:
 420     ligne = "Number of nodes: %d" % nbnodes
 421     print (ligne)
 422 #
 423   coords = meshMEDFileRead.getCoords()
 424   #print(coords[:,2])
 425   if verbose:
 426     nb_comp = coords.getNumberOfComponents()
 427     l_info = coords.getInfoOnComponents()
 428     ligne = ""
 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):
 433         if l_info[iaux]:
 434           saux = l_info[iaux]
 435         else:
 436           saux = l_info_0[iaux]
 437         ligne += " %s" % saux
 438         ligne += "=%f" % coords[id_node, iaux]
 439     print (ligne)
 440 #
 441 # 4. Exploration of every group of faces
 442 #
 443   tb_aux = np.zeros(nbnodes, dtype=np.bool)
 444 #
 445   bathy = np.zeros(nbnodes, dtype=np.float)
 446   bathy.fill(zUndef)
 447   if (coords.getNumberOfComponents() >2):
 448     bathy = coords[:,2].toNumPyArray()
 449   else:
 450     print("=== WARNING! === Mesh has no altitude component z, z will be filled with zUndef = %s outside the group!"%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 nid in np_aux[0]:
 482     nodeId = nid.item()
 483     vx.append(coords[nodeId, 0])
 484     vy.append(coords[nodeId, 1])
 485   #print ("vx:\n", vx)
 486   #print ("vy:\n", vy)
 487 #
 488   vz = bathy_obj.GetAltitudesForPoints(vx, vy, interp_method)
 489 #
 490   #print ("vz:\n", vz)
 491   minz = np.amin(vz)
 492   maxz = np.amax(vz)
 493   meanz = np.mean(vz)
 494   stdz = np.std(vz)
 495   v05z = np.percentile(vz, 0o5)
 496   v95z = np.percentile(vz, 95)
 497 #
 498   if verbose:
 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
 505     print (ligne)
 506 #
 507 # 4.4. Storage of the z and of the statistics for this region
 508 #
 509   statz_gr_face_name = (minz, maxz, meanz, stdz, v05z, v95z)
 510 #
 511   for iaux, nodeId in enumerate(np_aux[0]):
 512     bathy[nodeId] = vz[iaux]
 513 #
 514 # 5. Minimum:
 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.
 519 #
 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.
 522 #
 523   zUndefThreshold = zUndef + 10.
 524   if verbose:
 525     ligne = "zUndefThreshold: %f" % zUndefThreshold
 526     print (ligne)
 527 #
 528   #print "bathy :\n", bathy
 529   np_aux_z = (bathy < zUndefThreshold).nonzero()
 530   if verbose:
 531     ligne = ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
 532     print (ligne)
 533 #
 534 # 6. Option : xyz file
 535 #
 536   if xyzFile:
 537     if verbose:
 538       ligne = ".. Ecriture du champ de bathymétrie sur le fichier :\n%s" % fichierFonds
 539       print (ligne)
 540 #
 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])
 544         fo.write(ligne)
 545 #
 546 # 7. Final MED file
 547 # 7.1. Transformation of the bathymetry as a double array
 548 #
 549   bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
 550   bathy_dd.setInfoOnComponents(["Z [m]"])
 551 #
 552 # 7.2. If requested, modification of the z coordinate
 553 #
 554   if m3d:
 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)
 559 #
 560 # 7.3. Writes the mesh
 561 #
 562   if verbose:
 563     if m3d:
 564       saux = " 3D"
 565     else:
 566       saux = ""
 567     ligne = ".. Ecriture du maillage" + saux
 568     ligne += " sur le fichier :\n%s" % fichierFMaillage
 569     print (ligne)
 570 #
 571   meshMEDFileRead.write(fichierFMaillage, 2)
 572 #
 573 # 7.4. Writes the field
 574 #
 575   med_field_name = "BOTTOM"
 576   if verbose:
 577     ligne = ".. Ecriture du champ '%s'" % med_field_name
 578     print (ligne)
 579 #
 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
 586 # instant = 0.0
 587 # numero d'itération : 0
 588 # pas de numero d'ordre (-1)
 589   fieldOnNodes.setTime(0.0, 0, -1)
 590 #
 591   fMEDFile_ch_d = ml.MEDFileField1TS()
 592   fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
 593   fMEDFile_ch_d.write(fichierFMaillage, 0)
 594 #
 595   return statz_gr_face_name