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, 0, 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   fieldOnNodes.setTime(0.0,0,-1)
 140   mm=MEDFileMesh.New(fichierZMaillage)
 141   mm.write(fichierFMaillage,2)
 142   MEDLoader.WriteFieldUsingAlreadyWrittenMesh(fichierFMaillage,fieldOnNodes)
 143   print fichierFMaillage, " field BOTTOM OK"
 144
 145 # -----------------------------------------------------------------------------
 146
 147 import HYDROPy
 148
 149 class MyInterpolator( HYDROPy.HYDROData_IInterpolator ):
 150   def GetAltitudeForPoint( self, theCoordX, theCoordY ):
 151     alt_obj = HYDROPy.HYDROData_IInterpolator.GetAltitudeObject( self );
 152     z = alt_obj.GetAltitudeForPoint( theCoordX, theCoordY )    # Custom calculation takes the base value and changes it to test
 153     #z2 = (z - 74.0)*10
 154     z2 = z
 155     return z2
 156
 157 # -----------------------------------------------------------------------------
 158
 159 import SMESH
 160 import SALOMEDS
 161 from salome.smesh import smeshBuilder
 162
 163 smesh = smeshBuilder.New(theStudy)
 164
 165 def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod = 0):
 166   """
 167   interpolZ takes a 2D (x,y) mesh and calls the active instance of module HYDRO
 168   to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node.
 169   interpolZ must be followed by createZfield1 or createZfield2.
 170   nomCas: Calculation Case Name in module HYDRO
 171   fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case
 172   dicoGroupeRegion: python dictionary giving the coorespondance of mesh groups to HYDRO regions.
 173                     Key: face group name, value: region name in the HYDRO Case
 174   zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
 175   interpolMethod: integer value, default 0 = nearest point on bathymetry, 1 = linear interpolation
 176   return <fichierMaillage>Z.med : med file with Z value on nodes
 177   return <fichierMaillage>F.med : an exact copy of <fichierMaillage>Z.med
 178   """
 179
 180   doc = HYDROPy.HYDROData_Document.Document( theStudyId )
 181   cas = doc.FindObjectByName(nomCas)
 182   print cas
 183   custom_inter = MyInterpolator()
 184
 185   noms = string.split(fichierMaillage,'.')
 186   basename = string.join(noms[:-1], '.')
 187   fichierZMaillage = basename + 'Z.med'
 188   fichierFMaillage = basename + 'F.med'
 189   fichierFonds = basename + '.xyz'
 190
 191   print fichierMaillage
 192   print fichierZMaillage
 193   print fichierFMaillage
 194   print fichierFonds
 195
 196   regions = {}
 197   ([maillagePlat], status) = smesh.CreateMeshesFromMED(fichierMaillage)
 198   groups = maillagePlat.GetGroups()
 199
 200   grpns = [grp for grp in groups if grp.GetType() == SMESH.NODE]
 201   if len(grpns) == 0:
 202     print "Problem! There are no groups of nodes in the mesh!"
 203     print "Please create at least the groups of nodes corresponding to each region of the HYDRO case"
 204     return {}
 205
 206
 207   for grp in groups:
 208     if grp.GetType() == SMESH.NODE:
 209       grpName = grp.GetName()
 210       print grpName
 211       if grpName in dicoGroupeRegion.keys():
 212         regions[dicoGroupeRegion[grpName]] = grp
 213
 214   fo = open(fichierFonds, 'w')
 215   statz = {}
 216   for nomreg, grp in regions.iteritems():
 217     print "------- region : ", nomreg
 218     region  = doc.FindObjectByName(nomreg)
 219     #print region
 220     #region.SetInterpolator(custom_inter)
 221     nodesIds = grp.GetListOfID()
 222     #print nodesIds
 223     vx = []
 224     vy = []
 225     for nodeId in nodesIds:
 226       xyz = maillagePlat.GetNodeXYZ(nodeId)
 227       #print xyz
 228       vx.append(xyz[0])
 229       vy.append(xyz[1])
 230     vz = cas.GetAltitudesForPoints( vx, vy, region, interpolMethod )
 231     minz = np.amin(vz)
 232     maxz = np.amax(vz)
 233     meanz = np.mean(vz)
 234     stdz = np.std(vz)
 235     v05z = np.percentile(vz,05)
 236     v95z = np.percentile(vz,95)
 237
 238     statz[grp.GetName()] = (minz, maxz, meanz, stdz, v05z, v95z)
 239
 240
 241     for i,nodeId in enumerate(nodesIds):
 242       x = vx[i]
 243       y = vy[i]
 244       z = vz[i]
 245       if z < -9000:
 246         z = zUndef
 247       #print i, nodeId, x, y, z
 248       maillagePlat.MoveNode(nodeId, x, y, z)
 249       l = "%10.2f %10.2f %10.2f\n"%(x, y, z)
 250       fo.write(l)
 251
 252   fo.close()
 253   maillagePlat.ExportMED(fichierZMaillage , 0, SMESH.MED_V2_2, 1 )
 254   maillagePlat.ExportMED(fichierFMaillage , 0, SMESH.MED_V2_2, 1 )
 255   return statz
 256
 257 # -----------------------------------------------------------------------------
 258