Split Oracle Sdo_Geometry Linestring at a Known Point using PL/SQL (1)
Splitting a LineString at a known point
Superceded
The implementation in this article has been superseded since the SDO_LRS package being made available to non-Spatial licensed users (see Oracle statements about this to ensure you are able to use SDO_LRS).
Another article describes how to build a function called ST_Snap based on SDO_LRS procedures and functions.
Background
A correspondent emailed me asking:
Just wondering if you could shed some light on a problem I have. In oracle using PLSQL I need to cut a polyline at two know (sic) points (not necessarly vertices and the line may be not be straight) and copy the cut line into another table. I would prefer not to go to linear referencing to do this. I cannot seem to find a function that will cut polylines. I can find the concatinate (sic) functions but this does not help much.
I’m new to oracal (sic) spatial but have a very extensive GIS background mostly ESRI and was suprised that there was not a built in function to split a line.
Any help you could offer would be greatly appreciated and really help me out.
Response
It is true that there there is no function in Oracle Locator/Spatial that can split a linestring at a known point and return the halves as separate linestrings. The closest LRS function in LRS is:
SDO_LRS.SPLIT_GEOM_SEGMENT(
geom_segment IN SDO_GEOMETRY,
split_measure IN NUMBER,
segment_1 OUT SDO_GEOMETRY,
segment_2 OUT SDO_GEOMETRY);
But note that it does not take a point sdo_geometry but, surprise, surprise, a measure! (Also it is a procedure and not a function.)
NOTE: Later versions of Oracle Spatial provide free access to SDO_LRS which removes the need for this implementation.
My correspondent happened to be lucky as I had just been playing around with some of the vector functions in my GEOM PL/SQL package to find the nearest vector (segment of a linestring) closest to other objects eg:
select rownum as id,
b.startcoord.x as x1,
b.startcoord.y as y1,
b.endcoord.x as x2,
b.endcoord.y as y2,
sdo_geom.sdo_distance(mdsys.sdo_geometry(2002,NULL,NULL,
mdsys.sdo_elem_info_array(1,2,1),
mdsys.sdo_ordinate_array(b.startcoord.x,b.startcoord.y,b.endcoord.x,b.endcoord.y)),
mdsys.sdo_geometry(2001,null,mdsys.sdo_point_type(380326.792,5167489.29,NULL),NULL,NULL),
0.005) as linedist
from projline2d a,
table(codesys.geom.getpipedvector2d(a.geom)) b
where linetype = 'VERTEX'
/
ID X1 Y1 X2 Y2 LINEDIST
---------- ---------- ---------- ---------- ---------- ----------
1 380326.792 5167089.29 380326.792 5167889.29 0
2 380326.792 5167889.29 380826.792 5167889.29 399.996
3 380826.792 5167889.29 380126.792 5167489.29 99.2243147
I realised that this was the basis of what my correspondent wanted so I created a procedure called Split and have integrated it into my GEOM PL/SQL package. The procedure’s type signature is as follows:
procedure Split( p_line in mdsys.sdo_geometry,
p_point in mdsys.sdo_geometry,
p_tolerance in number,
p_out_line1 out mdsys.sdo_geometry,
p_out_line2 out mdsys.sdo_geometry )
(The source code for this function is at the end of this article.)
This procedure can be called as follows:
set serveroutput on size unlimited
declare
v_line mdsys.sdo_geometry := mdsys.sdo_geometry(2006,NULL,NULL, mdsys.sdo_elem_info_array(1,2,1,5,2,1), mdsys.sdo_ordinate_array(0,0,5,5,5,5,10,10));
v_point mdsys.sdo_geometry := mdsys.sdo_geometry(2001,null,mdsys.sdo_point_type(2.5,2.5,NULL),NULL,NULL);
v_oline1 mdsys.sdo_geometry;
v_oline2 mdsys.sdo_geometry;
begin
geom.Split(v_line, v_point,0.005,v_oline1,v_oline2);
dbms_output.put_line('Outline1: '||v_oline1.get_WKT());
dbms_output.put_line('Outline2: '||v_oline2.get_wkt());
end;
/
Outline1: LINESTRING (0.0 0.0, 2.5 2.5)
Outline2: MULTILINESTRING ((2.5 2.5, 5.0 5.0), (5.0 5.0, 10.0 10.0))
PL/SQL procedure successfully completed.
Now the Split procedure will handle single (2002) and multi-part (2006) linestring geometries. It will also handle situations where the point:
- Falls on a vertex
- Falls on the line string between two vertices
- Does not fall on the linestring at all (A position on the line at the nearest projection of the point to the line is computed and the linestring is split at this point. The position is computed by simple distance ratios and so will not be accurate for geodetic data.)
Function not Procedure
For those who want a function instead of a procedure so that the function can be used in SQL, I have added a Split function that returns a set of sdo_geometries as follows.
Function Split( p_line in mdsys.sdo_geometry,
p_point in mdsys.sdo_geometry,
p_tolerance in number )
Return codesys.GeometrySetType pipelined;
This function can be used with a table function in a SQL SELECT statement as follows:
select b.*
from table( codesys.geom.split( mdsys.sdo_geometry(2006,NULL,NULL,
mdsys.sdo_elem_info_array(1,2,1,5,2,1),
mdsys.sdo_ordinate_array(0,0,5,5,5,5,10,10)),
mdsys.sdo_geometry(2001,null,mdsys.sdo_point_type(2.5,2.5,NULL),NULL,NULL),
0.005) ) b
/
GEOMETRY(SDO_GTYPE, SDO_SRID, SDO_POINT(X, Y, Z), SDO_ELEM_INFO, SDO_ORDINATES)
-----------------------------------------------------------------------------------------------------------------------------------
SDO_GEOMETRY(2002, NULL, NULL, SDO_ELEM_INFO_ARRAY(1, 2, 1), SDO_ORDINATE_ARRAY(0, 0, 2.5, 2.5))
SDO_GEOMETRY(2006, NULL, NULL, SDO_ELEM_INFO_ARRAY(1, 2, 1, 5, 2, 1), SDO_ORDINATE_ARRAY(2.5, 2.5, 5, 5, 5, 5, 10,10))
Two caveats:
- Split does not handle compound linestrings with 3 point arcs. Any results from calling Split with such an arc will not be correct as no check is made as to arc existence.
- Results for geodetic data will be correct only where the point actually falls on the line: all other situations will result in a rough approximation.
I will describe the simple algorithm underlying this implementation in another posting.
Additionally, in my T_GEOMETRY Object PL/SQL code I have implemented a fuller version of split that handles circular strings and is better with geodetic calculations.
Finally, if there are any errors in the functions presented in this post please let me know and I will endeavour to correct them.
I hope the function is useful.
Here is the implementation
Procedure Split( p_line in mdsys.sdo_geometry,
p_point in mdsys.sdo_geometry,
p_tolerance in number,
p_out_line1 out nocopy mdsys.sdo_geometry,
p_out_line2 out nocopy mdsys.sdo_geometry )
As
cursor c_vectors(p_geometry in mdsys.sdo_geometry)
Is
select rownum as id,
b.startcoord.x as x1,
b.startcoord.y as y1,
b.endcoord.x as x2,
b.endcoord.y as y2
from table(geom.ST_Vectorize(p_geometry)) b;
v_gtype number;
v_part number;
v_element number;
v_num_elements number;
v_vector_id number;
v_start_distance number;
v_x1 number;
v_y1 number;
v_end_distance number;
v_x2 number;
v_y2 number;
v_line_distance number;
v_ratio number;
v_geometry mdsys.sdo_geometry;
v_extract_geom mdsys.sdo_geometry;
v_geom_part mdsys.sdo_geometry;
v_min_dist number;
v_dist number;
v_vector_1 dgs_retea.T_Vector := dgs_retea.T_Vector(null,null,null);
v_vector_2 dgs_retea.T_Vector := dgs_retea.T_Vector(null,null,null);
begin
-- Check inputs
If ( p_point is NULL or p_line is NULL ) Then
RETURN;
End If;
v_gtype := MOD(p_line.Sdo_GType,10);
If ( v_gtype not in (2,6) ) Then
RETURN;
End If;
v_gtype := MOD(p_point.Sdo_GType,10);
If ( MOD(p_point.Sdo_Gtype,10) <> 1 ) Then
RETURN;
End If;
if ( p_tolerance is null ) Then
RETURN;
End If;
-- Get nearest element in multilinestring to process
-- Check number of elements in input line
--
v_num_elements := mdsys.sdo_util.GetNumElem(p_line);
If ( v_num_elements = 1 ) Then
v_geometry := p_line;
v_part := 1;
Else
v_min_dist := 999999999999.99999999; -- All distances should be less than this
<<for_all_vertices>>
FOR v_element IN 1..v_num_elements LOOP
v_extract_geom := mdsys.sdo_util.Extract(p_line,v_element); -- Extract element with all sub-elements
v_dist := mdsys.sdo_geom.sdo_distance(v_extract_geom,p_point,p_tolerance);
If ( v_dist < v_min_dist ) Then
-- dbms_output.put_line('Assigning element || ' || v_element || ' to v_geometry');
v_geometry := v_extract_geom;
v_min_dist := v_dist;
v_part := v_element;
End If;
END LOOP for_all_elements;
End If;
-- We have the line geometry for splitting in v_geometry
-- Find the vector in v_geometry that will be split
--
select id,startdist,x1,y1,enddist,x2,y2,linedist,startdist/(startdist+enddist) as ratio
into v_vector_id,
v_start_distance,
v_x1,
v_y1,
v_end_distance,
v_x2,
v_y2,
v_line_distance,
v_ratio
from (select rownum as id,
b.startcoord.x as x1,
b.startcoord.y as y1,
b.endcoord.x as x2,
b.endcoord.y as y2,
mdsys.sdo_geom.sdo_distance(
mdsys.sdo_geometry(2001,p_point.sdo_srid,
mdsys.sdo_point_type(b.startcoord.x,b.startcoord.y,NULL),NULL,NULL),
p_point,p_tolerance) as startDist,
mdsys.sdo_geom.sdo_distance(
mdsys.sdo_geometry(2001,p_point.sdo_srid,
mdsys.sdo_point_type(b.endcoord.x,b.endcoord.y,NULL),NULL,NULL),
p_point,p_tolerance) as endDist,
mdsys.sdo_geom.sdo_distance(
mdsys.sdo_geometry(2002,p_line.sdo_srid,NULL,
mdsys.sdo_elem_info_array(1,2,1),
mdsys.sdo_ordinate_array(b.startcoord.x,b.startcoord.y,b.endcoord.x,b.endcoord.y)),
p_point,p_tolerance) as linedist
from table(geom.ST_Vectorize(v_geometry)) b
order by 8
)
where rownum < 2;
-- Now do the splitting
If ( v_line_distance = 0 ) Then
-- dbms_output.put('provided point is on the line.');
if ( v_start_distance = 0 ) then
-- dbms_output.put('point can only split line at the start of the vector');
v_vector_1 := NULL;
v_vector_2 := dgs_retea.T_Vector(2,mdsys.sdo_point_type(v_x1,v_y1,null),mdsys.sdo_point_type(v_x2,v_y2,null));
elsif ( v_end_distance = 0 ) then
-- dbms_output.put('point can only split line at the start of the vector');
v_vector_1 := dgs_retea.T_Vector(1,mdsys.sdo_point_type(v_x1,v_y1,null),mdsys.sdo_point_type(v_x2,v_y2,null));
v_vector_2 := NULL;
else
-- dbms_output.put('point is between start and end of vector');
v_vector_1 := dgs_retea.T_Vector(1,mdsys.sdo_point_type(v_x1,v_y1,null),mdsys.sdo_point_type(p_point.sdo_point.x,p_point.sdo_point.y,null));
v_vector_2 := dgs_retea.T_Vector(2,mdsys.sdo_point_type(p_point.sdo_point.x,p_point.sdo_point.y,null),mdsys.sdo_point_type(v_x2,v_y2,null));
end if;
else
If ( v_line_distance = v_start_distance ) then
-- dbms_output.put('point can only split line at the start of the vector');
v_vector_1 := NULL;
v_vector_2 := dgs_retea.T_Vector(2,mdsys.sdo_point_type(v_x1,v_y1,null),mdsys.sdo_point_type(v_x2,v_y2,null));
elsIf ( v_line_distance = v_end_distance ) then
-- dbms_output.put('point can only split line at the end of the vector');
v_vector_1 := dgs_retea.T_Vector(1,mdsys.sdo_point_type(v_x1,v_y1,null),mdsys.sdo_point_type(v_x2,v_y2,null));
v_vector_2 := NULL;
else
-- dbms_output.put('point is between first and last vertex so split point is ratio of start/end distances');
v_vector_1 := dgs_retea.T_Vector(1,mdsys.sdo_point_type(v_x1,v_y1,null),mdsys.sdo_point_type(v_x1+(v_x2-v_x1)*v_ratio,v_y1+(v_y2-v_y1)*v_ratio,null));
v_vector_2 := dgs_retea.T_Vector(2,mdsys.sdo_point_type(v_x1+(v_x2-v_x1)*v_ratio,v_y1+(v_y2-v_y1)*v_ratio,null),mdsys.sdo_point_type(v_x2,v_y2,null));
end if;
End If;
-- dbms_output.put_line('Vector1: ('||v_vector_1.startcoord.x||','||v_vector_1.startcoord.y||')('||v_vector_1.endcoord.x||','||v_vector_1.endcoord.y||')');
-- dbms_output.put_line('Vector2: ('||v_vector_2.startcoord.x||','||v_vector_2.startcoord.y||')('||v_vector_2.endcoord.x||','||v_vector_2.endcoord.y||')');
-- Construct the output geometries
-- Add elements in multi-part geometry before split element to first output line
--
FOR v_element IN 1..(v_part-1) LOOP
-- dbms_output.put_line('Adding element || ' || v_element || ' to out line 1');
v_extract_geom := mdsys.sdo_util.Extract(p_line,v_element); -- Extract element with all sub-elements
p_out_line1 := mdsys.sdo_util.Append(p_out_line1,v_extract_geom);
END LOOP;
-- Now add the vertexes of the split geometry (in v_geometry) to the output lines
--
FOR rec IN c_vectors(v_geometry) LOOP
-- dbms_output.put(rec.id || ': ');
If ( rec.id < v_vector_id ) Then
-- dbms_output.put_line('Add to first part');
if ( rec.id = 1 ) Then
-- dbms_output.put_line('Creating output line1 geometry');
v_geom_part := mdsys.sdo_geometry(2002,p_line.sdo_srid,NULL,
mdsys.sdo_elem_info_array(1,2,1),
mdsys.sdo_ordinate_array(rec.x1,rec.y1,rec.x2,rec.y2));
Else
-- dbms_output.put_line('Append vector to output line1 geometry');
v_geom_part := mdsys.sdo_util.Concat_Lines(v_geom_part,
mdsys.sdo_geometry(2002,p_line.sdo_srid,NULL,
mdsys.sdo_elem_info_array(1,2,1),
mdsys.sdo_ordinate_array(rec.x1,rec.y1,rec.x2,rec.y2)));
End If;
ElsIf ( rec.id = v_vector_id ) Then
If ( v_vector_1 is not NULL ) Then
if ( v_geom_part is null ) Then
-- dbms_output.put_line('Creating output line1 geometry');
v_geom_part := mdsys.sdo_geometry(2002,p_line.sdo_srid,NULL,
mdsys.sdo_elem_info_array(1,2,1),
mdsys.sdo_ordinate_array(v_vector_1.startcoord.x,v_vector_1.startcoord.y,
v_vector_1.endcoord.x, v_vector_1.endcoord.y));
Else
-- dbms_output.put_line('Appending v_vector_1 to first part');
v_geom_part := mdsys.sdo_util.Concat_Lines(v_geom_part,
mdsys.sdo_geometry(2002,p_line.sdo_srid,NULL,
mdsys.sdo_elem_info_array(1,2,1),
mdsys.sdo_ordinate_array(v_vector_1.startcoord.x,v_vector_1.startcoord.y,
v_vector_1.endcoord.x, v_vector_1.endcoord.y)));
End If;
End If;
p_out_line1 := MDSYS.SDO_UTIL.Append(p_out_line1,v_geom_part);
p_out_line2 := NULL;
If ( v_vector_2 is not NULL ) Then
-- dbms_output.put_line(' Add v_vector_2 to p_out_line2 ready to collect up remaining vectors and elements into output line 2');
p_out_line2 := mdsys.sdo_geometry(2002,p_line.sdo_srid,NULL,
mdsys.sdo_elem_info_array(1,2,1),
mdsys.sdo_ordinate_array(v_vector_2.startcoord.x,v_vector_2.startcoord.y,
v_vector_2.endcoord.x,v_vector_2.endcoord.y));
End If;
Else
-- dbms_output.put_line(' Add any remaining vectors to v_geom_part');
if ( p_out_line2 is null ) Then
p_out_line2 := mdsys.sdo_geometry(2002,p_line.sdo_srid,NULL,
mdsys.sdo_elem_info_array(1,2,1),
mdsys.sdo_ordinate_array(rec.x1,rec.y1,rec.x2,rec.y2));
Else
p_out_line2 := mdsys.sdo_util.Concat_Lines(p_out_line2,
mdsys.sdo_geometry(2002,p_line.sdo_srid,NULL,
mdsys.sdo_elem_info_array(1,2,1),
mdsys.sdo_ordinate_array(rec.x1,rec.y1,rec.x2,rec.y2)));
End If;
End If;
END LOOP;
-- Now append any remaining elements in p_line to p_out_line2
FOR v_element IN (v_part+1)..v_num_elements LOOP
-- dbms_output.put_line('Adding element || ' || v_element || ' to out line 2');
v_extract_geom := mdsys.sdo_util.Extract(p_line,v_element); -- Extract element with all sub-elements
p_out_line2 := mdsys.sdo_util.Append(p_out_line2,v_extract_geom);
END LOOP;
End Split;