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 # -----------------------------------------------------------------------------
  27
  28 import string
  29
  30 # -----------------------------------------------------------------------------
  31
  32 import sys
  33 import salome
  34
  35 salome.salome_init()
  36 theStudy = salome.myStudy
  37 theStudyId = salome.myStudyId
  38
  39 import numpy as np
  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   noms = string.split(fichierMaillage,'.')
  64   basename = string.join(noms[:-1], '.')
  65   fichierFMaillage = basename + 'F.med'
  66   print fichierFMaillage
  67
  68   # --- ouverture du fichier
  69   fid=medfile.MEDfileOpen(fichierFMaillage, medenum.MED_ACC_RDEXT)
  70
  71   maa, sdim, mdim, meshtype, desc, dtunit, sort, nstep,  repere, axisname, axisunit = medmesh.MEDmeshInfo(fid, 1)
  72   print "Maillage de nom : %s , de dimension : %d , et de type %s"%(maa,mdim,meshtype)
  73   print "   Dimension de l'espace : %d"%(sdim)
  74
  75   # --- Combien de noeuds a lire ?
  76   nnoe, chgt, trsf = medmesh.MEDmeshnEntity(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT,
  77                                              medenum.MED_NODE, medenum.MED_NONE,
  78                                              medenum.MED_COORDINATE, medenum.MED_NO_CMODE)
  79
  80   if nnoe > 0:
  81     # --- Allocations memoire
  82     # --- table des coordonnees flt : (dimension * nombre de noeuds )
  83     coords = medfile.MEDFLOAT(nnoe*sdim)
  84     # --- table des numeros des noeuds
  85     numnoe = medfile.MEDINT(nnoe)
  86
  87     # --- Lecture des composantes des coordonnees des noeuds
  88     medmesh.MEDmeshNodeCoordinateRd(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT,
  89                                     medenum.MED_FULL_INTERLACE, coords)
  90     #print "Valeur de coords : ",coords
  91     valz=medfile.MEDFLOAT([z for (i,z) in enumerate(coords) if i%3==2])
  92     #print "Valeur de z : ",valz
  93
  94     # --- creation du champ
  95     nomcha1 = "BOTTOM"
  96     ncomp1 = 1
  97           # --1234567890123456--
  98     comp1  = "z               "
  99     unit1  = "m               "
 100     dtunit = ""
 101     medfield.MEDfieldCr(fid, nomcha1, medfile.MED_FLOAT64,
 102                         ncomp1, comp1, unit1, dtunit, maa)
 103
 104     # --- ecriture du champ
 105
 106     medfield.MEDfieldValueWr(fid, nomcha1, medenum.MED_NO_DT, medenum.MED_NO_IT, 0.0,
 107                              medenum.MED_NODE, medenum.MED_NONE,
 108                              medenum.MED_FULL_INTERLACE, medenum.MED_ALL_CONSTITUENT, nnoe, valz)
 109   # --- fermeture du fichier
 110   medfile.MEDfileClose(fid)
 111   print fichierFMaillage, " field BOTTOM OK"
 112
 113 # -----------------------------------------------------------------------------
 114
 115 import MEDLoader
 116 from MEDLoader import MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh
 117
 118 def createZfield2(fichierMaillage):
 119   """
 120   Complete the mesh for Telemac.
 121   Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes.
 122   createZfield2 is used after interpolZ. createZfield1 is base on MEDLoader interface.
 123   There is an alternate method based on Med file, equivalent (createZfield1).
 124   fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
 125   return <fichierMaillage>F.med : med file containing the field "BOTTOM"
 126   """
 127
 128   noms = string.split(fichierMaillage,'.')
 129   basename = string.join(noms[:-1], '.')
 130   fichierZMaillage = basename + 'Z.med'
 131   fichierFMaillage = basename + 'F.med'
 132   print fichierFMaillage
 133
 134   mymesh = MEDLoader.ReadUMeshFromFile(fichierZMaillage,0)
 135   fieldOnNodes=MEDCouplingFieldDouble.New(ON_NODES)
 136   fieldOnNodes.setName("BOTTOM")
 137   fieldOnNodes.setMesh(mymesh)
 138   fieldOnNodes.setArray(mymesh.getCoords()[:,2])
 139   mm=MEDFileMesh.New(fichierZMaillage)
 140   mm.write(fichierFMaillage,2)
 141   MEDLoader.WriteFieldUsingAlreadyWrittenMesh(fichierFMaillage,fieldOnNodes)
 142   print fichierFMaillage, " field BOTTOM OK"
 143
 144 # -----------------------------------------------------------------------------
 145
 146 import HYDROPy
 147
 148 class MyInterpolator( HYDROPy.HYDROData_IInterpolator ):
 149   def GetAltitudeForPoint( self, theCoordX, theCoordY ):
 150     alt_obj = HYDROPy.HYDROData_IInterpolator.GetAltitudeObject( self );
 151     z = alt_obj.GetAltitudeForPoint( theCoordX, theCoordY )    # Custom calculation takes the base value and changes it to test
 152     #z2 = (z - 74.0)*10
 153     z2 = z
 154     return z2
 155
 156 # -----------------------------------------------------------------------------
 157
 158 import SMESH
 159 import SALOMEDS
 160 from salome.smesh import smeshBuilder
 161
 162 smesh = smeshBuilder.New(theStudy)
 163
 164 def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod = 0):
 165   """
 166   interpolZ takes a 2D (x,y) mesh and calls the active instance of module HYDRO
 167   to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node.
 168   interpolZ must be followed by createZfield1 or createZfield2.
 169   nomCas: Calculation Case Name in module HYDRO
 170   fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case
 171   dicoGroupeRegion: python dictionary giving the coorespondance of mesh groups to HYDRO regions.
 172                     Key: face group name, value: region name in the HYDRO Case
 173   zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
 174   interpolMethod: integer value, default 0 = nearest point on bathymetry, 1 = linear interpolation
 175   return <fichierMaillage>Z.med : med file with Z value on nodes
 176   return <fichierMaillage>F.med : an exact copy of <fichierMaillage>Z.med
 177   """
 178
 179   doc = HYDROPy.HYDROData_Document.Document( theStudyId )
 180   cas = doc.FindObjectByName(nomCas)
 181   print cas
 182   custom_inter = MyInterpolator()
 183
 184   noms = string.split(fichierMaillage,'.')
 185   basename = string.join(noms[:-1], '.')
 186   fichierZMaillage = basename + 'Z.med'
 187   fichierFMaillage = basename + 'F.med'
 188   fichierFonds = basename + '.xyz'
 189
 190   print fichierMaillage
 191   print fichierZMaillage
 192   print fichierFMaillage
 193   print fichierFonds
 194
 195   regions = {}
 196   ([maillagePlat], status) = smesh.CreateMeshesFromMED(fichierMaillage)
 197   groups = maillagePlat.GetGroups()
 198
 199   grpns = [grp for grp in groups if grp.GetType() == SMESH.NODE]
 200   if len(grpns) == 0:
 201     print "Problem! There are no groups of nodes in the mesh!"
 202     print "Please create at least the groups of nodes corresponding to each region of the HYDRO case"
 203     return {}
 204
 205
 206   for grp in groups:
 207     if grp.GetType() == SMESH.NODE:
 208       grpName = grp.GetName()
 209       print grpName
 210       if grpName in dicoGroupeRegion.keys():
 211         regions[dicoGroupeRegion[grpName]] = grp
 212
 213   fo = open(fichierFonds, 'w')
 214   statz = {}
 215   for nomreg, grp in regions.iteritems():
 216     print "------- region : ", nomreg
 217     region  = doc.FindObjectByName(nomreg)
 218     #print region
 219     #region.SetInterpolator(custom_inter)
 220     nodesIds = grp.GetListOfID()
 221     #print nodesIds
 222     vx = []
 223     vy = []
 224     for nodeId in nodesIds:
 225       xyz = maillagePlat.GetNodeXYZ(nodeId)
 226       #print xyz
 227       vx.append(xyz[0])
 228       vy.append(xyz[1])
 229     vz = cas.GetAltitudesForPoints( vx, vy, region, interpolMethod )
 230     minz = np.amin(vz)
 231     maxz = np.amax(vz)
 232     meanz = np.mean(vz)
 233     stdz = np.std(vz)
 234     v05z = np.percentile(vz,05)
 235     v95z = np.percentile(vz,95)
 236
 237     statz[grp.GetName()] = (minz, maxz, meanz, stdz, v05z, v95z)
 238
 239
 240     for i,nodeId in enumerate(nodesIds):
 241       x = vx[i]
 242       y = vy[i]
 243       z = vz[i]
 244       if z < -9000:
 245         z = zUndef
 246       #print i, nodeId, x, y, z
 247       maillagePlat.MoveNode(nodeId, x, y, z)
 248       l = "%10.2f %10.2f %10.2f\n"%(x, y, z)
 249       fo.write(l)
 250
 251   fo.close()
 252   maillagePlat.ExportMED(fichierZMaillage , 0, SMESH.MED_V2_2, 1 )
 253   maillagePlat.ExportMED(fichierFMaillage , 0, SMESH.MED_V2_2, 1 )
 254   return statz
 255
 256 # -----------------------------------------------------------------------------
 257