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, outputDirectory="", 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     outputDirectory: if empty, write the resulting shapefiles in their respective directory, with the suffix '_adj', otherwise write in the outputDirectory.
 234     splitFreeBorder: boolean default True
 235     splitShapeAdjusted: boolean default True
 236     """
 237     print(" === fitShapePointsToMesh", freeBorderShapefile, shapefileToAdjust)
 238
 239     # --- find domain freeBorder:  bounding box englobing all others
 240     #     TODO: This may not be always the case, when there is not a single domain with holes, but several non connected domains.
 241
 242     fb = shapefile.Reader(freeBorderShapefile)
 243     fbShapes = fb.shapes()
 244     maxbbox=[1.e30, 1.e30, -1.e30, -1.e30]
 245     outerBboxIndex = -1
 246     for i,s in enumerate(fbShapes):
 247         bb = s.bbox
 248         if (bb[0] < maxbbox[0] and bb[1] < maxbbox[1] and bb[2] > maxbbox[2] and bb[3] > maxbbox[3]):
 249             maxbbox = bb
 250             outerBboxIndex = i
 251     fbs = fbShapes[outerBboxIndex] # the domain free border shape
 252
 253     # --- find the intersections of the shapefile to adjust and the domain free border:
 254     #     the closests points (two pairs of points)
 255
 256     sf = shapefile.Reader(shapefileToAdjust)
 257     shapes = sf.shapes()
 258     sfta = sf.shape(0)
 259     pdist =[]
 260     x0 = -1.e30
 261     y0 = -1.e30
 262     discountLastSftaPoint = 0
 263     for i,p in enumerate(sfta.points):
 264         if i == 0:
 265             x00 = p[0] # keep first point
 266             y00 = p[1]
 267         d = math.sqrt((x0-p[0])*(x0-p[0]) + (y0-p[1])*(y0-p[1])) # distance to previous point
 268         x0 = p[0] # keep previous point
 269         y0 = p[1]
 270         if d < 1.e-5:
 271             print("two consecutives points of shapefile To adjust are superposed, OK")
 272             continue # do not take into account consecutive superposed points in shape to adjust
 273         if i == len(sfta.points) -1:
 274             d = math.sqrt((x00-p[0])*(x00-p[0]) + (y00-p[1])*(y00-p[1])) # distance between first and last point
 275             if d < 1.e-5:
 276                 discountLastSftaPoint = 1
 277                 print("last point of shapefile To adjust is superposed to first point, last point discarded, OK")
 278                 continue # do not take into account last point if superposed to first point, in shape to adjust
 279         dmin = 1.e30
 280         jmin = -1
 281         for j,pfb in enumerate(fbs.points):
 282             d = math.sqrt((pfb[0]-p[0])*(pfb[0]-p[0]) + (pfb[1]-p[1])*(pfb[1]-p[1]))
 283             if d < dmin:
 284                 dmin = d
 285                 jmin = j
 286         pdist.append((dmin, jmin, i)) # distance, index in freeBorder, index in shapeToAdjust
 287     pdist.sort() # the 2 closest points must be set on the mesh freeBorder
 288     print(pdist)
 289     i1 = min(pdist[0][2], pdist[1][2]) # index of first adjusted point in shapeToAdjust
 290     i2 = max(pdist[0][2], pdist[1][2]) # index of second adjusted point in shapeToAdjust
 291     if i1 == pdist[0][2]:
 292         ifb1 = pdist[0][1]             # index of first adjusted point in freeBorder
 293         ifb2 = pdist[1][1]             # index of second adjusted point in freeBorder
 294     else:
 295         ifb1 = pdist[1][1]
 296         ifb2 = pdist[0][1]
 297     print("i1, i2, len(sfta.points)", i1, i2, len(sfta.points))
 298     print("ifb1, ifb2, len(fbs.points)", ifb1, ifb2, len(fbs.points))
 299
 300     # --- write the adusted shapefile: free border closest points replaced with corresponding points
 301     #     on the free border. two polylines, one inside the domain, one outside
 302
 303     if outputDirectory == "":
 304         outputDirectory = os.path.dirname(shapefileToAdjust)
 305     a = os.path.splitext(os.path.basename(shapefileToAdjust))
 306     shapefileAdjusted = os.path.join(outputDirectory, a[0] + '_adj' + a[1])
 307     chainName = a[0]
 308
 309     w = shapefile.Writer(shapefileAdjusted)
 310     w.shapeType = 3
 311     w.field('name', 'C')
 312
 313     if splitShapeAdjusted:
 314         chaincoords = []
 315         chaincoords.append(fbs.points[ifb1])
 316         for i in range(i1+1, i2):
 317             chaincoords.append(sfta.points[i])
 318         chaincoords.append(fbs.points[ifb2])
 319         w.line([chaincoords])
 320         w.record(chainName + '_0')
 321
 322         chaincoords = []
 323         chaincoords.append(fbs.points[ifb2])
 324         if i2+1 < len(sfta.points):
 325             for i in range(i2+1, len(sfta.points) -discountLastSftaPoint):
 326                 chaincoords.append(sfta.points[i])
 327         for i in range(i1):
 328             chaincoords.append(sfta.points[i])
 329         chaincoords.append(fbs.points[ifb1])
 330         w.line([chaincoords])
 331         w.record(chainName + '_1')
 332     else:
 333         chaincoords = []
 334         chaincoords.append(fbs.points[ifb1])
 335         for i in range(i1+1, i2):
 336             chaincoords.append(sfta.points[i])
 337         chaincoords.append(fbs.points[ifb2])
 338         if i2+1 < len(sfta.points):
 339             for i in range(i2+1, len(sfta.points) -discountLastSftaPoint):
 340                 chaincoords.append(sfta.points[i])
 341         for i in range(i1):
 342             chaincoords.append(sfta.points[i])
 343         if discountLastSftaPoint:
 344             chaincoords.append(fbs.points[ifb1]) # close shape when first point if superposed with last
 345         w.line([chaincoords])
 346         w.record(chainName)
 347
 348     w.close()
 349
 350     if splitFreeBorder:
 351         # write the free border split in two polylines (cut by the adjusted shapefile)
 352
 353         if outputDirectory == "":
 354             outputDirectory = os.path.dirname(freeBorderShapefile)
 355         a = os.path.splitext(os.path.basename(freeBorderShapefile))
 356         freeBorderSplit = os.path.join(outputDirectory, a[0] + '_split' + a[1])
 357         chainName = a[0]
 358
 359         w = shapefile.Writer(freeBorderSplit)
 360         w.shapeType = 3
 361         w.field('name', 'C')
 362
 363         if (ifb1 > ifb2):
 364             i = ifb1; ifb1 = ifb2; ifb2 = i
 365
 366         chaincoords = []
 367         for i in range(ifb1, ifb2+1):
 368             chaincoords.append(fbs.points[i])
 369         w.line([chaincoords])
 370         w.record(chainName + '_0')
 371
 372         chaincoords = []
 373         for i in range(ifb2, len(fbs.points)):
 374             chaincoords.append(fbs.points[i])
 375         for i in range(ifb1+1):
 376             chaincoords.append(fbs.points[i])
 377         w.line([chaincoords])
 378         w.record(chainName + '_1')
 379         w.close()
 380