-
-# -----------------------------------------------------------------------------
-
import salome
salome.salome_init()
import medcoupling as mc
import shapefile
+import math
+import os
+import json
-def freeBordersGroup(meshFile):
- print(" === freeBordersGroup", meshFile)
+def freeBordersGroup(meshFileIn, meshFileOut=""):
+ """
+ In a mesh, create a group of Edges for the borders: domain, internal holes (isles)
+ parameters:
+ meshFileIn: full path of the input mesh file, format MED
+ meshFileOut: full path of the output mesh file, format MED (default="" : when "", output file is suffixed with "_brd.med"
+ return full path of the output mesh file
+ """
smesh = smeshBuilder.New()
- smesh.SetEnablePublish( False ) # Set to False to avoid publish in study if not needed
- ([MESH], status) = smesh.CreateMeshesFromMED(meshFile)
- groups = MESH.GetGroups()
+ if not salome.sg.hasDesktop():
+ smesh.SetEnablePublish( False ) # Set to False to avoid publish in study if not needed
+ ([MESH], status) = smesh.CreateMeshesFromMED(meshFileIn)
+
+ nbAdded, MESH, addedBnd = MESH.MakeBoundaryElements( SMESH.BND_1DFROM2D, '', '', 0, [])
+
aCriteria = []
aCriterion = smesh.GetCriterion(SMESH.EDGE,SMESH.FT_FreeBorders,SMESH.FT_Undefined,0)
aCriteria.append(aCriterion)
aFilter = smesh.GetFilterFromCriteria(aCriteria)
aFilter.SetMesh(MESH.GetMesh())
FreeBorders = MESH.GroupOnFilter( SMESH.EDGE, 'FreeBorders', aFilter )
- smesh.SetName(MESH, 'MESH')
- newMeshName = '/tmp/freeBorders.med'
+
+ a = os.path.splitext(meshFileIn)
+ smesh.SetName(MESH, os.path.basename(a[0]))
+
+ if meshFileOut == "":
+ a = os.path.splitext(meshFileIn)
+ smesh.SetName(MESH, os.path.basename(a[0]))
+ meshFileOut = a[0] + '_brd' + a[1]
+
try:
- MESH.ExportMED(newMeshName,auto_groups=0,minor=41,overwrite=1,meshPart=None,autoDimension=1)
+ MESH.ExportMED(meshFileOut,auto_groups=0,minor=40,overwrite=1,meshPart=None,autoDimension=1)
pass
except:
print('ExportMED() failed. Invalid file name?')
- return newMeshName
+ return meshFileOut
+
def explodeGroup(grp, grpName):
+ """
+ from a group of edges loaded with MEDCoupling, create ordered lists of nodes, one list for each set of connected edges.
+ parameters:
+ grp: MEDCoupling object for the group of edges
+ grpName: name of the group
+ return:
+ List of descriptors [(ordered list of nodeIds, name of the list, closed status)]
+ """
print(" === explodeGroup", grpName)
nbCells=grp.getNumberOfCells()
- dicReverse = {} # id noeud --> id edges
+ dicReverse = {} # id node --> id edges
for i in range(nbCells):
nodcell = grp.getNodeIdsOfCell(i)
for j in range(len(nodcell)):
print(chainDesc[1:])
return nodeChains
-def writeShapeLines(mcMesh, grpName, nodeChains, offsetX=0., offsetY=0.):
+
+def writeShapeLines(mcMesh, grpName, nodeChains, outputDirectory, offsetX=0., offsetY=0.):
+ """
+ from a mesh loaded in memory with MEDLoader, and a list of list of connected nodes associated to a group of edges, write a shapefile of type line, with one record per list of connected nodes.
+ parameters:
+ mcMesh: mesh loaded in memory with MEDLoader
+ grpName: name associated to the group of edges
+ nodeChains: List of descriptors corresponding to the group [(ordered list of nodeIds, name of the list, closed status)]
+ outputDirectory: directory for writing the shapefile
+ offsetX : offset of origin X in the mesh (default 0). The shapefile is always written without local origin to be ready for a direct load in Qgis.
+ offsetY : offset of origin Y in the mesh (default 0). The shapefile is always written without local origin to be ready for a direct load in Qgis.
+ """
print(" === writeShapeLines", grpName)
coords = mcMesh.getCoords()
- w = shapefile.Writer(grpName)
+ shapeFileName = os.path.join(outputDirectory, grpName)
+ w = shapefile.Writer(shapeFileName)
w.shapeType = 3
w.field('name', 'C')
for (chain, chainName, closed) in nodeChains:
for node in chain:
coord = coords[node].getValues()
coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
- #print(" ", coordLb93)
chaincoords.append(coordLb93)
w.line([chaincoords])
w.record(chainName)
w.close()
+ return shapeFileName + '.shp'
-def writeShapePoints(mcMesh, grpName, nodeChains, offsetX=0., offsetY=0.):
+
+def writeShapePoints(mcMesh, grpName, nodeChains, outputDirectory, offsetX=0., offsetY=0.):
+ """
+ from a mesh loaded in memory with MEDLoader, and a list of list of connected nodes associated to a group of edges, write a shapefile of type multi points, with one record per list of connected nodes.
+ parameters:
+ mcMesh: mesh loaded in memory with MEDLoader
+ grpName: name associated to the group of edges
+ nodeChains: List of descriptors corresponding to the group [(ordered list of nodeIds, name of the list, closed status)]
+ outputDirectory: directory for writing the shapefile
+ offsetX : offset of origin X in the mesh (default 0). The shapefile is always written without local origin to be ready for a direct load in Qgis.
+ offsetY : offset of origin Y in the mesh (default 0). The shapefile is always written without local origin to be ready for a direct load in Qgis.
+ """
print(" === writeShapePoints", grpName)
coords = mcMesh.getCoords()
- w = shapefile.Writer(grpName + '_pts')
+ shapeFileName = os.path.join(outputDirectory, grpName)
+ w = shapefile.Writer(shapeFileName + '_pts')
#w.shapeType = 8
w.field('name', 'C')
for (chain, chainName, closed) in nodeChains:
for node in chain:
coord = coords[node].getValues()
coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
- #print(" ", coordLb93)
chaincoords.append(coordLb93)
w.multipoint(chaincoords)
w.record(chainName)
w.close()
-def exploreEdgeGroups(meshFile, offsetX=0., offsetY=0.):
+
+def exploreEdgeGroups(meshFile, outputDirectory="", offsetX=0., offsetY=0.):
+ """
+ Find all the groups of edges in a mesh and, for each group, create one shapefile of lines and one of points. The shapefiles are created in the same system of coordinates as the mesh (For instance Lambert 93), but without origin offset (to be ready for a direct load in Qgis)
+ parameters:
+ meshFile: full path of the input mesh file, format MED
+ outputDirectory: directory in which the shapefiles are written (default "", if "" use the directory containing the mesh
+ offsetX: local X origin of the mesh
+ offsetY: local Y origin of the mesh
+ """
print(" === exploreEdgeGroups", meshFile)
+ if outputDirectory == "":
+ outputDirectory = os.path.dirname(meshFile)
+
+ a = os.path.splitext(meshFile)
+ prefix = os.path.basename(a[0]) # prefix = file name without extension
+
mcMesh = ml.MEDFileMesh.New(meshFile)
- dim = mcMesh.getSpaceDimension()
+ dim = mcMesh.getMeshDimension()
d1=-1 # when dimension 2, edges are dim -1
if dim == 3: # when dimension 3, edges are dim -2
d1=-2
grp_names = mcMesh.getGroupsOnSpecifiedLev(d1) #names of edges groups
- groups = [mcMesh.getGroup(d1, name) for name in grp_names]
+
+ groups = [mcMesh.getGroup(d1, name) for name in grp_names] # list of groups in their name order
+
+ filenames = []
for (grp, grpName) in zip(groups, grp_names):
- nodeChains = explodeGroup(grp, grpName)
- writeShapeLines(mcMesh, grpName, nodeChains, offsetX, offsetY)
- writeShapePoints(mcMesh, grpName, nodeChains, offsetX, offsetY)
+ fullGrpName = prefix + '_' + grpName
+ nodeChains = explodeGroup(grp, fullGrpName)
+ filename = writeShapeLines(mcMesh, fullGrpName, nodeChains, outputDirectory, offsetX, offsetY)
+ writeShapePoints(mcMesh, fullGrpName, nodeChains, outputDirectory, offsetX, offsetY)
+ filenames.append(filename)
+ shapesListFile = os.path.join(outputDirectory, "shapesList.json")
+ with open(shapesListFile, 'w') as f:
+ json.dump(filenames, f)
+
+def fitShapePointsToMesh(freeBorderShapefile, shapefileToAdjust, outputDirectory="", splitFreeBorder=True, splitShapeAdjusted=True):
+ """
+ shapeFileToAdjust must be a closed line or polygon crossing freeBorderShapefile in 2 points.
+ Find in shapeFileToAdjust and in freeBorderShapefile the two closest corresponding points and move the points in shapeFileToAdjust to correspond to the points found in freeBorderShapefile. Split shapeFileToAdjust in two parts (inside or outside freeBorder). If requested, split freeBorderShapefile in two parts. Same for shapeFileToAdjust.
+ parameters:
+ freeBorderShapefile: a set of free border lines, as generated by the functions freeBordersGroup and exploreEdgeGroups.
+ shapefileToAdjust: a closed line or polygon, supposed to be drawn in qgis to pass as close as possible to the points to be connected, on the free border.
+ outputDirectory: if empty, write the resulting shapefiles in their respective directory, with the suffix '_adj', otherwise write in the outputDirectory.
+ splitFreeBorder: boolean default True
+ splitShapeAdjusted: boolean default True
+ """
+ print(" === fitShapePointsToMesh", freeBorderShapefile, shapefileToAdjust)
+
+ # --- find domain freeBorder: bounding box englobing all others
+ # TODO: This may not be always the case, when there is not a single domain with holes, but several non connected domains.
+
+ fb = shapefile.Reader(freeBorderShapefile)
+ fbShapes = fb.shapes()
+ maxbbox=[1.e30, 1.e30, -1.e30, -1.e30]
+ outerBboxIndex = -1
+ for i,s in enumerate(fbShapes):
+ bb = s.bbox
+ if (bb[0] < maxbbox[0] and bb[1] < maxbbox[1] and bb[2] > maxbbox[2] and bb[3] > maxbbox[3]):
+ maxbbox = bb
+ outerBboxIndex = i
+ fbs = fbShapes[outerBboxIndex] # the domain free border shape
-# ---
+ # --- find the intersections of the shapefile to adjust and the domain free border:
+ # the closests points (two pairs of points)
+
+ sf = shapefile.Reader(shapefileToAdjust)
+ shapes = sf.shapes()
+ sfta = sf.shape(0)
+ pdist =[]
+ x0 = -1.e30
+ y0 = -1.e30
+ discountLastSftaPoint = 0
+ for i,p in enumerate(sfta.points):
+ if i == 0:
+ x00 = p[0] # keep first point
+ y00 = p[1]
+ d = math.sqrt((x0-p[0])*(x0-p[0]) + (y0-p[1])*(y0-p[1])) # distance to previous point
+ x0 = p[0] # keep previous point
+ y0 = p[1]
+ if d < 1.e-5:
+ print("two consecutives points of shapefile To adjust are superposed, OK")
+ continue # do not take into account consecutive superposed points in shape to adjust
+ if i == len(sfta.points) -1:
+ d = math.sqrt((x00-p[0])*(x00-p[0]) + (y00-p[1])*(y00-p[1])) # distance between first and last point
+ if d < 1.e-5:
+ discountLastSftaPoint = 1
+ print("last point of shapefile To adjust is superposed to first point, last point discarded, OK")
+ continue # do not take into account last point if superposed to first point, in shape to adjust
+ dmin = 1.e30
+ jmin = -1
+ for j,pfb in enumerate(fbs.points):
+ d = math.sqrt((pfb[0]-p[0])*(pfb[0]-p[0]) + (pfb[1]-p[1])*(pfb[1]-p[1]))
+ if d < dmin:
+ dmin = d
+ jmin = j
+ pdist.append((dmin, jmin, i)) # distance, index in freeBorder, index in shapeToAdjust
+ pdist.sort() # the 2 closest points must be set on the mesh freeBorder
+ print(pdist)
+ i1 = min(pdist[0][2], pdist[1][2]) # index of first adjusted point in shapeToAdjust
+ i2 = max(pdist[0][2], pdist[1][2]) # index of second adjusted point in shapeToAdjust
+ if i1 == pdist[0][2]:
+ ifb1 = pdist[0][1] # index of first adjusted point in freeBorder
+ ifb2 = pdist[1][1] # index of second adjusted point in freeBorder
+ else:
+ ifb1 = pdist[1][1]
+ ifb2 = pdist[0][1]
+ print("i1, i2, len(sfta.points)", i1, i2, len(sfta.points))
+ print("ifb1, ifb2, len(fbs.points)", ifb1, ifb2, len(fbs.points))
+
+ # --- write the adusted shapefile: free border closest points replaced with corresponding points
+ # on the free border. two polylines, one inside the domain, one outside
+
+ if outputDirectory == "":
+ outputDirectory = os.path.dirname(shapefileToAdjust)
+ a = os.path.splitext(os.path.basename(shapefileToAdjust))
+ shapefileAdjusted = os.path.join(outputDirectory, a[0] + '_adj' + a[1])
+ chainName = a[0] + '_adj'
+
+ w = shapefile.Writer(shapefileAdjusted)
+ w.shapeType = 3
+ w.field('name', 'C')
+
+ if splitShapeAdjusted:
+ chaincoords = []
+ chaincoords.append(fbs.points[ifb1])
+ for i in range(i1+1, i2):
+ chaincoords.append(sfta.points[i])
+ chaincoords.append(fbs.points[ifb2])
+ w.line([chaincoords])
+ w.record(chainName + '_0')
+
+ chaincoords = []
+ chaincoords.append(fbs.points[ifb2])
+ if i2+1 < len(sfta.points):
+ for i in range(i2+1, len(sfta.points) -discountLastSftaPoint):
+ chaincoords.append(sfta.points[i])
+ for i in range(i1):
+ chaincoords.append(sfta.points[i])
+ chaincoords.append(fbs.points[ifb1])
+ w.line([chaincoords])
+ w.record(chainName + '_1')
+ else:
+ chaincoords = []
+ chaincoords.append(fbs.points[ifb1])
+ for i in range(i1+1, i2):
+ chaincoords.append(sfta.points[i])
+ chaincoords.append(fbs.points[ifb2])
+ if i2+1 < len(sfta.points):
+ for i in range(i2+1, len(sfta.points) -discountLastSftaPoint):
+ chaincoords.append(sfta.points[i])
+ for i in range(i1):
+ chaincoords.append(sfta.points[i])
+ if discountLastSftaPoint:
+ chaincoords.append(fbs.points[ifb1]) # close shape when first point if superposed with last
+ w.line([chaincoords])
+ w.record(chainName)
+
+ w.close()
+
+ if splitFreeBorder:
+ # write the free border split in two polylines (cut by the adjusted shapefile)
+
+ if outputDirectory == "":
+ outputDirectory = os.path.dirname(freeBorderShapefile)
+ a = os.path.splitext(os.path.basename(freeBorderShapefile))
+ freeBorderSplit = os.path.join(outputDirectory, a[0] + '_split' + a[1])
+ chainName = a[0] + '_split'
+
+ w = shapefile.Writer(freeBorderSplit)
+ w.shapeType = 3
+ w.field('name', 'C')
+
+ if (ifb1 > ifb2):
+ i = ifb1; ifb1 = ifb2; ifb2 = i
+
+ chaincoords = []
+ for i in range(ifb1, ifb2+1):
+ chaincoords.append(fbs.points[i])
+ w.line([chaincoords])
+ w.record(chainName + '_0')
+
+ chaincoords = []
+ for i in range(ifb2, len(fbs.points)):
+ chaincoords.append(fbs.points[i])
+ for i in range(ifb1+1):
+ chaincoords.append(fbs.points[i])
+ w.line([chaincoords])
+ w.record(chainName + '_1')
+ w.close()
-#meshFile = freeBordersGroup('/home/paul/projets/hydro95/V9_5_BR/tests/Maill_SB1610_ALB1501_barr_eff.med')
-#exploreEdgeGroups(meshFile, 1000000., 6000000.)
