1 # -*- coding: utf-8 -*-
3 # -----------------------------------------------------------------------------
7 from scipy import spatial
11 def changeBathy(meshFileIn, meshFileOut, bathyFile, groupName, offsetX=0., offsetY=0., zUndef=-100., interpMethod=0, m3d=False, stats=False):
13 ChangeBathy takes a 2D (x,y) mesh with an already existing field of double ('BOTTOM') for altimetry.
14 The mesh may already have a z coordinate set with the z altimetry.
15 The new altimetry from the bathymetry file is set on the nodes of the selected group of face.
18 meshFileIn: med file name produced by SMESH, corresponding to the HYDRO case
19 bathyFile: bathymetry file
20 groupName: name of the group of faces to modify.
21 offsetX, offsetY: coordinates of origin in MESH
22 zUndef: Z value to use for nodes outside the bathymetry interpolation (should not happen: if so, algoritm problem).
23 default value is -100.
24 interpMethod: integer value
25 0 = nearest point on bathymetry
26 1 = linear interpolation
27 m3d: True/False to produce a 3D mesh. Default is False.
28 stats: compute statistics for SALOME test
30 statz: statistique for z
32 Value: (minz, maxz, meanz, stdz, v05z, v95z)
34 meshFileOut : med file with Z value in a field "BOTTOM"
35 Option: Z value also in Z coordinate if m3D is true
38 # --- load the bathymetry in a KD tree, working with the same origin coordinates as the mesh
40 fi=open(bathyFile, 'r')
55 x = float(vals[0]) - offsetX
56 y = float(vals[1]) - offsetY
68 print("bornes x", minx, maxx)
69 print("bornes y", miny, maxy)
70 print("bornes z", minz, maxz)
71 print("nombre points:", len(alti))
73 bathyPts = np.array(coords)
76 tree = spatial.KDTree(bathyPts)
78 # --- load the nodes coordinates of the mesh group
80 mcMesh = ml.MEDFileMesh.New(meshFileIn)
81 dim = mcMesh.getMeshDimension()
82 d0 = 0 # when dimension 2, faces are dim
83 if dim == 3: # when dimension 3, edges are dim -1
86 groupNames = mcMesh.getGroupsOnSpecifiedLev(d0) #names of face groups
88 if groupName not in groupNames:
89 print("there is no group of faces named %s in the mesh" %groupName)
92 group = mcMesh.getGroup(d0, groupName)
94 # nodes ids in use: number of nodes and order of the mesh, nodes not in the group marked -1,
95 # nodes indices = 0 to nbNodeInGroup, not their indice in mesh!
97 AllNodeIds = group.getNodeIdsInUse()[0].toNumPyArray()
98 print("AllNodeIds.size:", AllNodeIds.size)
99 cond = AllNodeIds != -1
100 #nodeIds = np.extract(cond, AllNodeIds)
101 nids = [i for i in range(AllNodeIds.size) if int(AllNodeIds[i]) >=0] # indices in mesh
102 nodeIds = np.array(nids)
103 print("nodeIds.size:", nodeIds.size)
104 print("nodeIds:", nodeIds)
106 allCoords = group.getCoords().toNumPyArray()
107 print("allCoords.shape:",allCoords.shape)
108 dimCoord = allCoords.shape[1]
110 condition = np.vstack((cond, cond, cond)).transpose()
112 condition = np.vstack((cond, cond)).transpose()
113 coords = np.extract(condition, allCoords)
114 nbRaws = coords.size // dimCoord
115 coords = coords.reshape(nbRaws, dimCoord)
116 print("coords.shape:",coords.shape)
118 if coords.shape[1] == 3:
119 coordsXY = coords[:, 0:2] # 2 first columns
120 coordsZ = coords[:, 2] # 3rd column
125 if interpMethod != 0:
126 print("TODO: the linear interpolation on bathymetry is not yet implemented, using nearest point")
128 # --- get nearest node Ids for nodes of group
130 (distances, ids) = tree.query(coordsXY)
131 print("ids.size", ids.size)
134 # --- check if there is a 'BOTTOM' Field
137 bottomField = ml.ReadField(meshFileIn, 'BOTTOM')
138 bottomArray = bottomField.getArray()
139 except ml.InterpKernelException:
140 print("no BOTTOM array")
143 if bottomField is None:
144 bottomField = ml.MEDCouplingFieldDouble(ml.ON_NODES)
145 bottomField.setName('BOTTOM')
146 bottomField.setMesh(mcMesh.getMeshAtLevel(0))
147 npField = np.zeros((mcMesh.getNumberOfNodes()))
148 bottomArray = ml.DataArrayDouble(npField)
150 # --- set altitude on field
153 for i,id in enumerate(ids):
154 #print(i, "---", nodeIds[i], "---", id, alti[id])
156 bottomArray[j] = alti[id]
160 bottomField.setArray(bottomArray)
161 bottomField.setTime(0.0, 0, -1)
163 if (coordsZ is not None) or m3d:
164 print("TODO: coordsZ")
165 mlCoords = group.getCoords()
166 coords3D = ml.DataArrayDouble.Meld([mlCoords[:,0:2], bottomArray])
167 coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
168 mcMesh.setCoords(coords3D)
170 mcMesh.write(meshFileOut, 2)
171 ml.WriteFieldUsingAlreadyWrittenMesh(meshFileOut, bottomField)
179 v05z = np.percentile(vz, 0o5)
180 v95z = np.percentile(vz, 95)
181 ligne = ".. Minimum: %f" % minz
182 ligne += ", maximum: %f" % maxz
183 ligne += ", mean: %f\n" % meanz
184 ligne += ".. stdeviation: %f" % stdz
185 ligne += ", v05z: %f" % v05z
186 ligne += ", v95z: %f" % v95z
188 statz[groupName] = (minz, maxz, meanz, stdz, v05z, v95z)
192 if __name__=='__main__':
193 meshFileIn = 'garonne_1F.med'
194 meshFileOut = 'garonne_1FF.med'
195 bathyFile = 'newBathy.xyz'
196 #groupName = 'litMineur'
197 groupName = 'riveGauche'
198 #groupName = 'riveDroite'
201 stats = changeBathy(meshFileIn, meshFileOut, bathyFile, groupName, offsetX, offsetY, stats=True)