src/HYDROTools/interpolZ.py

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