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., interpolMethod=0, m3d=False, 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
  80                     0 = nearest point on bathymetry (default)
  81                     1 = linear interpolation
  82     m3d: True/False to produce a 3D mesh. Default is False.
  83     xyzFile: True/False to write an ascii file with xyz for every node. Default is False.
  84   Out:
  85     statz: statistique for z
  86                       Key: face group name
  87                       Value: (minz, maxz, meanz, stdz, v05z, v95z)
  88   Out:
  89     return <fichierMaillage>F.med : med file with Z value in a field "BOTTOM"
  90                                     Option: Z value also in Z coordinate if m3D is true
  91     return <fichierMaillage>.xyz : text file with X, Y, Z values
  92   """
  93   statz = dict()
  94   erreur = 0
  95   message = ""
  96
  97   while not erreur:
  98
  99     if verbose:
 100       ligne = "nomCas: %s" % nomCas
 101       ligne += "\ninterpolMethod: %d" % interpolMethod
 102       if (zUndef != None ):
 103         ligne += "\nzUndef: %f" % zUndef
 104       ligne += "\nm3d: %d" % m3d
 105       print (ligne)
 106
 107     doc = HYDROPy.HYDROData_Document.Document(theStudyId)
 108     cas = doc.FindObjectByName(nomCas)
 109     print ( "cas : ", cas)
 110     custom_inter = MyInterpolator()
 111
 112     basename = fichierMaillage[:-4]
 113     fichierFMaillage = basename + 'F.med'
 114
 115     print ("dicoGroupeRegion = ", dicoGroupeRegion)
 116     ligne = "fichierMaillage  = %s" % fichierMaillage
 117     ligne += "\nfichierFMaillage = %s" % fichierFMaillage
 118     if xyzFile:
 119       fichierFonds = basename + '.xyz'
 120       ligne += "\nfichierFonds     = %s" % fichierFonds
 121     print (ligne)
 122 #
 123 # 1. Reads the mesh
 124 #
 125     meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
 126     if verbose:
 127       print (meshMEDFileRead)
 128 #
 129 # 2. Checks the names of the groups of faces
 130 #
 131     t_group_n = meshMEDFileRead.getGroupsNames()
 132     dicoGroupeRegion_0 = dict()
 133     nb_pb = 0
 134     for gr_face_name in dicoGroupeRegion:
 135       saux = gr_face_name.strip()
 136       if saux not in t_group_n:
 137         message += "Group: '" + gr_face_name + "'\n"
 138         nb_pb += 1
 139       else :
 140         dicoGroupeRegion_0[saux] =  dicoGroupeRegion[gr_face_name]
 141     if verbose:
 142       ligne = "Number of problems: %d" % nb_pb
 143       print (ligne)
 144 #
 145     if nb_pb > 0:
 146       if nb_pb == 1:
 147         message += "This group does"
 148       else:
 149         message += "These %d groups do" % nb_pb
 150       message += " not belong to the mesh.\n"
 151       message += "Please check the names of the group(s) of faces corresponding to each region of the HYDRO case.\n"
 152       message += "Groups for this file:\n"
 153       for group_n in t_group_n :
 154         message += "'%s'\n" % group_n
 155       erreur = 2
 156       break
 157 #
 158 # 3. Gets the information about the nodes
 159 #
 160     nbnodes = meshMEDFileRead.getNumberOfNodes()
 161     if verbose:
 162       ligne = "Number of nodes: %d" % nbnodes
 163       print (ligne)
 164 #
 165     coords = meshMEDFileRead.getCoords()
 166     #print coords
 167     if verbose:
 168       nb_comp = coords.getNumberOfComponents()
 169       l_info = coords.getInfoOnComponents()
 170       ligne = ""
 171       l_info_0=["X", "Y", "Z"]
 172       for id_node in (0, 1, nbnodes-1):
 173         ligne += "\nNode #%6d:" % id_node
 174         for iaux in range(nb_comp):
 175           if l_info[iaux]:
 176             saux = l_info[iaux]
 177           else:
 178             saux = l_info_0[iaux]
 179           ligne += " %s" % saux
 180           ligne += "=%f" % coords[id_node, iaux]
 181       print (ligne)
 182 #
 183 # 4. Exploration of every group of faces
 184 #
 185     tb_aux = np.zeros(nbnodes, dtype=np.bool)
 186 #
 187     bathy = np.zeros(nbnodes, dtype=np.float)
 188     bathy.fill(zUndef)
 189 #
 190     for gr_face_name in dicoGroupeRegion_0:
 191 #
 192 #     4.1. Region connected to the group
 193 #
 194       nomreg = dicoGroupeRegion_0[gr_face_name]
 195       ligne = "------- Region: '%s'" % nomreg
 196       ligne += ", connected to group '%s'" % gr_face_name
 197       print (ligne)
 198       region = doc.FindObjectByName(nomreg)
 199 #
 200 #     4.2. Mesh of the group
 201 #
 202       mesh_of_the_group = meshMEDFileRead.getGroup(0, gr_face_name, False)
 203       nbr_cells = mesh_of_the_group.getNumberOfCells()
 204       if verbose:
 205         ligne = "\t. Number of cells: %d" % nbr_cells
 206         print (ligne)
 207 #
 208 #     4.3. Nodes of the meshes of the group
 209 #          Every node is flagged in tb_aux
 210 #
 211       tb_aux.fill(False)
 212       for id_elem in range(nbr_cells):
 213         l_nodes = mesh_of_the_group.getNodeIdsOfCell(id_elem)
 214         #print l_nodes
 215         for id_node in l_nodes:
 216           tb_aux[id_node] = True
 217       np_aux = tb_aux.nonzero()
 218       if len(np_aux[0]):
 219         if verbose:
 220           ligne = "\t. Number of nodes for this group: %d" % len(np_aux[0])
 221           print (ligne)
 222       #print ("np_aux:", np_aux)
 223 #
 224 #     4.4. Interpolation over the nodes of the meshes of the group
 225 #
 226       vx = list()
 227       vy = list()
 228       for nodeId in np_aux[0]:
 229         vx.append(coords[nodeId, 0])
 230         vy.append(coords[nodeId, 1])
 231       #print ("vx:\n", vx)
 232       #print ("vy:\n", vy)
 233 #
 234       vz = cas.GetAltitudesForPoints(vx, vy, region, interpolMethod)
 235 #
 236       #print ("vz:\n", vz)
 237       minz = np.amin(vz)
 238       maxz = np.amax(vz)
 239       meanz = np.mean(vz)
 240       stdz = np.std(vz)
 241       v05z = np.percentile(vz, 05)
 242       v95z = np.percentile(vz, 95)
 243 #
 244       if verbose:
 245         ligne = ".. Minimum: %f" % minz
 246         ligne += ", maximum: %f" % maxz
 247         ligne += ", mean: %f\n" % meanz
 248         ligne += ".. stdeviation: %f" % stdz
 249         ligne += ", v05z: %f" % v05z
 250         ligne += ", v95z: %f" % v95z
 251         print (ligne)
 252 #
 253 #     4.5. Storage of the z and of the statistics for this region
 254 #
 255       statz[gr_face_name] = (minz, maxz, meanz, stdz, v05z, v95z)
 256 #
 257       for iaux, nodeId in enumerate(np_aux[0]):
 258         bathy[nodeId] = vz[iaux]
 259 #
 260 # 5. Minimum:
 261 #    During the interpolation, if no value is available over a node, a default value
 262 #    is set: -9999. It has no importance for the final computation, but if the field
 263 #    or the mesh is displayed, it makes huge gap. To prevent this artefact, a more
 264 #    convenient "undefined" value is set. This new undefined value is given by the user.
 265 #
 266 #    zUndefThreshold: the default value is -9000. It is tied with the value -9999. given
 267 #                     by the interpolation when no value is defined.
 268 #
 269     zUndefThreshold = -9000.
 270     if verbose:
 271       ligne = "zUndefThreshold: %f" % zUndefThreshold
 272       print (ligne)
 273 #
 274     #print "bathy :\n", bathy
 275     np_aux_z = (bathy < zUndefThreshold).nonzero()
 276     if verbose:
 277       ligne = ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
 278       print (ligne)
 279     if len(np_aux_z[0]):
 280       for iaux in np_aux_z[0]:
 281         bathy[iaux] = zUndef
 282 #
 283 # 6. Option : xyz file
 284 #
 285     if xyzFile:
 286 #
 287       if verbose:
 288         ligne = ".. Ecriture du champ de bathymétrie sur le fichier :\n%s" % fichierFonds
 289         print (ligne)
 290 #
 291       with open(fichierFonds, "w") as fo :
 292         for nodeId in range(nbnodes):
 293           ligne = "%11.3f %11.3f %11.3f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
 294           fo.write(ligne)
 295 #
 296 # 7. Final MED file
 297 # 7.1. Transformation of the bathymetry as a double array
 298 #
 299     bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
 300     bathy_dd.setInfoOnComponents(["Z [m]"])
 301 #
 302 # 7.2. If requested, modification of the z coordinate
 303 #
 304     if m3d:
 305       coords3D = ml.DataArrayDouble.Meld([coords[:,0:2], bathy_dd])
 306       coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
 307       #print "coords3D =\n", coords3D
 308       meshMEDFileRead.setCoords(coords3D)
 309 #
 310 # 7.3. Writes the mesh
 311 #
 312     if verbose:
 313       if m3d:
 314         saux = " 3D"
 315       else:
 316         saux = ""
 317       ligne = ".. Ecriture du maillage" + saux
 318       ligne += " sur le fichier :\n%s" % fichierFMaillage
 319       print (ligne)
 320 #
 321     meshMEDFileRead.write(fichierFMaillage, 2)
 322 #
 323 # 7.4. Writes the field
 324 #
 325     med_field_name = "BOTTOM"
 326     if verbose:
 327       ligne = ".. Ecriture du champ '%s'" % med_field_name
 328       print (ligne)
 329 #
 330     #print "bathy_dd =\n", bathy_dd
 331     fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES)
 332     fieldOnNodes.setName(med_field_name)
 333     fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0))
 334     fieldOnNodes.setArray(bathy_dd)
 335 #   Ces valeurs d'instants sont mises pour assurer la lecture par TELEMAC
 336 #   instant = 0.0
 337 #   numero d'itération : 0
 338 #   pas de numero d'ordre (-1)
 339     fieldOnNodes.setTime(0.0, 0, -1)
 340 #
 341     fMEDFile_ch_d = ml.MEDFileField1TS()
 342     fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
 343     fMEDFile_ch_d.write(fichierFMaillage, 0)
 344 #
 345     break
 346 #
 347   if erreur:
 348     print message
 349 #
 350   return statz
 351