From 51af016ce2c42266bfcd51bac66ed86f1ce284f9 Mon Sep 17 00:00:00 2001 From: isn Date: Mon, 28 Jan 2019 16:02:16 +0300 Subject: [PATCH] lot 19 : interpolZ_B method for interpolZ on bathymetry --- src/HYDROTools/interpolZ.py | 213 ++++++++++++++++++++++++++++++++++++ 1 file changed, 213 insertions(+) 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 -- 2.39.2