Merge remote-tracking branch 'origin/pre/V8_2_BR' into BR_PORTING_OCCT_7

[modules/hydro.git] / src / HYDROTools / interpolZ.py

diff --git a/src/HYDROTools/interpolZ.py b/src/HYDROTools/interpolZ.py
index 0321ca9db35f17528012804e09d9cad8567cf4f8..bf6a48eedfea40fb08e4d407d76f44a3e94a0c42 100644 (file)
--- a/src/HYDROTools/interpolZ.py
+++ b/src/HYDROTools/interpolZ.py
@@ -14,7 +14,7 @@ dicoGroupeRegion= dict(litMineur          = 'inondation1_litMineur',
                        litMajeurDroite    = 'inondation1_riveDroite',
                        litMajeurGauche    = 'inondation1_riveGauche',
                        )
-# --- value to use for Z when the node is not in a region (used to detect problems)                       
+# --- value to use for Z when the node is not in a region (used to detect problems)
 zUndef = 90
 
 # --- Z interpolation on the bathymety/altimetry on the mesh nodes
@@ -36,6 +36,8 @@ salome.salome_init()
 theStudy = salome.myStudy
 theStudyId = salome.myStudyId
 
+import numpy as np
+
 # -----------------------------------------------------------------------------
 
 from med import medfile
@@ -57,7 +59,7 @@ def createZfield1(fichierMaillage):
   fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
   return <fichierMaillage>F.med : med file containing the field "BOTTOM"
   """
-  
+
   noms = string.split(fichierMaillage,'.')
   basename = string.join(noms[:-1], '.')
   fichierFMaillage = basename + 'F.med'
@@ -82,21 +84,21 @@ def createZfield1(fichierMaillage):
     # --- table des numeros des noeuds
     numnoe = medfile.MEDINT(nnoe)
 
-    # --- Lecture des composantes des coordonnees des noeuds 
+    # --- Lecture des composantes des coordonnees des noeuds
     medmesh.MEDmeshNodeCoordinateRd(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT,
                                     medenum.MED_FULL_INTERLACE, coords)
     #print "Valeur de coords : ",coords
     valz=medfile.MEDFLOAT([z for (i,z) in enumerate(coords) if i%3==2])
     #print "Valeur de z : ",valz
 
-    # --- creation du champ 
+    # --- creation du champ
     nomcha1 = "BOTTOM"
     ncomp1 = 1
           # --1234567890123456--
     comp1  = "z               "
     unit1  = "m               "
     dtunit = ""
-    medfield.MEDfieldCr(fid, nomcha1, medfile.MED_FLOAT64, 
+    medfield.MEDfieldCr(fid, nomcha1, medfile.MED_FLOAT64,
                         ncomp1, comp1, unit1, dtunit, maa)
 
     # --- ecriture du champ
@@ -106,6 +108,7 @@ def createZfield1(fichierMaillage):
                              medenum.MED_FULL_INTERLACE, medenum.MED_ALL_CONSTITUENT, nnoe, valz)
   # --- fermeture du fichier
   medfile.MEDfileClose(fid)
+  print fichierFMaillage, " field BOTTOM OK"
 
 # -----------------------------------------------------------------------------
 
@@ -121,7 +124,7 @@ def createZfield2(fichierMaillage):
   fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
   return <fichierMaillage>F.med : med file containing the field "BOTTOM"
   """
-  
+
   noms = string.split(fichierMaillage,'.')
   basename = string.join(noms[:-1], '.')
   fichierZMaillage = basename + 'Z.med'
@@ -136,6 +139,7 @@ def createZfield2(fichierMaillage):
   mm=MEDFileMesh.New(fichierZMaillage)
   mm.write(fichierFMaillage,2)
   MEDLoader.WriteFieldUsingAlreadyWrittenMesh(fichierFMaillage,fieldOnNodes)
+  print fichierFMaillage, " field BOTTOM OK"
 
 # -----------------------------------------------------------------------------
 
@@ -148,7 +152,7 @@ class MyInterpolator( HYDROPy.HYDROData_IInterpolator ):
     #z2 = (z - 74.0)*10
     z2 = z
     return z2
-    
+
 # -----------------------------------------------------------------------------
 
 import SMESH
@@ -187,25 +191,25 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod
   print fichierZMaillage
   print fichierFMaillage
   print fichierFonds
-  
+
   regions = {}
   ([maillagePlat], status) = smesh.CreateMeshesFromMED(fichierMaillage)
   groups = maillagePlat.GetGroups()
-  
+
   grpns = [grp for grp in groups if grp.GetType() == SMESH.NODE]
   if len(grpns) == 0:
     print "Problem! There are no groups of nodes in the mesh!"
-    print "Please create at least the groups of nodes corresponding to each region of the HYDRO case" 
+    print "Please create at least the groups of nodes corresponding to each region of the HYDRO case"
     return {}
-    
-  
+
+
   for grp in groups:
     if grp.GetType() == SMESH.NODE:
       grpName = grp.GetName()
       print grpName
       if grpName in dicoGroupeRegion.keys():
         regions[dicoGroupeRegion[grpName]] = grp
-      
+
   fo = open(fichierFonds, 'w')
   statz = {}
   for nomreg, grp in regions.iteritems():
@@ -223,11 +227,16 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod
       vx.append(xyz[0])
       vy.append(xyz[1])
     vz = cas.GetAltitudesForPoints( vx, vy, region, interpolMethod )
-    minz = min(vz)
-    maxz = max(vz)
-    statz[grp.GetName()] = (minz, maxz)
-
+    minz = np.amin(vz)
+    maxz = np.amax(vz)
+    meanz = np.mean(vz)
+    stdz = np.std(vz)
+    v05z = np.percentile(vz,05)
+    v95z = np.percentile(vz,95)
     
+    statz[grp.GetName()] = (minz, maxz, meanz, stdz, v05z, v95z)
+
+
     for i,nodeId in enumerate(nodesIds):
       x = vx[i]
       y = vy[i]