From: isn <isn@opencascade.com>
Date: Mon, 28 Jan 2019 13:02:16 +0000 (+0300)
Subject: lot 19 : interpolZ_B method for interpolZ on bathymetry
X-Git-Tag: SH_V2_2_0~99^2~37
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=51af016ce2c42266bfcd51bac66ed86f1ce284f9;p=modules%2Fhydro.git

lot 19 : interpolZ_B method for interpolZ on bathymetry
---

diff --git a/src/HYDROTools/interpolZ.py b/src/HYDROTools/interpolZ.py
index 1d12cfe0..d1454d5f 100644
--- a/src/HYDROTools/interpolZ.py
+++ b/src/HYDROTools/interpolZ.py
@@ -360,3 +360,216 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., regions_int
 #
   return statz
 
+
+def interpolZ_B(bathyName, fichierMaillage, gr_face_name, zUndef=90., interp_method=0, m3d=False, xyzFile=False, verbose=False):
+
+  doc = HYDROPy.HYDROData_Document.Document(theStudyId)
+  bathy_obj = doc.FindObjectByName(bathyName)
+  print ( "bathy : ", bathy_obj)
+  if bathy_obj is None:
+    print ( "bathy is None")
+    return False
+
+  basename = fichierMaillage[:-4]
+  fichierFMaillage = basename + 'F.med'
+
+  ligne = "fichierMaillage  = %s" % fichierMaillage
+  ligne += "\nfichierFMaillage = %s" % fichierFMaillage
+  
+  if xyzFile:
+    fichierFonds = basename + '.xyz'
+    ligne += "\nfichierFonds     = %s" % fichierFonds
+  print (ligne)
+#
+# 1. Reads the mesh
+#
+  meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
+  if verbose:
+    print (meshMEDFileRead)
+#
+# 2. Checks the names of the groups of faces
+#
+  t_group_n = meshMEDFileRead.getGroupsNames()
+  gr_face_name_tr = gr_face_name.strip()
+  if gr_face_name_tr not in t_group_n:
+    print "Group not found"
+    return False	  
+#
+# 3. Gets the information about the nodes
+#
+  nbnodes = meshMEDFileRead.getNumberOfNodes()
+  if verbose:
+    ligne = "Number of nodes: %d" % nbnodes
+    print (ligne)
+#
+  coords = meshMEDFileRead.getCoords()
+  #print coords
+  if verbose:
+    nb_comp = coords.getNumberOfComponents()
+    l_info = coords.getInfoOnComponents()
+    ligne = ""
+    l_info_0=["X", "Y", "Z"]
+    for id_node in (0, 1, nbnodes-1):
+      ligne += "\nNode #%6d:" % id_node
+      for iaux in range(nb_comp):
+        if l_info[iaux]:
+          saux = l_info[iaux]
+        else:
+          saux = l_info_0[iaux]
+        ligne += " %s" % saux
+        ligne += "=%f" % coords[id_node, iaux]
+    print (ligne)
+#
+# 4. Exploration of every group of faces
+#
+  tb_aux = np.zeros(nbnodes, dtype=np.bool)
+#
+  bathy = np.zeros(nbnodes, dtype=np.float)
+  bathy.fill(zUndef)
+
+#
+# 4.1. Mesh of the group
+#
+  mesh_of_the_group = meshMEDFileRead.getGroup(0, gr_face_name_tr, False)
+  nbr_cells = mesh_of_the_group.getNumberOfCells()
+  if verbose:
+    ligne = "\t. Number of cells: %d" % nbr_cells
+    print (ligne)
+#
+# 4.2. Nodes of the meshes of the group
+#      Every node is flagged in tb_aux
+#
+  tb_aux.fill(False)
+  for id_elem in range(nbr_cells):
+    l_nodes = mesh_of_the_group.getNodeIdsOfCell(id_elem)
+    #print l_nodes
+    for id_node in l_nodes:
+      tb_aux[id_node] = True
+  np_aux = tb_aux.nonzero()
+  if len(np_aux[0]):
+    if verbose:
+      ligne = "\t. Number of nodes for this group: %d" % len(np_aux[0])
+      print (ligne)
+  #print ("np_aux:", np_aux)
+#
+# 4.3. Interpolation over the nodes of the meshes of the group
+#
+  vx = list()
+  vy = list()
+  for nodeId in np_aux[0]:
+    vx.append(coords[nodeId, 0])
+    vy.append(coords[nodeId, 1])
+  #print ("vx:\n", vx)
+  #print ("vy:\n", vy)
+#
+  vz = bathy_obj.GetAltitudesForPoints(vx, vy, interp_method)
+#
+  #print ("vz:\n", vz)
+  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)
+#
+  if verbose:
+    ligne = ".. Minimum: %f" % minz
+    ligne += ", maximum: %f" % maxz
+    ligne += ", mean: %f\n" % meanz
+    ligne += ".. stdeviation: %f" % stdz
+    ligne += ", v05z: %f" % v05z
+    ligne += ", v95z: %f" % v95z
+    print (ligne)
+#
+# 4.4. Storage of the z and of the statistics for this region
+#
+  statz_gr_face_name = (minz, maxz, meanz, stdz, v05z, v95z)
+#
+  for iaux, nodeId in enumerate(np_aux[0]):
+    bathy[nodeId] = vz[iaux]
+#
+# 5. Minimum:
+#    During the interpolation, if no value is available over a node, a default value
+#    is set: -9999. It has no importance for the final computation, but if the field
+#    or the mesh is displayed, it makes huge gap. To prevent this artefact, a more
+#    convenient "undefined" value is set. This new undefined value is given by the user.
+#
+#    zUndefThreshold: the default value is -9000. It is tied with the value -9999. given
+#                     by the interpolation when no value is defined.
+#
+  zUndefThreshold = -9000.
+  if verbose:
+    ligne = "zUndefThreshold: %f" % zUndefThreshold
+    print (ligne)
+#
+  #print "bathy :\n", bathy
+  np_aux_z = (bathy < zUndefThreshold).nonzero()
+  if verbose:
+    ligne = ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
+    print (ligne)
+  if len(np_aux_z[0]):
+    for iaux in np_aux_z[0]:
+      bathy[iaux] = zUndef
+#
+# 6. Option : xyz file
+#
+  if xyzFile:
+    if verbose:
+      ligne = ".. Ecriture du champ de bathymétrie sur le fichier :\n%s" % fichierFonds
+      print (ligne)
+#
+    with open(fichierFonds, "w") as fo :
+      for nodeId in range(nbnodes):
+        ligne = "%11.3f %11.3f %11.3f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
+        fo.write(ligne)
+#
+# 7. Final MED file
+# 7.1. Transformation of the bathymetry as a double array
+#
+  bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
+  bathy_dd.setInfoOnComponents(["Z [m]"])
+#
+# 7.2. If requested, modification of the z coordinate
+#
+  if m3d:
+    coords3D = ml.DataArrayDouble.Meld([coords[:,0:2], bathy_dd])
+    coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
+    #print "coords3D =\n", coords3D
+    meshMEDFileRead.setCoords(coords3D)
+#
+# 7.3. Writes the mesh
+#
+  if verbose:
+    if m3d:
+      saux = " 3D"
+    else:
+      saux = ""
+    ligne = ".. Ecriture du maillage" + saux
+    ligne += " sur le fichier :\n%s" % fichierFMaillage
+    print (ligne)
+#
+  meshMEDFileRead.write(fichierFMaillage, 2)
+#
+# 7.4. Writes the field
+#
+  med_field_name = "BOTTOM"
+  if verbose:
+    ligne = ".. Ecriture du champ '%s'" % med_field_name
+    print (ligne)
+#
+  #print "bathy_dd =\n", bathy_dd
+  fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES)
+  fieldOnNodes.setName(med_field_name)
+  fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0))
+  fieldOnNodes.setArray(bathy_dd)
+# Ces valeurs d'instants sont mises pour assurer la lecture par TELEMAC
+# instant = 0.0
+# numero d'itération : 0
+# pas de numero d'ordre (-1)
+  fieldOnNodes.setTime(0.0, 0, -1)
+#
+  fMEDFile_ch_d = ml.MEDFileField1TS()
+  fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
+  fMEDFile_ch_d.write(fichierFMaillage, 0)
+#
+  return statz_gr_face_name