src/HYDROTools/interpolZ.py

   1 # -*- coding: utf-8 -*-
   2
   3 """
   4 example of use case of interpolZ:
   5
   6 # --- case name in HYDRO
   7 nomCas = 'inondation1'
   8
   9 # --- med file 2D(x,y) of the case produced by SMESH
  10 fichierMaillage = '/home/B27118/projets/LNHE/garonne/inondation1.med'
  11
  12 # --- dictionary [med group name] = region name
  13 dicoGroupeRegion= dict(litMineur          = 'inondation1_litMineur',
  14                        litMajeurDroite    = 'inondation1_riveDroite',
  15                        litMajeurGauche    = 'inondation1_riveGauche',
  16                        )
  17 # --- value to use for Z when the node is not in a region (used to detect problems)
  18 zUndef = 90
  19
  20 # --- Z interpolation on the bathymety/altimetry on the mesh nodes
  21 interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod)
  22 """
  23 __revision__ = "V2.04"
  24 # -----------------------------------------------------------------------------
  25
  26 # -----------------------------------------------------------------------------
  27
  28 import salome
  29
  30 salome.salome_init()
  31 theStudy = salome.myStudy
  32 theStudyId = salome.myStudyId
  33
  34 import numpy as np
  35 import MEDLoader as ml
  36 import MEDCoupling as mc
  37
  38 # -----------------------------------------------------------------------------
  39
  40 from med import medfile
  41 from med import medmesh
  42 #from med import medfilter
  43 from med import medfield
  44 from med import medenum
  45 #from med import medprofile
  46 #from med import medlocalization
  47 #from med import medlink
  48
  49 # La fonction createZfield1 ne sert plus à partir du 12/07/2017
  50 def createZfield1(fichierMaillage):
  51   """
  52   Complete the mesh for Telemac.
  53   Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes.
  54   createZfield1 is used after interpolZ. createZfield1 is base on med file interface.
  55   There is an alternate method based on MEDLoader, equivalent (createZfield2).
  56   The file <fichierMaillage>F.med produced by interpolz must exist, and is modified.
  57   fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
  58   return <fichierMaillage>F.med : med file containing the field "BOTTOM"
  59   """
  60
  61   basename = fichierMaillage[:-4]
  62   fichierFMaillage = basename + 'F.med'
  63   print fichierFMaillage
  64
  65   # --- ouverture du fichier
  66   fid=medfile.MEDfileOpen(fichierFMaillage, medenum.MED_ACC_RDEXT)
  67
  68   maa, sdim, mdim, meshtype, desc, dtunit, sort, nstep,  repere, axisname, axisunit = medmesh.MEDmeshInfo(fid, 1)
  69   print "Maillage de nom : %s , de dimension : %d , et de type %s"%(maa,mdim,meshtype)
  70   print "   Dimension de l'espace : %d"%(sdim)
  71
  72   # --- Combien de noeuds a lire ?
  73   nnoe, chgt, trsf = medmesh.MEDmeshnEntity(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT,
  74                                             medenum.MED_NODE, medenum.MED_NONE,
  75                                             medenum.MED_COORDINATE, medenum.MED_NO_CMODE)
  76
  77   if nnoe > 0:
  78     # --- Allocations memoire
  79     # --- table des coordonnees flt : (dimension * nombre de noeuds )
  80     coords = medfile.MEDFLOAT(nnoe*sdim)
  81     # --- table des numeros des noeuds
  82     numnoe = medfile.MEDINT(nnoe)
  83
  84     # --- Lecture des composantes des coordonnees des noeuds
  85     medmesh.MEDmeshNodeCoordinateRd(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT,
  86                                     medenum.MED_FULL_INTERLACE, coords)
  87     #print "Valeur de coords : ",coords
  88     valz=medfile.MEDFLOAT([z for (i,z) in enumerate(coords) if i%3==2])
  89     #print "Valeur de z : ",valz
  90
  91     # --- creation du champ
  92     nomcha1 = "BOTTOM"
  93     ncomp1 = 1
  94           # --1234567890123456--
  95     comp1  = "z               "
  96     unit1  = "m               "
  97     dtunit = ""
  98     medfield.MEDfieldCr(fid, nomcha1, medfile.MED_FLOAT64,
  99                         ncomp1, comp1, unit1, dtunit, maa)
 100
 101     # --- ecriture du champ
 102
 103     medfield.MEDfieldValueWr(fid, nomcha1, 0, medenum.MED_NO_IT, 0.0,
 104                              medenum.MED_NODE, medenum.MED_NONE,
 105                              medenum.MED_FULL_INTERLACE, medenum.MED_ALL_CONSTITUENT, nnoe, valz)
 106   # --- fermeture du fichier
 107   medfile.MEDfileClose(fid)
 108   print fichierFMaillage, " field BOTTOM OK"
 109
 110 # -----------------------------------------------------------------------------
 111
 112 # La fonction createZfield2 ne sert plus à partir du 12/07/2017
 113 def createZfield2(fichierMaillage):
 114   """
 115   Complete the mesh for Telemac.
 116   Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes.
 117   createZfield2 is used after interpolZ. createZfield1 is base on MEDLoader interface.
 118   There is an alternate method based on Med file, equivalent (createZfield1).
 119   fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
 120   return <fichierMaillage>F.med : med file containing the field "BOTTOM"
 121   """
 122
 123   import MEDLoader
 124   from MEDLoader import MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh
 125
 126   basename = fichierMaillage[:-4]
 127   fichierZMaillage = basename + 'Z.med'
 128   fichierFMaillage = basename + 'F.med'
 129   print fichierFMaillage
 130
 131   mymesh = MEDLoader.ReadUMeshFromFile(fichierZMaillage,0)
 132   fieldOnNodes=MEDCouplingFieldDouble.New(ON_NODES)
 133   fieldOnNodes.setName("BOTTOM")
 134   fieldOnNodes.setMesh(mymesh)
 135   fieldOnNodes.setArray(mymesh.getCoords()[:,2])
 136   fieldOnNodes.setTime(0.0,0,-1)
 137   mm=MEDFileMesh.New(fichierZMaillage)
 138   mm.write(fichierFMaillage,2)
 139   MEDLoader.WriteFieldUsingAlreadyWrittenMesh(fichierFMaillage,fieldOnNodes)
 140   print fichierFMaillage, " field BOTTOM OK"
 141
 142 # -----------------------------------------------------------------------------
 143
 144 import HYDROPy
 145
 146 class MyInterpolator( HYDROPy.HYDROData_IInterpolator ):
 147   def GetAltitudeForPoint( self, theCoordX, theCoordY ):
 148     alt_obj = HYDROPy.HYDROData_IInterpolator.GetAltitudeObject( self );
 149     z = alt_obj.GetAltitudeForPoint( theCoordX, theCoordY )    # Custom calculation takes the base value and changes it to test
 150     #z2 = (z - 74.0)*10
 151     z2 = z
 152     return z2
 153
 154 # -----------------------------------------------------------------------------
 155
 156
 157 def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMethod=0, xyzFile=False, verbose=False):
 158   """
 159   interpolZ takes a 2D (x,y) mesh and calls the active instance of module HYDRO
 160   to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node.
 161
 162   In:
 163     nomCas: Calculation Case Name in module HYDRO
 164     fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case
 165     dicoGroupeRegion: python dictionary giving the correspondance of mesh groups to HYDRO regions.
 166                       Key: face group name
 167                       Value: region name in the HYDRO Case
 168     zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
 169                      default value is 90.
 170     interpolMethod: integer value
 171                     0 = nearest point on bathymetry (default)
 172                     1 = linear interpolation
 173     zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
 174     xyzFile: True/False to write an ascii file with xyz for every node. Default is False.
 175   Out:
 176     statz: statistique for z
 177                       Key: face group name
 178                       Value: (minz, maxz, meanz, stdz, v05z, v95z)
 179   Out:
 180     return <fichierMaillage>F.med : med file with Z value on nodes and in a field "BOTTOM"
 181   """
 182   statz = dict()
 183   erreur = 0
 184   message = ""
 185
 186   while not erreur:
 187
 188     if verbose:
 189       print "nomCas:", nomCas
 190       print "interpolMethod: %d" % interpolMethod
 191       print "zUndef:", zUndef
 192
 193     doc = HYDROPy.HYDROData_Document.Document(theStudyId)
 194     cas = doc.FindObjectByName(nomCas)
 195     print cas
 196     custom_inter = MyInterpolator()
 197
 198     basename = fichierMaillage[:-4]
 199     fichierFMaillage = basename + 'F.med'
 200
 201     print "dicoGroupeRegion =", dicoGroupeRegion
 202     print "fichierMaillage  =", fichierMaillage
 203     print "fichierFMaillage =", fichierFMaillage
 204     if xyzFile:
 205       fichierFonds = basename + '.xyz'
 206       print "fichierFonds     =", fichierFonds
 207 #
 208 # 1. Reads the mesh
 209 #
 210     meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
 211 #
 212 # 2. Checks the names of the groups of faces
 213 #
 214     l_gr_faces = list(dicoGroupeRegion.keys())
 215     t_groupe_n = meshMEDFileRead.getGroupsNames()
 216     nb_pb = 0
 217     for gr_face_name in l_gr_faces:
 218       if gr_face_name not in t_groupe_n:
 219         message += "Group: '" + gr_face_name + "'\n"
 220         nb_pb += 1
 221     if verbose:
 222       print "Number of problems: %d" % nb_pb
 223 #
 224     if nb_pb > 0:
 225       if nb_pb == 1:
 226         message += "This group does"
 227       else:
 228         message += "That %d groups do" % nb_pb
 229       message += " not belongs to the mesh.\n"
 230       message += "Please check the names of the group(s) of faces corresponding to each region of the HYDRO case"
 231       erreur = 2
 232       break
 233 #
 234 # 3. Gets the information about the nodes
 235 #
 236     nbnodes = meshMEDFileRead.getNumberOfNodes()
 237     if verbose:
 238       print "Number of nodes: %d" % nbnodes
 239 #
 240     coords = meshMEDFileRead.getCoords()
 241     #print coords
 242     #print coords[0,0]
 243     #print coords[0,1]
 244 #
 245 # 4. Exploration of every group of faces
 246 #
 247     tb_aux = np.zeros(nbnodes, dtype=np.bool)
 248 #
 249     bathy = np.zeros(nbnodes, dtype=np.float)
 250     bathy.fill(zUndef)
 251 #
 252     for gr_face_name in l_gr_faces:
 253 #
 254 #     4.1. Region connected to the group
 255 #
 256       nomreg = dicoGroupeRegion[gr_face_name]
 257       line = "------- Region: '" + nomreg + "'"
 258       line += ", connected to group '" + gr_face_name + "'"
 259       print line
 260       region = doc.FindObjectByName(nomreg)
 261       #print region
 262       #region.SetInterpolator(custom_inter)
 263 #
 264 #     4.2. Mesh of the group
 265 #
 266       mesh_of_the_group = meshMEDFileRead.getGroup(0, gr_face_name, False)
 267       nbr_cells = mesh_of_the_group.getNumberOfCells()
 268       if verbose:
 269         print "\t. Number of cells: %d" % nbr_cells
 270 #
 271 #     4.3. Nodes of the meshes of the group
 272 #          Every node is flagged in tb_aux
 273 #
 274       tb_aux.fill(False)
 275       for id_elem in range(nbr_cells):
 276         l_nodes = mesh_of_the_group.getNodeIdsOfCell(id_elem)
 277         #print l_nodes
 278         for id_node in l_nodes:
 279           tb_aux[id_node] = True
 280       np_aux = tb_aux.nonzero()
 281       if len(np_aux[0]):
 282         if verbose:
 283           print "\t. Number of nodes for this group: %d" % len(np_aux[0])
 284       #print "np_aux:", np_aux
 285 #
 286 #     4.4. Interpolation over the nodes of the meshes of the group
 287 #
 288       vx = []
 289       vy = []
 290       for nodeId in np_aux[0]:
 291         vx.append(coords[nodeId, 0])
 292         vy.append(coords[nodeId, 1])
 293       #print "vx:\n", vx
 294       #print "vy:\n", vy
 295       vz = cas.GetAltitudesForPoints(vx, vy, region, interpolMethod)
 296       #print "vz:\n", vz
 297       minz = np.amin(vz)
 298       maxz = np.amax(vz)
 299       meanz = np.mean(vz)
 300       stdz = np.std(vz)
 301       v05z = np.percentile(vz, 05)
 302       v95z = np.percentile(vz, 95)
 303       if verbose:
 304         ligne = ".. Minimum: %f" % minz
 305         ligne += ", maximum: %f" % maxz
 306         ligne += ", mean: %f\n" % meanz
 307         ligne += ".. stdeviation: %f" % stdz
 308         ligne += ", v05z: %f" % v05z
 309         ligne += ", v95z: %f" % v95z
 310         print ligne
 311 #
 312 #     4.5. Storage of the z and of the statistics for this region
 313 #
 314       statz[gr_face_name] = (minz, maxz, meanz, stdz, v05z, v95z)
 315 #
 316       for iaux, nodeId in enumerate(np_aux[0]):
 317         bathy[nodeId] = vz[iaux]
 318 #
 319 # 5. Minimum:
 320 #    During the interpolation, if no value is available over a node, a default value
 321 #    is set: -9999. It has no importance for the final computation, but if the field
 322 #    or the mesh is displayed, it makes huge gap. To prevent this artefact, a more
 323 #    convenient "undefined" value is set. This new undefinde value is given by the user.
 324 #
 325 #    zUndefThreshold: the default value is -9000. It is tied with the value -9999. given
 326 #                     by the interpolation when no value is defined.
 327 #
 328     zUndefThreshold = -9000.
 329     if verbose:
 330       print "zUndefThreshold: %f" % zUndefThreshold#
 331 #
 332     #print "bathy :\n", bathy
 333     np_aux_z = (bathy < zUndefThreshold).nonzero()
 334     if verbose:
 335       print ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
 336     if len(np_aux_z[0]):
 337       for iaux in np_aux_z[0]:
 338         bathy[iaux] = zUndef
 339 #
 340 # 6. xyz file
 341 #
 342     if xyzFile:
 343 #
 344       if verbose:
 345         print ".. Ecriture du champ de bathymétrie sur le fichier :\n", fichierFonds
 346       fo = open(fichierFonds, 'w')
 347       for nodeId in range(nbnodes):
 348         line = "%10.2f %10.2f %10.2f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
 349         fo.write(line)
 350       fo.close()
 351 #
 352 # 7. Final MED file
 353 # 7.1. Modification of the z coordinates
 354 #
 355     bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
 356     bathy_dd.setInfoOnComponents(["Z [m]"])
 357 #
 358     coords3D = ml.DataArrayDouble.Meld([coords, bathy_dd])
 359     coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
 360     #print "coords3D =\n", coords3D
 361 #
 362     meshMEDFileRead.setCoords(coords3D)
 363 #
 364 # 7.2. Writes the 3D mesh
 365 #
 366     if verbose:
 367       print ".. Ecriture du maillage 3D sur le fichier :\n", fichierFMaillage
 368     meshMEDFileRead.write(fichierFMaillage, 2)
 369 #
 370 # 7.3. Writes the field
 371 #
 372     med_field_name = "BOTTOM"
 373     if verbose:
 374       print ".. Ecriture du champ '"+med_field_name+"'"
 375     #print "bathy_dd =\n", bathy_dd
 376     fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES)
 377     fieldOnNodes.setName(med_field_name)
 378     fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0))
 379     fieldOnNodes.setArray(bathy_dd)
 380 #   Ces valeurs d'instants sont mises pour assurer la lecture par TELEMAC
 381 #   instant = 0.0
 382 #   numero d'itération : 0
 383 #   pas de numero d'ordre (-1)
 384     fieldOnNodes.setTime(0.0, 0, -1)
 385 #
 386     fMEDFile_ch_d = ml.MEDFileField1TS()
 387     fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
 388     fMEDFile_ch_d.write(fichierFMaillage, 0)
 389 #
 390     break
 391
 392   if erreur:
 393     print message
 394
 395   return statz
 396