X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FHYDROTools%2FinterpolZ.py;h=9f1c978fb5cb67d4b68ebfda37667fbc329a14ab;hb=ef44f612e7baa53c452b3ceb2b7202b8330eb234;hp=268c0e350cf57ad4f65061f7b565d51340f4d2f0;hpb=43f06c83878240e7c9bdc9aa33b84115c6a56022;p=modules%2Fhydro.git diff --git a/src/HYDROTools/interpolZ.py b/src/HYDROTools/interpolZ.py index 268c0e35..9f1c978f 100644 --- a/src/HYDROTools/interpolZ.py +++ b/src/HYDROTools/interpolZ.py @@ -1,7 +1,6 @@ # -*- coding: utf-8 -*- - """ -example of use case of interpolZ and createZfield methods: +Example of use case of interpolZ with the default values: # --- case name in HYDRO nomCas = 'inondation1' @@ -13,234 +12,340 @@ fichierMaillage = '/home/B27118/projets/LNHE/garonne/inondation1.med' 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) -zUndef = 90 # --- Z interpolation on the bathymety/altimetry on the mesh nodes -interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef) - -# --- add a field on nodes of type double with z values, named "BOTTOM" -createZfield2(fichierMaillage) +interpolZ(nomCas, fichierMaillage, dicoGroupeRegion) """ +__revision__ = "V3.01" # ----------------------------------------------------------------------------- -import string - # ----------------------------------------------------------------------------- -import sys import salome salome.salome_init() theStudy = salome.myStudy theStudyId = salome.myStudyId +import numpy as np +import MEDLoader as ml +import MEDCoupling as mc + # ----------------------------------------------------------------------------- from med import medfile from med import medmesh -#from med import medfilter from med import medfield from med import medenum -#from med import medprofile -#from med import medlocalization -#from med import medlink - -def createZfield1(fichierMaillage): - """ - Complete the mesh for Telemac. - Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes. - createZfield1 is used after interpolZ. createZfield1 is base on med file interface. - There is an alternate method based on MEDLoader, equivalent (createZfield2). - The file F.med produced by interpolz must exist, and is modified. - fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ. - return F.med : med file containing the field "BOTTOM" - """ - - noms = string.split(fichierMaillage,'.') - basename = string.join(noms[:-1], '.') - fichierFMaillage = basename + 'F.med' - print fichierFMaillage - - # --- ouverture du fichier - fid=medfile.MEDfileOpen(fichierFMaillage, medenum.MED_ACC_RDEXT) - - maa, sdim, mdim, meshtype, desc, dtunit, sort, nstep, repere, axisname, axisunit = medmesh.MEDmeshInfo(fid, 1) - print "Maillage de nom : %s , de dimension : %d , et de type %s"%(maa,mdim,meshtype) - print " Dimension de l'espace : %d"%(sdim) - - # --- Combien de noeuds a lire ? - nnoe, chgt, trsf = medmesh.MEDmeshnEntity(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT, - medenum.MED_NODE, medenum.MED_NONE, - medenum.MED_COORDINATE, medenum.MED_NO_CMODE) - - if nnoe > 0: - # --- Allocations memoire - # --- table des coordonnees flt : (dimension * nombre de noeuds ) - coords = medfile.MEDFLOAT(nnoe*sdim) - # --- table des numeros des noeuds - numnoe = medfile.MEDINT(nnoe) - - # --- 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 - nomcha1 = "BOTTOM" - ncomp1 = 1 - # --1234567890123456-- - comp1 = "z " - unit1 = "m " - dtunit = "" - medfield.MEDfieldCr(fid, nomcha1, medfile.MED_FLOAT64, - ncomp1, comp1, unit1, dtunit, maa) - - # --- ecriture du champ - - medfield.MEDfieldValueWr(fid, nomcha1, medenum.MED_NO_DT, medenum.MED_NO_IT, 0.0, - medenum.MED_NODE, medenum.MED_NONE, - medenum.MED_FULL_INTERLACE, medenum.MED_ALL_CONSTITUENT, nnoe, valz) - # --- fermeture du fichier - medfile.MEDfileClose(fid) - -# ----------------------------------------------------------------------------- - -from MEDLoader import MEDLoader, MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh - -def createZfield2(fichierMaillage): - """ - Complete the mesh for Telemac. - Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes. - createZfield2 is used after interpolZ. createZfield1 is base on MEDLoader interface. - There is an alternate method based on Med file, equivalent (createZfield1). - fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ. - return F.med : med file containing the field "BOTTOM" - """ - - noms = string.split(fichierMaillage,'.') - basename = string.join(noms[:-1], '.') - fichierZMaillage = basename + 'Z.med' - fichierFMaillage = basename + 'F.med' - print fichierFMaillage - - mymesh = MEDLoader.ReadUMeshFromFile(fichierZMaillage,0) - fieldOnNodes=MEDCouplingFieldDouble.New(ON_NODES) - fieldOnNodes.setName("BOTTOM") - fieldOnNodes.setMesh(mymesh) - fieldOnNodes.setArray(mymesh.getCoords()[:,2]) - mm=MEDFileMesh.New(fichierZMaillage) - mm.write(fichierFMaillage,2) - MEDLoader.WriteFieldUsingAlreadyWrittenMesh(fichierFMaillage,fieldOnNodes) # ----------------------------------------------------------------------------- import HYDROPy class MyInterpolator( HYDROPy.HYDROData_IInterpolator ): + """ + Class MyInterpolator + """ + def __init__ (self) : + """ +Constructor + """ def GetAltitudeForPoint( self, theCoordX, theCoordY ): - alt_obj = HYDROPy.HYDROData_IInterpolator.GetAltitudeObject( self ); + """ + Function + """ + alt_obj = HYDROPy.HYDROData_IInterpolator.GetAltitudeObject( self ) z = alt_obj.GetAltitudeForPoint( theCoordX, theCoordY ) # Custom calculation takes the base value and changes it to test #z2 = (z - 74.0)*10 z2 = z return z2 - -# ----------------------------------------------------------------------------- -import SMESH -import SALOMEDS -from salome.smesh import smeshBuilder +# ----------------------------------------------------------------------------- -smesh = smeshBuilder.New(theStudy) -def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef): +def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMethod=0, m3d=False, xyzFile=False, verbose=False): """ interpolZ takes a 2D (x,y) mesh and calls the active instance of module HYDRO to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node. - interpolZ must be followed by createZfield1 or createZfield2. - nomCas: Calculation Case Name in module HYDRO - fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case - dicoGroupeRegion: python dictionary giving the coorespondance of mesh groups to HYDRO regions. - Key: face group name, value: region name in the HYDRO Case - zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct). - return Z.med : med file with Z value on nodes - return F.med : an exact copy of Z.med - """ - doc = HYDROPy.HYDROData_Document.Document( theStudyId ) - cas = doc.FindObjectByName(nomCas) - print cas - custom_inter = MyInterpolator() - - noms = string.split(fichierMaillage,'.') - basename = string.join(noms[:-1], '.') - fichierZMaillage = basename + 'Z.med' - fichierFMaillage = basename + 'F.med' - fichierFonds = basename + '.xyz' - - print fichierMaillage - 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" - 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(): - print "------- region : ", nomreg - region = doc.FindObjectByName(nomreg) - #print region - #region.SetInterpolator(custom_inter) - nodesIds = grp.GetListOfID() - #print nodesIds - vx = [] - vy = [] - for nodeId in nodesIds: - xyz = maillagePlat.GetNodeXYZ(nodeId) - #print xyz - vx.append(xyz[0]) - vy.append(xyz[1]) - vz = cas.GetAltitudesForPoints( vx, vy, region ) - minz = min(vz) - maxz = max(vz) - statz[grp.GetName()] = (minz, maxz) - - - for i,nodeId in enumerate(nodesIds): - x = vx[i] - y = vy[i] - z = vz[i] - if z < -9000: - z = zUndef - #print i, nodeId, x, y, z - maillagePlat.MoveNode(nodeId, x, y, z) - l = "%10.2f %10.2f %10.2f\n"%(x, y, z) - fo.write(l) - - fo.close() - maillagePlat.ExportMED(fichierZMaillage , 0, SMESH.MED_V2_2, 1 ) - maillagePlat.ExportMED(fichierFMaillage , 0, SMESH.MED_V2_2, 1 ) + In: + nomCas: Calculation Case Name in module HYDRO + fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case + dicoGroupeRegion: python dictionary giving the correspondance of mesh groups to HYDRO regions. + Key: face group name + Value: region name in the HYDRO Case + zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct). + default value is 90. + interpolMethod: integer value + 0 = nearest point on bathymetry (default) + 1 = linear interpolation + m3d: True/False to produce a 3D mesh. Default is False. + xyzFile: True/False to write an ascii file with xyz for every node. Default is False. + Out: + statz: statistique for z + Key: face group name + Value: (minz, maxz, meanz, stdz, v05z, v95z) + Out: + return F.med : med file with Z value in a field "BOTTOM" + Option: Z value also in Z coordinate if m3D is true + return .xyz : text file with X, Y, Z values + """ + statz = dict() + erreur = 0 + message = "" + + while not erreur: + + if verbose: + ligne = "nomCas: %s" % nomCas + ligne += "\ninterpolMethod: %d" % interpolMethod + if (zUndef != None ): + ligne += "\nzUndef: %f" % zUndef + ligne += "\nm3d: %d" % m3d + print (ligne) + + doc = HYDROPy.HYDROData_Document.Document(theStudyId) + cas = doc.FindObjectByName(nomCas) + print ( "cas : ", cas) + custom_inter = MyInterpolator() + + basename = fichierMaillage[:-4] + fichierFMaillage = basename + 'F.med' + + print ("dicoGroupeRegion = ", dicoGroupeRegion) + 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() + dicoGroupeRegion_0 = dict() + nb_pb = 0 + for gr_face_name in dicoGroupeRegion: + saux = gr_face_name.strip() + if saux not in t_group_n: + message += "Group: '" + gr_face_name + "'\n" + nb_pb += 1 + else : + dicoGroupeRegion_0[saux] = dicoGroupeRegion[gr_face_name] + if verbose: + ligne = "Number of problems: %d" % nb_pb + print (ligne) +# + if nb_pb > 0: + if nb_pb == 1: + message += "This group does" + else: + message += "These %d groups do" % nb_pb + message += " not belong to the mesh.\n" + message += "Please check the names of the group(s) of faces corresponding to each region of the HYDRO case.\n" + message += "Groups for this file:\n" + for group_n in t_group_n : + message += "'%s'\n" % group_n + erreur = 2 + break +# +# 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) +# + for gr_face_name in dicoGroupeRegion_0: +# +# 4.1. Region connected to the group +# + nomreg = dicoGroupeRegion_0[gr_face_name] + ligne = "------- Region: '%s'" % nomreg + ligne += ", connected to group '%s'" % gr_face_name + print (ligne) + region = doc.FindObjectByName(nomreg) +# +# 4.2. Mesh of the group +# + mesh_of_the_group = meshMEDFileRead.getGroup(0, gr_face_name, False) + nbr_cells = mesh_of_the_group.getNumberOfCells() + if verbose: + ligne = "\t. Number of cells: %d" % nbr_cells + print (ligne) +# +# 4.3. 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.4. 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 = cas.GetAltitudesForPoints(vx, vy, region, interpolMethod) +# + #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.5. 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) +# + break +# + if erreur: + print message +# return statz -# ----------------------------------------------------------------------------- -