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