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