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