+ # --- 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()