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=90., 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 90.
  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 -9000. It is tied with the value -9999. given
 277 #                     by the interpolation when no value is defined.
 278 #
 279     zUndefThreshold = -9000.
 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     if len(np_aux_z[0]):
 290       for iaux in np_aux_z[0]:
 291         bathy[iaux] = zUndef
 292 #
 293 # 6. Option : xyz file
 294 #
 295     if xyzFile:
 296 #
 297       if verbose:
 298         ligne = ".. Ecriture du champ de bathymétrie sur le fichier :\n%s" % fichierFonds
 299         print (ligne)
 300 #
 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])
 304           fo.write(ligne)
 305 #
 306 # 7. Final MED file
 307 # 7.1. Transformation of the bathymetry as a double array
 308 #
 309     bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
 310     bathy_dd.setInfoOnComponents(["Z [m]"])
 311 #
 312 # 7.2. If requested, modification of the z coordinate
 313 #
 314     if m3d:
 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)
 319 #
 320 # 7.3. Writes the mesh
 321 #
 322     if verbose:
 323       if m3d:
 324         saux = " 3D"
 325       else:
 326         saux = ""
 327       ligne = ".. Ecriture du maillage" + saux
 328       ligne += " sur le fichier :\n%s" % fichierFMaillage
 329       print (ligne)
 330 #
 331     meshMEDFileRead.write(fichierFMaillage, 2)
 332 #
 333 # 7.4. Writes the field
 334 #
 335     med_field_name = "BOTTOM"
 336     if verbose:
 337       ligne = ".. Ecriture du champ '%s'" % med_field_name
 338       print (ligne)
 339 #
 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
 346 #   instant = 0.0
 347 #   numero d'itération : 0
 348 #   pas de numero d'ordre (-1)
 349     fieldOnNodes.setTime(0.0, 0, -1)
 350 #
 351     fMEDFile_ch_d = ml.MEDFileField1TS()
 352     fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
 353     fMEDFile_ch_d.write(fichierFMaillage, 0)
 354 #
 355     break
 356 #
 357   if erreur:
 358     print(message)
 359 #
 360   return statz
 361
 362
 363 def interpolZ_B(bathyName, fichierMaillage, gr_face_name, zUndef=90., interp_method=0, m3d=True, xyzFile=False, verbose=False):
 364   """
 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.
 367
 368   In:
 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).
 373                      default value is 90.
 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.
 379   Out:
 380     if OK : statz_gr_face_name: statistic for z: (minz, maxz, meanz, stdz, v05z, v95z)
 381         else FALSE
 382   Out:
 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
 386   """
 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")
 392     return False
 393
 394   basename = fichierMaillage[:-4]
 395   fichierFMaillage = basename + 'F.med'
 396
 397   ligne = "fichierMaillage  = %s" % fichierMaillage
 398   ligne += "\nfichierFMaillage = %s" % fichierFMaillage
 399
 400   if xyzFile:
 401     fichierFonds = basename + '.xyz'
 402     ligne += "\nfichierFonds     = %s" % fichierFonds
 403   print (ligne)
 404 #
 405 # 1. Reads the mesh
 406 #
 407   meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
 408   if verbose:
 409     print (meshMEDFileRead)
 410 #
 411 # 2. Checks the names of the groups of faces
 412 #
 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")
 417     return False
 418 #
 419 # 3. Gets the information about the nodes
 420 #
 421   nbnodes = meshMEDFileRead.getNumberOfNodes()
 422   if verbose:
 423     ligne = "Number of nodes: %d" % nbnodes
 424     print (ligne)
 425 #
 426   coords = meshMEDFileRead.getCoords()
 427   #print(coords[:,2])
 428   if verbose:
 429     nb_comp = coords.getNumberOfComponents()
 430     l_info = coords.getInfoOnComponents()
 431     ligne = ""
 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):
 436         if l_info[iaux]:
 437           saux = l_info[iaux]
 438         else:
 439           saux = l_info_0[iaux]
 440         ligne += " %s" % saux
 441         ligne += "=%f" % coords[id_node, iaux]
 442     print (ligne)
 443 #
 444 # 4. Exploration of every group of faces
 445 #
 446   tb_aux = np.zeros(nbnodes, dtype=np.bool)
 447 #
 448   #bathy = np.zeros(nbnodes, dtype=np.float)
 449   #bathy.fill(zUndef)
 450   bathy = coords[:,2].toNumPyArray()
 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 -9000. It is tied with the value -9999. given
 521 #                     by the interpolation when no value is defined.
 522 #
 523   zUndefThreshold = -9000.
 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   if len(np_aux_z[0]):
 534     for iaux in np_aux_z[0]:
 535       bathy[iaux] = zUndef
 536 #
 537 # 6. Option : xyz file
 538 #
 539   if xyzFile:
 540     if verbose:
 541       ligne = ".. Ecriture du champ de bathymétrie sur le fichier :\n%s" % fichierFonds
 542       print (ligne)
 543 #
 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])
 547         fo.write(ligne)
 548 #
 549 # 7. Final MED file
 550 # 7.1. Transformation of the bathymetry as a double array
 551 #
 552   bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
 553   bathy_dd.setInfoOnComponents(["Z [m]"])
 554 #
 555 # 7.2. If requested, modification of the z coordinate
 556 #
 557   if m3d:
 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)
 562 #
 563 # 7.3. Writes the mesh
 564 #
 565   if verbose:
 566     if m3d:
 567       saux = " 3D"
 568     else:
 569       saux = ""
 570     ligne = ".. Ecriture du maillage" + saux
 571     ligne += " sur le fichier :\n%s" % fichierFMaillage
 572     print (ligne)
 573 #
 574   meshMEDFileRead.write(fichierFMaillage, 2)
 575 #
 576 # 7.4. Writes the field
 577 #
 578   med_field_name = "BOTTOM"
 579   if verbose:
 580     ligne = ".. Ecriture du champ '%s'" % med_field_name
 581     print (ligne)
 582 #
 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
 589 # instant = 0.0
 590 # numero d'itération : 0
 591 # pas de numero d'ordre (-1)
 592   fieldOnNodes.setTime(0.0, 0, -1)
 593 #
 594   fMEDFile_ch_d = ml.MEDFileField1TS()
 595   fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
 596   fMEDFile_ch_d.write(fichierFMaillage, 0)
 597 #
 598   return statz_gr_face_name