Intersection Point between 3D Line and 3D polygon

Hi,

The LineOnAreaOverlayer does work on 3D features and so does the TopologyBuilder but you need to overlay the resulting node with the original line to force the node (point) into 3D.

Hi @faifai1214, I think you will have to perform geometric operations to achieve that. This is an experimental prototype implemented with Python 2.7 script. Hope this helps.

3d-intersection-experiment.fmw (FME 2017.0.1)

# PythonCaller Script Example (Python 2.7 only)# Create the intersection point between a plane containing the first three vertices# of 3D polygon and a straight line containing the first segment of 3D line.# Note: The resulting intersection could be outside of the source polygon# and/or the source line segment.import fmeobjectsclass FeatureProcessor(object):    def __init__(self):        self.polygons = []        self.lines = []            def input(self, feature):        # Collect 3D polygons and 3D lines.        if feature.getDimension() == fmeobjects.FME_THREE_D:            geom = feature.getGeometry()            if isinstance(geom, fmeobjects.FMEPolygon):                self.polygons.append(feature)            elif isinstance(geom, fmeobjects.FMELine):                self.lines.append(feature)                    def close(self):        # Add ID to source polygons and lines, then output them.        for i, area in enumerate(self.polygons):            area.setAttribute('_area_id', i)            self.pyoutput(area)        for i, line in enumerate(self.lines):            line.setAttribute('_line_id', i)            self.pyoutput(line)                    # Create and output intersections for each pair of plane and line.        for area in self.polygons:            p = [area.getCoordinate(i) for i in range(3)] # First three vertices.            area_id = area.getAttribute('_area_id')            for line in self.lines:                q = [line.getCoordinate(i) for i in range(2)] # First two vertices.                coord = intersection(p[0], p[1], p[2], q[0], q[1])                if coord != None:                    point = fmeobjects.FMEFeature()                    point.setGeometry(fmeobjects.FMEPoint(*coord))                    point.setAttribute('_area_id', area_id)                    point.setAttribute('_line_id', line.getAttribute('_line_id'))                    self.pyoutput(point)                    # Return coordinates of the intersection between a plane and a straight line in 3D.# Parameters: 3 coordinates specify the plane, 2 coordinates specify the line.def intersection((xp0, yp0, zp0), (xp1, yp1, zp1), (xp2, yp2, zp2),    (xl0, yl0, zl0), (xl1, yl1, zl1)):    # Normal vector of the plane.    x0, y0, z0 = float(xp1 - xp0), float(yp1 - yp0), float(zp1 - zp0)    x1, y1, z1 = float(xp2 - xp0), float(yp2 - yp0), float(zp2 - zp0)    nx, ny, nz = y0 * z1 - z0 * y1, z0 * x1 - x0 * z1, x0 * y1 - y0 * x1        # Direction (vector) of the line.    vx, vy, vz = float(xl1- xl0), float(yl1 - yl0), float(zl1 - zl0)    # Calculate coordinates of the intersection.        d = nx * vx + ny * vy + nz * vz    if d != 0:        t = (nx * (xp0- xl0) + ny * (yp0 - yl0) + nz * (zp0 - zl0)) / d        return (xl0 + vx * t, yl0 + vy * t, zl0 + vz * t)    else:        return None

@faifai1214 I needed to dissect ~10k complex 3d lines where they intersect 8 complex 3d planes. My 3d planes started off as xyz grids which I converted to DEM Raster’s using SurfaceModeller. I converted my 3d lines to 3D points using LineDivider with a spacing of 1ft increments. Once I had the 3D points and a DEM Raster I perform a PointOnRasterOverlay which gives me the raster cell value directly under each point. With this information I can then calculate the difference between each point and it’s corresponding raster cell value. Weather the calculated difference is positive or negative tells me if that point is above or below the 3d surface and I can assign an attribute to identify that relationship. I compared each point to one surface at a time because the lines represented by the series of points were complex. A line could intersect one, many, all or none of the planes. A line could also intersect a plane multiple times. I had to use a conditional statement and the adjacent feature attributes option to assign a segment number or grouping to each point which could be used later to help reassemble the points into lines using LineBuilder. In LineBuilder you specify a connection break attribute to define where to start and stop each line assembly. In fact I had to specify three attributes to define the segment grouping (Line ID, Segment, Zone). In the finished process my ~10k lines were dissected into ~60k segments. This process is a bit imprecise because I divided my lines into points at 1ft increments. When I reassemble the points into lines based on LineID, Segement, and Zone, I end up with a small gap of 1 ft between each dissected line. In the overall scale 1ft gaps will not bother me and in fact I can’t even see them in my maps. Your coordinate system units and LineDivider length parameter determine your final resolution or gap. The process is very intensive and I continue to look for more efficient approaches. A LineOnAreaOverlayer which works in 3D might be very valuable. Am I overlooking a more efficient or more precise approach?