arcgis-10.0arcgis-enginearcpysnapping
I have a layer containing a vector 'Line 1', that originates at point 'A' on a second layer. I also have point 'B' on a third layer. I want to create 'Line 2' that has the same angle and length of 'Line 1', but originates from point 'B'. And I want to repeat this for all vectors on the layer containing 'Line 1'.
i'm working with ArcGIS Engine 10 for .NET. I'm open to solutions working with just ArcGIS 10, so I can create the same geoprocessing tool in ArcGIS Engine.
I've managed to solve this, without using the mentioned GRASS tools or topological functions.
Basically I take all start- and endnodes, put them in a new, temporary table, put a buffer around them, union the buffer objects, and move all found nodes in each buffer to the centroid of the buffer.
When that's done I move the original begin and end points to the new location.
Easier than expected, and still fast, but I expected PostGIS to have some built-in function for this - that would be even quicker.
Edit: in the interest of giving back to the community, this is my (quite crappy) code for now.
drop table if exists nodes;
drop table if exists nodes2;
drop table if exists buffers;
-- Get Start and End nodes
select ST_StartPoint(wkb_geometry) startnode, ST_EndPoint(wkb_geometry) endnode, ogc_fid into nodes from sourceTable;
-- Combine all nodes into one table for easier queries
select startnode node, ogc_fid into nodes2 from nodes;
insert into nodes2 select endnode node, ogc_fid from nodes;
-- Some indexes to speed everything up
CREATE INDEX nodesstart_idx ON nodes USING gist (startnode);
CREATE INDEX nodesend_idx ON nodes USING gist (endnode);
CREATE INDEX nodes2_idx ON nodes2 USING gist (node);
CREATE INDEX nodes_ogcfid_idx ON nodes USING btree (ogc_fid ASC NULLS LAST);
-- Create buffers, combine them, split combined objects again
select (ST_Dump(ST_Union(ST_Buffer(node, 1)))).geom geom into buffers from nodes2;
CREATE INDEX buffers_idx ON buffers USING gist (geom);
-- Update start/end nodes table
UPDATE nodes SET startnode = ST_Centroid((select geom from buffers WHERE geom && startnode));
UPDATE nodes SET endnode = ST_Centroid((select geom from buffers WHERE geom && endnode));
-- Update original points
update sourceTable set wkb_geometry = ST_SetPoint(
ST_SetPoint(wkb_geometry, 0, (select startnode from nodes where ogc_fid=sourceTable.ogc_fid)),
ST_NumPoints(wkb_geometry) - 1, (select endnode from nodes where ogc_fid=sourceTable.ogc_fid));
DROP TABLE nodes;
DROP TABLE nodes2;
DROP TABLE buffers;
Here is a solution for MultiLineStrings and MultiPoints. It moves the start- and endpoints of every part of the MultiLineString to the closest point.
In case you do not want to move start- and enpoints of every part, but only of every feature, you may want to change these two lines accordingly (simply comment/uncomment in the code below):
if part.startPoint() == geom.vertexAt(vert): to if feat.geometry().vertexAt(0) == geom.vertexAt(vert):elif part.endPoint() == geom.vertexAt(vert): to elif feat.geometry().vertexAt([i for i, f in enumerate(feat.geometry().vertices())][-1]) == geom.vertexAt(vert): (not very convenient, but simply doing feat.geometry().vertexAt(-1) directly doesnt work...)Here is the full code with further explanations as comments:
points = QgsProject.instance().mapLayersByName('points')[0]
lines = QgsProject.instance().mapLayersByName('lines')[0]
tolerance = 10000 # snapping tolerance in CRS units
if QgsWkbTypes.isMultiType(points.wkbType()): # if point-input is of type multigeometry
res = processing.run("native:multiparttosingleparts",{'INPUT':points,'OUTPUT':'TEMPORARY_OUTPUT'}) # convert to singlepoint
singlepoints = res['OUTPUT']
else: # if point-input is of type singlegeometry
singlepoints = points # leave the inputlayer as it is
# create the spatial index of the pointlayer
points_idx = QgsSpatialIndex(singlepoints.getFeatures())
with edit(lines):
for feat in lines.getFeatures():
geom = feat.geometry()
vert = 0 # reset vertices-count on every new feature
vert_dict = {} # clear vertices dictionary on every new feature
for part in geom.parts(): # iterate through all parts of each feature
for vertex in part.vertices(): # iterate through all vertices of each part
# if feat.geometry().vertexAt(0) == geom.vertexAt(vert): # if the current vertex is a startpoint of the feature
if part.startPoint() == geom.vertexAt(vert): # if the current vertex is a startpoint of the features part, then:
nearestneighbor_start = points_idx.nearestNeighbor(QgsPointXY(part.startPoint()), neighbors=1, maxDistance=tolerance) # find the closest point to the current startpoint
if len(nearestneighbor_start) > 0:# only adjust the startpoint if there is a nearest point within the maxdistance
nearpoint_startgeom = singlepoints.getFeature(nearestneighbor_start[0]).geometry() # get the geometry of the nearest point in pointlayer
vert_dict[vert] = nearpoint_startgeom.asPoint() # add the index of the current vertex as key and the geometry of the closest point as value to the dictionary
# elif feat.geometry().vertexAt([i for i, f in enumerate(feat.geometry().vertices())][-1]) == geom.vertexAt(vert): # if the current vertex is an endpoint of the current feature
elif part.endPoint() == geom.vertexAt(vert): # if the current vertex is an endpoint of the current features part
nearestneighbor_end = points_idx.nearestNeighbor(QgsPointXY(part.endPoint()), neighbors=1, maxDistance=tolerance)
if len(nearestneighbor_end) > 0:
nearpoint_endgeom = singlepoints.getFeature(nearestneighbor_end[0]).geometry()
vert_dict[vert] = nearpoint_endgeom.asPoint()
else:
pass # if current vertex is not a start or end point skip it...
vert += 1 # increase vertices-counter
for vertindex, newpoint in vert_dict.items(): # for every feature iterate over the just created dictionary (vertindex (=dict key) is the start or endpoint we want to move and newpoint (=dict value) the position we want to move it to)
QgsVectorLayerEditUtils(lines).moveVertex(newpoint.x(),newpoint.y(),feat.id(),vertindex) # use QgsVectorLayerEditUtils to edit the vertex, for more details see https://qgis.org/pyqgis/3.4/core/QgsVectorLayerEditUtils.html#qgis.core.QgsVectorLayerEditUtils
# https://gis.stackexchange.com/questions/410768/snapping-lines-start-point-end-points-to-the-closest-point-features-to-them/410778#410778
For further background you may want to read Modifying specific vertices of MultiLineString using PyQGIS
Example for snapping start- and endpoints of all parts of all feature:
Example for snapping start- and enpoints only for all features:
Best Answer
I have made an ArcPy code to virtually snap the lines, so the script calculates the bearing of the lines, snaps the origin of the vector, then build new lines from new snapped points using the bearing.
But I have noticed a weird error using the function
BearingDistanceToLine_managementin ArcPy. It adds an angle to the bearings. Without understanding the cause of the error, I have eliminated the difference by180 - 57.935