src/HYDROTools/shapesGroups.py

   1
   2 import salome
   3 salome.salome_init()
   4
   5 import  SMESH, SALOMEDS
   6 from salome.smesh import smeshBuilder
   7
   8 import numpy as np
   9 import MEDLoader as ml
  10 import medcoupling as mc
  11
  12 import shapefile
  13 import math
  14 import os
  15
  16
  17 def freeBordersGroup(meshFileIn, meshFileOut=""):
  18     """
  19     In a mesh, create a group of Edges for the borders: domain, internal holes (isles)
  20     parameters:
  21     meshFileIn: full path of the input mesh file, format MED
  22     meshFileOut: full path of the output mesh file, format MED (default="" : when "", output file is suffixed with "_brd.med"
  23     return full path of the output mesh file
  24     """
  25     smesh = smeshBuilder.New()
  26     if not salome.sg.hasDesktop():
  27         smesh.SetEnablePublish( False ) # Set to False to avoid publish in study if not needed
  28     ([MESH], status) = smesh.CreateMeshesFromMED(meshFileIn)
  29
  30     nbAdded, MESH, addedBnd = MESH.MakeBoundaryElements( SMESH.BND_1DFROM2D, '', '', 0, [])
  31
  32     aCriteria = []
  33     aCriterion = smesh.GetCriterion(SMESH.EDGE,SMESH.FT_FreeBorders,SMESH.FT_Undefined,0)
  34     aCriteria.append(aCriterion)
  35     aFilter = smesh.GetFilterFromCriteria(aCriteria)
  36     aFilter.SetMesh(MESH.GetMesh())
  37     FreeBorders = MESH.GroupOnFilter( SMESH.EDGE, 'FreeBorders', aFilter )
  38
  39     a = os.path.splitext(meshFileIn)
  40     smesh.SetName(MESH, os.path.basename(a[0]))
  41
  42     if meshFileOut == "":
  43         a = os.path.splitext(meshFileIn)
  44         smesh.SetName(MESH, os.path.basename(a[0]))
  45         meshFileOut = a[0] + '_brd' + a[1]
  46
  47     try:
  48         MESH.ExportMED(meshFileOut,auto_groups=0,minor=40,overwrite=1,meshPart=None,autoDimension=1)
  49         pass
  50     except:
  51         print('ExportMED() failed. Invalid file name?')
  52     return meshFileOut
  53
  54
  55 def explodeGroup(grp, grpName):
  56     """
  57     from a group of edges loaded with MEDCoupling, create ordered lists of nodes, one list for each set of connected edges.
  58     parameters:
  59     grp: MEDCoupling object for the group of edges
  60     grpName: name of the group
  61     return:
  62     List of descriptors [(ordered list of nodeIds, name of the list, closed status)]
  63     """
  64     print(" === explodeGroup", grpName)
  65     nbCells=grp.getNumberOfCells()
  66
  67     dicReverse = {} # id node --> id edges
  68     for i in range(nbCells):
  69         nodcell = grp.getNodeIdsOfCell(i)
  70         for j in range(len(nodcell)):
  71             if nodcell[j] in dicReverse:
  72                 dicReverse[nodcell[j]].append(i)
  73             else:
  74                 dicReverse[nodcell[j]] = [i]
  75
  76     nodeChains = []
  77     usedCells = [False] * nbCells
  78     while False in usedCells:
  79         icell = usedCells.index(False)
  80         usedCells[icell] = True
  81         nodcell = grp.getNodeIdsOfCell(icell)
  82         closed = False
  83         chain = [nodcell[0], nodcell[1]]
  84         nextnode = nodcell[1]
  85         prevnode = nodcell[0]
  86         while nextnode in dicReverse:
  87             nextcells = dicReverse[nextnode]
  88             if len(nextcells) != 2:             # end of chain(1) or "edges connector"(>2): stop
  89                 closed = False
  90                 nextnode = -1                   # stop
  91             else:
  92                 newcell =False
  93                 for i in range(len(nextcells)):
  94                     ncell = nextcells[i]
  95                     if not usedCells[ncell]:       # the chain of nodes grows
  96                         usedCells[ncell] = True
  97                         newcell = True
  98                         nodcell = grp.getNodeIdsOfCell(ncell)
  99                         if nodcell[0] == nextnode:
 100                             nextnode = nodcell[1]  # forward edge
 101                         else:
 102                             nextnode = nodcell[0]  # reversed edge ?
 103                         chain.append(nextnode)
 104                 if not newcell:                    # end of chain, closed
 105                     closed =True
 106                     nextnode = -1
 107         while prevnode in dicReverse:
 108             prevcells = dicReverse[prevnode]
 109             if len(prevcells) != 2:             # end of chain(1) or "edges connector"(>2): stop
 110                 closed = False
 111                 prevnode = -1                   # stop
 112             else:
 113                 newcell =False
 114                 for i in range(len(prevcells)):
 115                     ncell = prevcells[i]
 116                     if not usedCells[ncell]:       # the chain of nodes grows
 117                         usedCells[ncell] = True
 118                         newcell = True
 119                         nodcell = grp.getNodeIdsOfCell(ncell)
 120                         if nodcell[1] == prevnode:
 121                             prevnode = nodcell[0]  # forward edge
 122                         else:
 123                             prevnode = nodcell[1]  # reversed edge ?
 124                         chain.insert(0, prevnode)
 125                 if not newcell:                    # end of chain, closed
 126                     closed =True
 127                     prevnode = -1
 128
 129         chainDesc = (chain, grpName +"_%s" % len(nodeChains), closed)
 130         nodeChains.append(chainDesc)
 131         print(chainDesc[1:])
 132     return nodeChains
 133
 134
 135 def writeShapeLines(mcMesh, grpName, nodeChains, outputDirectory, offsetX=0., offsetY=0.):
 136     """
 137     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.
 138     parameters:
 139     mcMesh: mesh loaded in memory with MEDLoader
 140     grpName: name associated to the group of edges
 141     nodeChains: List of descriptors corresponding to the group [(ordered list of nodeIds, name of the list, closed status)]
 142     outputDirectory: directory for writing the shapefile
 143     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.
 144     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.
 145     """
 146     print(" === writeShapeLines", grpName)
 147     coords = mcMesh.getCoords()
 148     shapeFileName = os.path.join(outputDirectory, grpName)
 149     w = shapefile.Writer(shapeFileName)
 150     w.shapeType = 3
 151     w.field('name', 'C')
 152     for (chain, chainName, closed) in nodeChains:
 153         print("   --- ", chainName)
 154         chaincoords = []
 155         for node in chain:
 156             coord = coords[node].getValues()
 157             coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
 158             chaincoords.append(coordLb93)
 159         w.line([chaincoords])
 160         w.record(chainName)
 161     w.close()
 162
 163
 164 def writeShapePoints(mcMesh, grpName, nodeChains, outputDirectory, offsetX=0., offsetY=0.):
 165     """
 166     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.
 167     parameters:
 168     mcMesh: mesh loaded in memory with MEDLoader
 169     grpName: name associated to the group of edges
 170     nodeChains: List of descriptors corresponding to the group [(ordered list of nodeIds, name of the list, closed status)]
 171     outputDirectory: directory for writing the shapefile
 172     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.
 173     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.
 174     """
 175     print(" === writeShapePoints", grpName)
 176     coords = mcMesh.getCoords()
 177     shapeFileName = os.path.join(outputDirectory, grpName)
 178     w = shapefile.Writer(shapeFileName + '_pts')
 179     #w.shapeType = 8
 180     w.field('name', 'C')
 181     for (chain, chainName, closed) in nodeChains:
 182         print("   --- ", chainName)
 183         chaincoords = []
 184         for node in chain:
 185             coord = coords[node].getValues()
 186             coordLb93=[coord[0] + offsetX, coord[1] + offsetY]
 187             chaincoords.append(coordLb93)
 188         w.multipoint(chaincoords)
 189         w.record(chainName)
 190     w.close()
 191
 192
 193 def exploreEdgeGroups(meshFile, outputDirectory="", offsetX=0., offsetY=0.):
 194     """
 195     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)
 196     parameters:
 197     meshFile: full path of the input mesh file, format MED
 198     outputDirectory: directory in which the shapefiles are written (default "", if "" use the directory containing the mesh
 199     offsetX: local X origin of the mesh
 200     offsetY: local Y origin of the mesh
 201     """
 202     print(" === exploreEdgeGroups", meshFile)
 203     if outputDirectory == "":
 204         outputDirectory = os.path.dirname(meshFile)
 205
 206     a = os.path.splitext(meshFile)
 207     prefix = os.path.basename(a[0]) # prefix = file name without extension
 208
 209     mcMesh = ml.MEDFileMesh.New(meshFile)
 210     dim = mcMesh.getMeshDimension()
 211     d1=-1        # when dimension 2, edges are dim -1
 212     if dim == 3: # when dimension 3, edges are dim -2
 213         d1=-2
 214
 215     grp_names = mcMesh.getGroupsOnSpecifiedLev(d1) #names of edges groups
 216
 217     groups = [mcMesh.getGroup(d1, name) for name in grp_names] # list of groups in their name order
 218
 219     for (grp, grpName) in zip(groups, grp_names):
 220         fullGrpName = prefix + '_' + grpName
 221         nodeChains = explodeGroup(grp, fullGrpName)
 222         writeShapeLines(mcMesh, fullGrpName, nodeChains, outputDirectory, offsetX, offsetY)
 223         writeShapePoints(mcMesh, fullGrpName, nodeChains, outputDirectory, offsetX, offsetY)
 224
 225
 226 def fitShapePointsToMesh(freeBorderShapefile, shapefileToAdjust, splitFreeBorder=True, splitShapeAdjusted=True):
 227     """
 228     shapeFileToAdjust must be a closed line or polygon crossing freeBorderShapefile in 2 points.
 229     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.
 230     parameters:
 231     freeBorderShapefile: a set of free border lines, as generated by the functions freeBordersGroup and exploreEdgeGroups.
 232     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.
 233     splitFreeBorder: boolean default True
 234     splitShapeAdjusted: boolean default True
 235     """
 236     print(" === fitShapePointsToMesh", freeBorderShapefile, shapefileToAdjust)
 237
 238     # --- find domain freeBorder:  bounding box englobing all others
 239     #     TODO: This may not be always the case, when there is not a single domain with holes, but several non connected domains.
 240
 241     fb = shapefile.Reader(freeBorderShapefile)
 242     fbShapes = fb.shapes()
 243     maxbbox=[1.e30, 1.e30, -1.e30, -1.e30]
 244     outerBboxIndex = -1
 245     for i,s in enumerate(fbShapes):
 246         bb = s.bbox
 247         if (bb[0] < maxbbox[0] and bb[1] < maxbbox[1] and bb[2] > maxbbox[2] and bb[3] > maxbbox[3]):
 248             maxbbox = bb
 249             outerBboxIndex = i
 250     fbs = fbShapes[outerBboxIndex] # the domain free border shape
 251
 252     # --- find the intersections of the shapefile to adjust and the domain free border:
 253     #     the closests points (two pairs of points)
 254
 255     sf = shapefile.Reader(shapefileToAdjust)
 256     shapes = sf.shapes()
 257     sfta = sf.shape(0)
 258     pdist =[]
 259     x0 = -1.e30
 260     y0 = -1.e30
 261     for i,p in enumerate(sfta.points):
 262         d = math.sqrt((x0-p[0])*(x0-p[0]) + (y0-p[1])*(y0-p[1]))
 263         x0 = p[0]
 264         y0 = p[1]
 265         if d < 1.e-5:
 266             continue # do not take into account superposed points in shape to adjust
 267         dmin = 1.e30
 268         jmin = -1
 269         for j,pfb in enumerate(fbs.points):
 270             d = math.sqrt((pfb[0]-p[0])*(pfb[0]-p[0]) + (pfb[1]-p[1])*(pfb[1]-p[1]))
 271             if d < dmin:
 272                 dmin = d
 273                 jmin = j
 274         pdist.append((dmin, jmin, i)) # distance, index in freeBorder, index in shapeToAdjust
 275     pdist.sort() # the 2 closest points must be set on the mesh freeBorder
 276     print(pdist)
 277     i1 = min(pdist[0][2], pdist[1][2]) # index of first adjusted point in shapeToAdjust
 278     i2 = max(pdist[0][2], pdist[1][2]) # index of second adjusted point in shapeToAdjust
 279     if i1 == pdist[0][2]:
 280         ifb1 = pdist[0][1]             # index of first adjusted point in freeBorder
 281         ifb2 = pdist[1][1]             # index of second adjusted point in freeBorder
 282     else:
 283         ifb1 = pdist[1][1]
 284         ifb2 = pdist[0][1]
 285     print("i1, i2, len(sfta.points)", i1, i2, len(sfta.points))
 286     print("ifb1, ifb2, len(fbs.points)", ifb1, ifb2, len(fbs.points))
 287
 288     # --- write the adusted shapefile: free border closest points replaced with corresponding points
 289     #     on the free border. two polylines, one inside the domain, one outside
 290
 291     a = os.path.splitext(shapefileToAdjust)
 292     shapefileAdjusted = a[0] + '_adj' + a[1]
 293     chainName = os.path.basename(a[0])
 294
 295     w = shapefile.Writer(shapefileAdjusted)
 296     w.shapeType = 3
 297     w.field('name', 'C')
 298
 299     if splitShapeAdjusted:
 300         chaincoords = []
 301         chaincoords.append(fbs.points[ifb1])
 302         for i in range(i1+1, i2):
 303             chaincoords.append(sfta.points[i])
 304         chaincoords.append(fbs.points[ifb2])
 305         w.line([chaincoords])
 306         w.record(chainName + '_0')
 307
 308         chaincoords = []
 309         chaincoords.append(fbs.points[ifb2])
 310         if i2+1 < len(sfta.points):
 311             for i in range(i2+1, len(sfta.points)):
 312                 chaincoords.append(sfta.points[i])
 313         for i in range(i1):
 314             chaincoords.append(sfta.points[i])
 315         chaincoords.append(fbs.points[ifb1])
 316         w.line([chaincoords])
 317         w.record(chainName + '_1')
 318     else:
 319         chaincoords = []
 320         chaincoords.append(fbs.points[ifb1])
 321         for i in range(i1+1, i2):
 322             chaincoords.append(sfta.points[i])
 323         chaincoords.append(fbs.points[ifb2])
 324         if i2+1 < len(sfta.points):
 325             for i in range(i2+1, len(sfta.points)):
 326                 chaincoords.append(sfta.points[i])
 327         for i in range(i1):
 328             chaincoords.append(sfta.points[i])
 329         w.line([chaincoords])
 330         w.record(chainName)
 331
 332     w.close()
 333
 334     if splitFreeBorder:
 335         # write the free border split in two polylines (cut by the adjusted shapefile)
 336
 337         a = os.path.splitext(freeBorderShapefile)
 338         freeBorderSplit = a[0] + '_split' + a[1]
 339         chainName = os.path.basename(a[0])
 340
 341         w = shapefile.Writer(freeBorderSplit)
 342         w.shapeType = 3
 343         w.field('name', 'C')
 344
 345         if (ifb1 > ifb2):
 346             i = ifb1; ifb1 = ifb2; ifb2 = i
 347
 348         chaincoords = []
 349         for i in range(ifb1, ifb2+1):
 350             chaincoords.append(fbs.points[i])
 351         w.line([chaincoords])
 352         w.record(chainName + '_0')
 353
 354         chaincoords = []
 355         for i in range(ifb2, len(fbs.points)):
 356             chaincoords.append(fbs.points[i])
 357         for i in range(ifb1+1):
 358             chaincoords.append(fbs.points[i])
 359         w.line([chaincoords])
 360         w.record(chainName + '_1')
 361         w.close()
 362