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   statz = dict()
 385   doc = HYDROPy.HYDROData_Document.Document()
 386   bathy_obj = doc.FindObjectByName(bathyName)
 387   print( "bathy : ", bathy_obj)
 388   if bathy_obj is None:
 389     print ( "bathy is None")
 390     return False
 391
 392   basename = fichierMaillage[:-4]
 393   fichierFMaillage = basename + 'F.med'
 394
 395   ligne = "fichierMaillage  = %s" % fichierMaillage
 396   ligne += "\nfichierFMaillage = %s" % fichierFMaillage
 397
 398   if xyzFile:
 399     fichierFonds = basename + '.xyz'
 400     ligne += "\nfichierFonds     = %s" % fichierFonds
 401   print (ligne)
 402 #
 403 # 1. Reads the mesh
 404 #
 405   meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
 406   if verbose:
 407     print (meshMEDFileRead)
 408 #
 409 # 2. Checks the names of the groups of faces
 410 #
 411   t_group_n = meshMEDFileRead.getGroupsNames()
 412   gr_face_name_tr = gr_face_name.strip()
 413   if gr_face_name_tr not in t_group_n:
 414     print("Group not found")
 415     return False
 416 #
 417 # 3. Gets the information about the nodes
 418 #
 419   nbnodes = meshMEDFileRead.getNumberOfNodes()
 420   if verbose:
 421     ligne = "Number of nodes: %d" % nbnodes
 422     print (ligne)
 423 #
 424   coords = meshMEDFileRead.getCoords()
 425   #print(coords[:,2])
 426   if verbose:
 427     nb_comp = coords.getNumberOfComponents()
 428     l_info = coords.getInfoOnComponents()
 429     ligne = ""
 430     l_info_0=["X", "Y", "Z"]
 431     for id_node in (0, 1, nbnodes-1):
 432       ligne += "\nNode #%6d:" % id_node
 433       for iaux in range(nb_comp):
 434         if l_info[iaux]:
 435           saux = l_info[iaux]
 436         else:
 437           saux = l_info_0[iaux]
 438         ligne += " %s" % saux
 439         ligne += "=%f" % coords[id_node, iaux]
 440     print (ligne)
 441 #
 442 # 4. Exploration of every group of faces
 443 #
 444   tb_aux = np.zeros(nbnodes, dtype=np.bool)
 445 #
 446   bathy = np.zeros(nbnodes, dtype=np.float)
 447   bathy.fill(zUndef)
 448   if (coords.getNumberOfComponents() >2):
 449     bathy = coords[:,2].toNumPyArray()
 450   else:
 451     print("=== WARNING! === Mesh has no altitude component z, z will be filled with zUndef = %s outside the group!"%zUndef)
 452
 453 #
 454 # 4.1. Mesh of the group
 455 #
 456   mesh_of_the_group = meshMEDFileRead.getGroup(0, gr_face_name_tr, False)
 457   nbr_cells = mesh_of_the_group.getNumberOfCells()
 458   if verbose:
 459     ligne = "\t. Number of cells: %d" % nbr_cells
 460     print (ligne)
 461 #
 462 # 4.2. Nodes of the meshes of the group
 463 #      Every node is flagged in tb_aux
 464 #
 465   tb_aux.fill(False)
 466   for id_elem in range(nbr_cells):
 467     l_nodes = mesh_of_the_group.getNodeIdsOfCell(id_elem)
 468     #print l_nodes
 469     for id_node in l_nodes:
 470       tb_aux[id_node] = True
 471   np_aux = tb_aux.nonzero()
 472   if len(np_aux[0]):
 473     if verbose:
 474       ligne = "\t. Number of nodes for this group: %d" % len(np_aux[0])
 475       print (ligne)
 476   #print ("np_aux:", np_aux)
 477 #
 478 # 4.3. Interpolation over the nodes of the meshes of the group
 479 #
 480   vx = list()
 481   vy = list()
 482   for nid in np_aux[0]:
 483     nodeId = nid.item()
 484     vx.append(coords[nodeId, 0])
 485     vy.append(coords[nodeId, 1])
 486   #print ("vx:\n", vx)
 487   #print ("vy:\n", vy)
 488 #
 489   vz = bathy_obj.GetAltitudesForPoints(vx, vy, interp_method)
 490 #
 491   #print ("vz:\n", vz)
 492   minz = np.amin(vz)
 493   maxz = np.amax(vz)
 494   meanz = np.mean(vz)
 495   stdz = np.std(vz)
 496   v05z = np.percentile(vz, 0o5)
 497   v95z = np.percentile(vz, 95)
 498 #
 499   if verbose:
 500     ligne = ".. Minimum: %f" % minz
 501     ligne += ", maximum: %f" % maxz
 502     ligne += ", mean: %f\n" % meanz
 503     ligne += ".. stdeviation: %f" % stdz
 504     ligne += ", v05z: %f" % v05z
 505     ligne += ", v95z: %f" % v95z
 506     print (ligne)
 507 #
 508 # 4.4. Storage of the z and of the statistics for this region
 509 #
 510   statz[gr_face_name] = (minz, maxz, meanz, stdz, v05z, v95z)
 511 #
 512   for iaux, nodeId in enumerate(np_aux[0]):
 513     bathy[nodeId] = vz[iaux]
 514 #
 515 # 5. Minimum:
 516 #    During the interpolation, if no value is available over a node, a default value
 517 #    is set: -9999. It has no importance for the final computation, but if the field
 518 #    or the mesh is displayed, it makes huge gap. To prevent this artefact, a more
 519 #    convenient "undefined" value is set. This new undefined value is given by the user.
 520 #
 521 #    zUndefThreshold: the default value is zUndef +10. It is tied with the value -100. given
 522 #                     by the interpolation when no value is defined.
 523 #
 524   zUndefThreshold = zUndef + 10.
 525   if verbose:
 526     ligne = "zUndefThreshold: %f" % zUndefThreshold
 527     print (ligne)
 528 #
 529   #print "bathy :\n", bathy
 530   np_aux_z = (bathy < zUndefThreshold).nonzero()
 531   if verbose:
 532     ligne = ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
 533     print (ligne)
 534 #
 535 # 6. Option : xyz file
 536 #
 537   if xyzFile:
 538     if verbose:
 539       ligne = ".. Ecriture du champ de bathymétrie sur le fichier :\n%s" % fichierFonds
 540       print (ligne)
 541 #
 542     with open(fichierFonds, "w") as fo :
 543       for nodeId in range(nbnodes):
 544         ligne = "%11.3f %11.3f %11.3f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
 545         fo.write(ligne)
 546 #
 547 # 7. Final MED file
 548 # 7.1. Transformation of the bathymetry as a double array
 549 #
 550   bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
 551   bathy_dd.setInfoOnComponents(["Z [m]"])
 552 #
 553 # 7.2. If requested, modification of the z coordinate
 554 #
 555   if m3d:
 556     coords3D = ml.DataArrayDouble.Meld([coords[:,0:2], bathy_dd])
 557     coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
 558     #print "coords3D =\n", coords3D
 559     meshMEDFileRead.setCoords(coords3D)
 560 #
 561 # 7.3. Writes the mesh
 562 #
 563   if verbose:
 564     if m3d:
 565       saux = " 3D"
 566     else:
 567       saux = ""
 568     ligne = ".. Ecriture du maillage" + saux
 569     ligne += " sur le fichier :\n%s" % fichierFMaillage
 570     print (ligne)
 571 #
 572   meshMEDFileRead.write(fichierFMaillage, 2)
 573 #
 574 # 7.4. Writes the field
 575 #
 576   med_field_name = "BOTTOM"
 577   if verbose:
 578     ligne = ".. Ecriture du champ '%s'" % med_field_name
 579     print (ligne)
 580 #
 581   #print "bathy_dd =\n", bathy_dd
 582   fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES)
 583   fieldOnNodes.setName(med_field_name)
 584   fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0))
 585   fieldOnNodes.setArray(bathy_dd)
 586 # Ces valeurs d'instants sont mises pour assurer la lecture par TELEMAC
 587 # instant = 0.0
 588 # numero d'itération : 0
 589 # pas de numero d'ordre (-1)
 590   fieldOnNodes.setTime(0.0, 0, -1)
 591 #
 592   fMEDFile_ch_d = ml.MEDFileField1TS()
 593   fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
 594   fMEDFile_ch_d.write(fichierFMaillage, 0)
 595 #
 596   return statz