# -*- coding: utf-8 -*-
-
"""
-example of use case of interpolZ:
+Example of use case of interpolZ with the default values:
# --- case name in HYDRO
nomCas = 'inondation1'
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, interpolMethod)
+interpolZ(nomCas, fichierMaillage, dicoGroupeRegion)
"""
-__revision__ = "V2.04"
+__revision__ = "V3.01"
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
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
-
-# La fonction createZfield1 ne sert plus à partir du 12/07/2017
-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 <fichierMaillage>F.med produced by interpolz must exist, and is modified.
- fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
- return <fichierMaillage>F.med : med file containing the field "BOTTOM"
- """
-
- basename = fichierMaillage[:-4]
- 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, 0, 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)
- print fichierFMaillage, " field BOTTOM OK"
-
-# -----------------------------------------------------------------------------
-
-# La fonction createZfield2 ne sert plus à partir du 12/07/2017
-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 <fichierMaillage>F.med : med file containing the field "BOTTOM"
- """
-
- import MEDLoader
- from MEDLoader import MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh
-
- basename = fichierMaillage[:-4]
- 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])
- fieldOnNodes.setTime(0.0,0,-1)
- mm=MEDFileMesh.New(fichierZMaillage)
- mm.write(fichierFMaillage,2)
- MEDLoader.WriteFieldUsingAlreadyWrittenMesh(fichierFMaillage,fieldOnNodes)
- print fichierFMaillage, " field BOTTOM OK"
# -----------------------------------------------------------------------------
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
# -----------------------------------------------------------------------------
-def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMethod=0, xyzFile=False, verbose=False):
+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.
interpolMethod: integer value
0 = nearest point on bathymetry (default)
1 = linear interpolation
- zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
+ 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 <fichierMaillage>F.med : med file with Z value on nodes and in a field "BOTTOM"
+ return <fichierMaillage>F.med : med file with Z value in a field "BOTTOM"
+ Option: Z value also in Z coordinate if m3D is true
+ return <fichierMaillage>.xyz : text file with X, Y, Z values
"""
statz = dict()
erreur = 0
while not erreur:
if verbose:
- print "nomCas:", nomCas
- print "interpolMethod: %d" % interpolMethod
- print "zUndef:", zUndef
+ 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
+ print ( "cas : ", cas)
custom_inter = MyInterpolator()
basename = fichierMaillage[:-4]
fichierFMaillage = basename + 'F.med'
- print "dicoGroupeRegion =", dicoGroupeRegion
- print "fichierMaillage =", fichierMaillage
- print "fichierFMaillage =", fichierFMaillage
+ print ("dicoGroupeRegion = ", dicoGroupeRegion)
+ ligne = "fichierMaillage = %s" % fichierMaillage
+ ligne += "\nfichierFMaillage = %s" % fichierFMaillage
if xyzFile:
fichierFonds = basename + '.xyz'
- print "fichierFonds =", fichierFonds
+ 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
#
- l_gr_faces = list(dicoGroupeRegion.keys())
- t_groupe_n = meshMEDFileRead.getGroupsNames()
+ t_group_n = meshMEDFileRead.getGroupsNames()
+ dicoGroupeRegion_0 = dict()
nb_pb = 0
- for gr_face_name in l_gr_faces:
- if gr_face_name not in t_groupe_n:
+ 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:
- print "Number of problems: %d" % nb_pb
+ ligne = "Number of problems: %d" % nb_pb
+ print (ligne)
#
if nb_pb > 0:
if nb_pb == 1:
message += "This group does"
else:
- message += "That %d groups do" % nb_pb
- message += " not belongs to the mesh.\n"
- message += "Please check the names of the group(s) of faces corresponding to each region of the HYDRO case"
+ 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
#
#
nbnodes = meshMEDFileRead.getNumberOfNodes()
if verbose:
- print "Number of nodes: %d" % nbnodes
+ ligne = "Number of nodes: %d" % nbnodes
+ print (ligne)
#
coords = meshMEDFileRead.getCoords()
#print coords
- #print coords[0,0]
- #print coords[0,1]
+ 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
#
bathy = np.zeros(nbnodes, dtype=np.float)
bathy.fill(zUndef)
#
- for gr_face_name in l_gr_faces:
+ for gr_face_name in dicoGroupeRegion_0:
#
# 4.1. Region connected to the group
#
- nomreg = dicoGroupeRegion[gr_face_name]
- line = "------- Region: '" + nomreg + "'"
- line += ", connected to group '" + gr_face_name + "'"
- print line
+ nomreg = dicoGroupeRegion_0[gr_face_name]
+ ligne = "------- Region: '%s'" % nomreg
+ ligne += ", connected to group '%s'" % gr_face_name
+ print (ligne)
region = doc.FindObjectByName(nomreg)
- #print region
- #region.SetInterpolator(custom_inter)
#
# 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:
- print "\t. Number of cells: %d" % nbr_cells
+ 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
np_aux = tb_aux.nonzero()
if len(np_aux[0]):
if verbose:
- print "\t. Number of nodes for this group: %d" % len(np_aux[0])
- #print "np_aux:", np_aux
+ 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 = []
- vy = []
+ 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
+ #print ("vx:\n", vx)
+ #print ("vy:\n", vy)
+#
vz = cas.GetAltitudesForPoints(vx, vy, region, interpolMethod)
- #print "vz:\n", vz
+#
+ #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 += ".. stdeviation: %f" % stdz
ligne += ", v05z: %f" % v05z
ligne += ", v95z: %f" % v95z
- print ligne
+ print (ligne)
#
# 4.5. Storage of the z and of the statistics for this region
#
# 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 undefinde value is given by the user.
+# 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:
- print "zUndefThreshold: %f" % zUndefThreshold#
+ ligne = "zUndefThreshold: %f" % zUndefThreshold
+ print (ligne)
#
#print "bathy :\n", bathy
np_aux_z = (bathy < zUndefThreshold).nonzero()
if verbose:
- print ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
+ 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. xyz file
+# 6. Option : xyz file
#
if xyzFile:
#
if verbose:
- print ".. Ecriture du champ de bathymétrie sur le fichier :\n", fichierFonds
- fo = open(fichierFonds, 'w')
- for nodeId in range(nbnodes):
- line = "%10.2f %10.2f %10.2f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
- fo.write(line)
- fo.close()
+ 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. Modification of the z coordinates
+# 7.1. Transformation of the bathymetry as a double array
#
bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
bathy_dd.setInfoOnComponents(["Z [m]"])
#
- coords3D = ml.DataArrayDouble.Meld([coords, bathy_dd])
- coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
- #print "coords3D =\n", coords3D
+# 7.2. If requested, modification of the z coordinate
#
- meshMEDFileRead.setCoords(coords3D)
+ 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.2. Writes the 3D mesh
+# 7.3. Writes the mesh
#
if verbose:
- print ".. Ecriture du maillage 3D sur le fichier :\n", fichierFMaillage
+ 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.3. Writes the field
+# 7.4. Writes the field
#
med_field_name = "BOTTOM"
if verbose:
- print ".. Ecriture du champ '"+med_field_name+"'"
+ 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)
fMEDFile_ch_d.write(fichierFMaillage, 0)
#
break
-
+#
if erreur:
print message
-
+#
return statz
