Retrieve the start and end points of a polyline in Oracle

Strangely enough no such function exists (or at least I couldn’t find it).. PostGIS has ST_StartPoint and ST_EndPoint, SQL Server 2008 Spatial has STStartPoint and STEndPoint but not Oracle. It does however has SDO_LRS.GEOM_SEGMENT_START_PT and SDO_LRS.GEOM_SEGMENT_END_PT but this applies (it seems) only for measured polylines – and you need to have the full Oracle Spatial installed. With versions that only have Locator you are out of luck.

So putting my thrifty hat again, I put together 2 simple functions (with the very imaginative names get_start_point and get_end_point) that do just that and included them in the package. Functions should work with both 2D and 3D data. You can download the update “thrifty” package here (or from the box widget on the right).

SQL Server 2008 Spatial and Oracle Spatial comparison and cheat sheet

This is a brief comparison between SQL Server 2008 and Oracle’s spatial functions along with relevant links and simple examples. It is by no means exhaustive but it should help people when migrating.

Note that all Oracle examples and links here refer to the 10g R2 release.

UPDATED: Added link for updated version of SQLSpatialTools from Geographica.

Spatial Data

Feature Oracle 10g SQL Server 2008
Spatial Data types

Defined by implementing the SDO_GEOMETRY type:

CREATE TYPE sdo_geometry AS OBJECT (
SDO_GTYPE NUMBER,
SDO_SRID NUMBER,
SDO_POINT SDO_POINT_TYPE,
SDO_ELEM_INFO SDO_ELEM_INFO_ARRAY,
SDO_ORDINATES SDO_ORDINATE_ARRAY);

Defines two types of spatial data depending whether they contain geographic (round-earth coordinate system) or geometric (Euclidean/flat coordinate system) data: Geometry and Geography. Both types are implemented as a .NET common language runtime (CLR) data types.

Spatial Tables Statement to create a spatial table:
CREATE TABLE sp_table(sp_id NUMBER PRIMARY KEY,
  shape SDO_GEOMETRY);

Statement to create a spatial table (geometry):

CREATE TABLE sp_table 
(sp_id int IDENTITY (1,1), shape geometry);

Statement to create a spatial table (geography):

CREATE TABLE sp_table 
(sp_id int IDENTITY (1,1),shape geography);

Geometry Metadata For spatial tables to perform correctly, Oracle requires the update of the geometry metadata views. These views describe the dimensions, lower and upper bounds, and tolerance in each dimension for each spatial table and are stored in a global table owned by the MDSYS user.

The main view is USER_SDO_GEOM_METADATA which contains metadata information for all spatial tables owned by the user (schema).

Example statement to insert data into the USER_SDO_GEOM_METADATA view:

INSERT INTO user_sdo_geom_metadata
(TABLE_NAME,COLUMN_NAME,
  DIMINFO, SRID)
VALUES (‘SP_TABLE’, ‘SHAPE’, 
SDO_DIM_ARRAY(    SDO_DIM_ELEMENT(‘X’, 0, 20, 0.005), 
SDO_DIM_ELEMENT(‘Y’, 0, 20, 0.005)), NULL);

SQL Server does not have an equivalent of Oracle’s geometry metadata views
Spatial Indexes Supports R-tree and Quadtree indexing. However, Oracle states that:

the use of quadtree indexes is discouraged, and you are strongly encouraged to use R-tree indexing

For more information about Oracle’s spatial indexing click here.

Creating a spatial (R-tree) index:

CREATE INDEX sp_table_spidx   ON sp_table(shape) INDEXTYPE IS MDSYS.SPATIAL_INDEX;

Uses B-tree indexing. For more information on how spatial indexing is handled, click here.

Creating spatial indexing syntax varies depending whether the column is of type Geometry or Geography.

Creating a spatial index on a Geography column:

CREATE SPATIAL INDEX sp_table_spidx    ON sp_table(shape);

Creating a spatial index on a Geometry column:

CREATE SPATIAL INDEX sp_table_spidx ON sp_table(shape) WITH (
BOUNDING_BOX = ( 0, 0, 500, 200 ), 
GRIDS = ( LEVEL_4 = HIGH, LEVEL_3 = MEDIUM ) );

More examples can be found here

OGG Compliance Conforms to the Open Geospatial Consortium (OGC) Simple Features for SQL Specification version 1.1.0. Conforms to the Open Geospatial Consortium (OGC) Simple Features for SQL Specification version 1.1.0

Spatial object

Constructors

Examples
More geometry examples here
Examples
       Points

(Notice that 4326 is the SRID)

1. Using the SDO_GEOMETRY type:

insert into sdo_points (sp_id, shape)
values (1, SDO_GEOMETRY(
      2001,
      4326,
      SDO_POINT_TYPE(40, 22, NULL),  NULL,   NULL));

2. Using the Well-Known text (WKT) syntax:

SELECT SDO_GEOMETRY(‘POINT(-79 37)’) FROM DUAL;

(Notice that 4326 is the SRID)

1. Using the Well-Known text (WKT) syntax (STPointFromText function):

INSERT INTO GeomTable (geom)

VALUES (geography::STPointFromText(‘POINT(40 22)’, 4326))

2. Using the geometry::Point or geography::Point syntax:

INSERT INTO GeogTable (geog) VALUES (geography::Point(40, 22, 4326))

3. Using the GML syntax (GeomFromGml function)

INSERT INTO GeomTable (geom)
     VALUES (geography::GeomFromGml(
‘<Point xmlns="
http://www.opengis.net/gml">
<pos>40 22</pos></Point>’
, 4326))

       Simple Lines

Using the SDO_GEOMETRY type:

INSERT INTO sdo_lines VALUES(
  1, SDO_GEOMETRY(
    2002,
    NULL, — NULL SRID
    NULL,
    SDO_ELEM_INFO_ARRAY(1,2,1),
    SDO_ORDINATE_ARRAY(411392.088118, 4493608.698554, 411436.431582, 4493649.698236)));

(You can also create line geometries using the WKT syntax as above)

Using the Well-Known text (WKT) syntax (STLineFromText function):

INSERT INTO GeogTable (geog)
VALUES (geography::STLineFromText(‘LINESTRING(-122.360 47.656,
-122.343 47.656 )’, 4326));

(You can also create line geometries using the WKT and GML syntax as above)

       Simple Polygons

Using the SDO_GEOMETRY type:

INSERT INTO sdo_polys VALUES(
2,
MDSYS.SDO_GEOMETRY(
2003, — 2-dimensional polygon
NULL,
NULL,
MDSYS.SDO_ELEM_INFO_ARRAY(1,1003,1), — one polygon (exterior polygon ring)
MDSYS.SDO_ORDINATE_ARRAY(5,1, 8,1, 8,6, 5,7, 5,1)));

(You can also create polygon geometries using the WKT syntax as above)

Using the Well-Known text (WKT) syntax (STPolyFromText function):

INSERT INTO GeogTable (geog) VALUES (geography::STPolyFromText(‘POLYGON((-122.358 47.653,
-122.348 47.649, -122.348 47.658, -122.358 47.658, -122.358 47.653))’
, 4326))

(You can also create polygon geometries using the WKT and GML syntax as above)

3D/4D support

Both DBs support 3D and 4D (X,Y,Z,M) shapes

Loader and conversion  Tools Oracle makes available for download from its website two utilities to convert shapefile to oracle spatial tables. The first one is a command line utility shp2sdo.exe,and the second a Java shapefile converter. You can find more information about these utilities here. There is a free utility called Shape2Sql available from the SharpGIS website

Spatial Functions

The table list some commonly used spatial functions.

Function Description Oracle 10g

Spatial Functions see also the
SDO_GEOM package
SQL Server 2008

All Geography functions
All Geometry functions

Retrieve the spatial extent of a given table(Minimum Bounding Rectangle (MBR)) Use the SDO_AGGR_MBR function, e.g.:

SELECT SDO_AGGR_MBR(shape) FROM cola_markets;

Does not contain a similar function in the standard install. You will have to use the GeographyUnionAggregate and GeometryEnvelopeAggregate functions available in the SQL Server Spatial Tools CRL package available from Codeplex.
e.g.

select state, dbo.
GeometryEnvelopeAggregate(shape_geom) from zipcodes

Identify objects within a distance

Use the SDO_WITHIN_DISTANCE function e.g.:

SELECT c.name FROM cola_markets c WHERE SDO_WITHIN_DISTANCE(c.shape,
  SDO_GEOMETRY(2003, NULL, NULL, SDO_ELEM_INFO_ARRAY(1,1003,3),
    SDO_ORDINATE_ARRAY(4,6, 8,8)),
  ‘distance=10’) = ‘TRUE’;

Use the STDistance (Geography) or STDistance (Geometry) e.g.:

DECLARE @g geometry;
DECLARE @h geometry;
SET @g = geometry::STGeomFromText(‘POLYGON((0 0, 2 0, 2 2, 0 2, 0 0))’, 0);
SET @h = geometry::STGeomFromText(‘POINT(10 10)’, 0);
SELECT @g.STDistance(@h);

Create a buffer around a spatial object Use the SDO_GEOM.SDO_BUFFER function e.g.:

SELECT c.name, SDO_GEOM.SDO_BUFFER(c.shape, m.diminfo, 1)
  FROM cola_markets c, user_sdo_geom_metadata m
  WHERE m.table_name = ‘COLA_MARKETS’  AND m.column_name = ‘SHAPE’
  AND c.name = ‘cola_a’;

Use the STBuffer (Geography) or STBuffer (geometry) e.g.:

DECLARE @g geometry;
SET @g = geometry::STGeomFromText(‘LINESTRING(0 0, 4 0)’, 0);
SELECT @g.STBuffer(1).ToString();

Get the length of an object Use the SDO_GEOM.SDO_LENGTH function, e.g.:

SELECT c.name, SDO_GEOM.SDO_LENGTH(c.shape, m.diminfo)
  FROM cola_markets c, user_sdo_geom_metadata m
  WHERE m.table_name = ‘COLA_MARKETS’ AND m.column_name = ‘SHAPE’;

Use the STLength (Geography) or STLength (Geometry) e.g.:

 

DECLARE @g geography;
SET @g = geography::STGeomFromText(‘LINESTRING(-122.360 47.656, -122.343 47.656)’, 4326);
SELECT @g.STLength();

Get the area of a polygon   Use the STArea (Geography) or STArea (Geometry) e.g.:

DECLARE @g geography;
SET @g = geography::STGeomFromText(‘POLYGON((-122.358 47.653, -122.348 47.649, -122.348 47.658, -122.358 47.658, -122.358 47.653))’, 4326);
SELECT @g.STArea();

Check the validity of an object Use the SDO_GEOM. VALIDATE_GEOMETRY_WITH_CONTEXT function e.g.:

SELECT c.name, SDO_GEOM. VALIDATE_GEOMETRY_WITH_CONTEXT(c.shape, 0.005)
   FROM cola_markets c WHERE c.name = ‘cola_invalid_geom’;

Use the STIsValid (Geometry) e.g.:

DECLARE @g geometry;
SET @g = geometry::STGeomFromText(‘LINESTRING(0 0, 2 2, 1 0)’, 0);
SELECT @g.STIsValid();

The Geography data type does not –for reasons unknown to me- include a STIsValid function.

Linear Referencing Contains a large number of functions for linear referencing. Lots of examples can be found here SQL Server does not support any true Linear referencing functions. SQL Server Spatial Tools however, provides some rudimentary linear referencing functions such as LocateAlongGeog and LocateAlongGeom which will return the point at a given instance along a linear object.
UPDATED: Geographica provides an updated version of SQLSpatialTools which includes a new function to display a linear event as well (a “LocateLineAlongGeom” function but called CreateLinearReferenceFeature. More details can be found at this address: http://geographika.co.uk/linear-referencing-in-sql-server-2008

Create point geometries dynamically using dynamic segmentation in SQLServer 2008

UPDATED: Change trigger to allow for multiple linear geometries

I have started experimenting with the new spatial functionality in SQL Server 2008. Having used Oracle Spatial for a quite a few years this is my first attempt with SQL Server so my apologies in advance if there is an easier way of doing this! 🙂

So the scenario is as follows:- Consider you have a layer containing A SINGLE multiple linear geometries with the following structure called ROUTES:

image

The structure can be different- for now all we care about is the [Shape] and [ID] fields.

We also have another table called EVENTS- this is the table we will be updating using the trigger. The table structure is:

image

Where ROUTE_ID is the ID the Route the event is  located along, (as a FK to ID field in the ROUTES table) and  START_MP is the distance from the start of the Route.

The task at hand was to create a trigger so any changes in the START_MP value would update the shape column in the EVENTS table, creating a Point geometry at that location.

That was easy enough using the SQL Server Spatial tools available from http://sqlspatialtools.codeplex.com. After registering the SQLSpatialTools.dll in my database I was able to use the LocateAlongGeom function in a trigger which would deliver the point geometries. The trigger source is:

ALTER TRIGGER [dbo].[upd_shape_trg]
   ON  [dbo].[EVENTS]
   AFTER UPDATE
AS

BEGIN 
  DECLARE @lineshape geometry
    DECLARE @routeid int
    --SELECT @lineshape = r.shape from marathonia_diadromi_spatial r
    --, events e
    --where r.id= e.route_id
    SELECT @routeid=route_id from deleted
    SELECT @lineshape = shape from marathonia_diadromi_spatial where id=@routeid
    -- SET NOCOUNT ON added to prevent extra result sets from
    -- interfering with SELECT statements.
    SET NOCOUNT ON;
   UPDATE [mpoint].[dbo].[EVENTS]
   SET [SHAPE] = [mpoint].[dbo].[LocateAlongGeom](@lineshape, start_mp).ToString()

END

Note that right now this would only work if the ROUTES table contains a single linear geometry. As said- first attempt only this! More later…

Oracle Spatial for the thrifty developer- Part 3: The basics of dynamic segmentation…

In this post I will show you a function that will return a linear segment along a line between a start and end length. This function together with the PointAlong() function as discussed in the previous post could be used as the basic blocks for building dynamic segmentation or address geocoding routines. At the end of the day these two processes share the same basic principle: to create a geometry (point or line) for a feature (known as an “event”), based on descriptive information only such as start/end address or start/end linear reference plus the underlying linear geometry the feature is located against.

The code for the LineAlong function is shown below:


FUNCTION lineAlong (pi_geom IN mdsys.sdo_geometry
,pi_from IN NUMBER
,pi_to IN NUMBER
,pi_srid IN NUMBER DEFAULT NULL)
RETURN mdsys.sdo_geometry IS
--
last1 INTEGER;
st_set BOOLEAN := TRUE;
i INTEGER;
i_v INTEGER;
n_srid NUMBER;
l_segment mdsys.sdo_geometry;
x_current NUMBER;
y_current NUMBER;
x_next NUMBER;
y_next NUMBER;
l_geom_type NUMBER;
n_total_length NUMBER;
n_seg_length NUMBER;
n_addr_length1 NUMBER;
n_addr_length2 NUMBER;
n_run_length NUMBER;
n_dist NUMBER;
n_grad NUMBER;
found_start_pnt BOOLEAN;
n_dist_to_point NUMBER;
l_ords1 mdsys.sdo_ordinate_array; -- Holds the xy of the from address point
l_ords2 mdsys.sdo_ordinate_array; -- Holds the xy of the to address point
l_ords_all mdsys.sdo_ordinate_array; -- Holds the xy of the address range
n_ords_count NUMBER;
inval_geom_type EXCEPTION;
--
BEGIN
-- Do some validation
IF pi_from > pi_to THEN
raise_application_error(-20000,'[From] value must be greater than the [To] value');
END IF;
n_total_length := get_line_length(pi_geom); -- Get total length of shape
IF pi_to > n_total_length THEN
raise_application_error(-20000,'[To] value cannot be greater than the total line length');
END IF;
--
-- Initialise vars
n_grad := NULL;
n_dist := NULL;
l_geom_type := pi_geom.sdo_gtype;
n_srid := pi_srid;
found_start_pnt := FALSE;
--
--Check geometry type. Must be a 2d line
IF l_geom_type = 2002 THEN
-- Get number of vertices
last1 := pi_geom.sdo_ordinates.LAST;
--dbms_output.put_line('Total vertices:' || to_char(last1 / 2));
-- Initialise counter
i := 1;
-- Calculate absolute distances of from address point from segment start
n_addr_length1 := pi_from;
n_addr_length2 := pi_to;
--
-- Start looping through vertices
WHILE TRUE LOOP
IF st_set THEN
i_v := 1;
x_current := pi_geom.sdo_ordinates(i);
y_current := pi_geom.sdo_ordinates(i + 1);
x_next := pi_geom.sdo_ordinates(i + 2);
y_next := pi_geom.sdo_ordinates(i + 3);
-- Greate segment geometry
l_segment := get_line_from_xy(x_current
,y_current
,x_next
,y_next
,n_srid);
-- Get segment length
n_seg_length := get_line_length(l_segment);
n_run_length := n_seg_length;
--Find xy for start point address
IF found_start_pnt = FALSE THEN
IF n_run_length >= n_addr_length1 THEN
n_grad := get_grad_from_xy(x_current
,y_current
,x_next
,y_next);
n_dist := get_length_from_xy(x_current
,y_current
,x_next
,y_next);
n_dist_to_point := n_addr_length1;
l_ords1 := get_point_at_distance_as_ords(x_current
,y_current
,n_dist_to_point
,n_grad
,n_srid);
-- Insert into ordinate array
l_ords_all := mdsys.sdo_ordinate_array(l_ords1(1),
l_ords1(2));
found_start_pnt := TRUE;
END IF;
END IF;
--Find xy for end point address
IF n_run_length >= n_addr_length2 THEN
n_grad := get_grad_from_xy(x_current,
y_current,
x_next,
y_next);
n_dist := get_length_from_xy(x_current,
y_current,
x_next,
y_next);
n_dist_to_point := n_addr_length2;
l_ords2 := get_point_at_distance_as_ords(x_current
,y_current
,n_dist_to_point
,n_grad);
--
/*dbms_output.put_line('End X: ' || to_char(l_ords2(1)) ||
' End Y: ' || to_char(l_ords2(2)));
*/
IF found_start_pnt = FALSE THEN
raise_application_error(-20000,'Invalid start length');
END IF;
-- Insert into ordinate array
n_ords_count := l_ords_all.COUNT;
l_ords_all.EXTEND(2);
l_ords_all(n_ords_count + 1) := l_ords2(1);
l_ords_all(n_ords_count + 2) := l_ords2(2);
EXIT;
END IF;
--dbms_output.put_line('segment counter:' || to_char(i_v));
st_set := FALSE;
i := i + 2;
i_v := i_v + 1; --Increment segment counter
ELSE
-- Remaining segments
IF i + 1 >= last1 THEN
-- We are on the last pair
IF n_grad IS NULL OR n_dist IS NULL THEN
raise_application_error(-20000,'Could not create line');
END IF;
EXIT;
END IF;
--
x_current := pi_geom.sdo_ordinates(i);
y_current := pi_geom.sdo_ordinates(i + 1);
x_next := pi_geom.sdo_ordinates(i + 2);
y_next := pi_geom.sdo_ordinates(i + 3);
--
IF found_start_pnt = TRUE THEN
-- Insert into ordinate array
n_ords_count := l_ords_all.COUNT;
l_ords_all.EXTEND(2);
l_ords_all(n_ords_count + 1) := x_current;
l_ords_all(n_ords_count + 2) := y_current;
--l_ords_all(n_ords_count+3) := x_next;
--l_ords_all(n_ords_count+4) := y_next;
END IF;
-- Greate segment geometry
l_segment := get_line_from_xy(x_current
,y_current
,x_next
,y_next
,n_srid);
-- Get segment length
n_seg_length := get_line_length(l_segment);
n_run_length := n_run_length + n_seg_length;
--Find xy for start point address
IF found_start_pnt = FALSE THEN
IF n_run_length >= n_addr_length1 THEN
n_grad := get_grad_from_xy(x_current,
y_current,
x_next,
y_next);
n_dist := get_length_from_xy(x_current,
y_current,
x_next,
y_next);
n_dist_to_point := n_seg_length -
(n_run_length - n_addr_length1);
l_ords1 := get_point_at_distance_as_ords(x_current
,y_current
,n_dist_to_point
,n_grad);
-- Insert into ordinate array
l_ords_all := mdsys.sdo_ordinate_array(l_ords1(1),
l_ords1(2));
found_start_pnt := TRUE;
END IF;
END IF;
i := i + 2;
--dbms_output.put_line('segment counter:' || to_char(i_v));
END IF;
--
--Find xy for end point address
IF n_run_length >= n_addr_length2 THEN
n_grad := get_grad_from_xy(x_current,
y_current,
x_next,
y_next);
n_dist := get_length_from_xy(x_current,
y_current,
x_next,
y_next);
n_dist_to_point := n_seg_length - (n_run_length - n_addr_length2);
l_ords2 := get_point_at_distance_as_ords(x_current,
y_current,
n_dist_to_point,
n_grad);
IF found_start_pnt = FALSE THEN
raise_application_error(-20000,'Invalid start length');
END IF;
-- Insert into ordinate array
n_ords_count := l_ords_all.COUNT;
l_ords_all.EXTEND(2);
l_ords_all(n_ords_count + 1) := l_ords2(1);
l_ords_all(n_ords_count + 2) := l_ords2(2);
END IF;
END LOOP;
/* -- DEBUG INFO
FOR i IN 1 .. l_ords_all.COUNT LOOP
dbms_output.put_line(to_char(l_ords_all(i)));
END LOOP;
*/
RETURN mdsys.sdo_geometry(2002,
n_srid,
NULL,
mdsys.sdo_elem_info_array(1,2,1),l_ords_all);
END IF;
EXCEPTION
WHEN inval_geom_type THEN
raise_application_error(-20003, 'Function does not support this spatial type (' || to_char(l_geom_type) || ')');
END lineAlong;

You can download the whole package with all the functions here.

Oracle Spatial for the thrifty developer- Part 2: Calculating line lengths and more…

In this post we will talk about calulating line lengths as well as finding a point at a percentage along a line. This second function will prove quite useful as will be seen in next posts to implement simple address geocoding or building your own version of dynamic segmentation routines.
You can download the whole package here.

Calculate Line Lengths
This will return the length of a single-part polyline

FUNCTION get_line_length(pi_geom IN mdsys.sdo_geometry) RETURN NUMBER IS
last1 INTEGER;
st_set BOOLEAN := TRUE;
i INTEGER;
x_current NUMBER;
y_current NUMBER;
x_next NUMBER;
y_next NUMBER;
l_geom_type NUMBER;
n_total_length NUMBER;
ntol NUMBER;
inval_geom_type EXCEPTION;
BEGIN
--Initialize vars
l_geom_type := pi_geom.sdo_gtype;
n_total_length := 0;
i := 1;
last1 := pi_geom.sdo_ordinates.LAST;
ntol:= c_tol;
-- Validate input params
IF l_geom_type = 2002 THEN
-- Start looping through vertices
WHILE TRUE LOOP
IF st_set THEN
x_current := pi_geom.sdo_ordinates(i);
y_current := pi_geom.sdo_ordinates(i + 1);
x_next := pi_geom.sdo_ordinates(i + 2);
y_next := pi_geom.sdo_ordinates(i + 3);
n_total_length := n_total_length +
sdo_geom.sdo_distance(
get_pt_from_xy(x_current,
y_current),
get_pt_from_xy(x_next,
y_next),
to_number(ntol));
st_set := FALSE;
i := i + 2;
ELSE
IF i + 1 >= last1 THEN
-- We are on the last pair
RETURN n_total_length;
END IF;
--
x_current := pi_geom.sdo_ordinates(i);
y_current := pi_geom.sdo_ordinates(i + 1);
x_next := pi_geom.sdo_ordinates(i + 2);
y_next := pi_geom.sdo_ordinates(i + 3);
i := i + 2;
n_total_length := n_total_length +
sdo_geom.sdo_distance(
get_pt_from_xy(
x_current,
y_current),
get_pt_from_xy(
x_next,
y_next),
to_number(ntol));
END IF;
END LOOP;
RETURN n_total_length;
ELSE
RETURN NULL;
END IF;
EXCEPTION
WHEN inval_geom_type THEN
raise_application_error(-20000,
'Funtion does not support this spatial type
(' || to_char(l_geom_type) || ')');
END get_line_length;

Get Point Along Line

This function will return a point at a percentage along a Polyline:

FUNCTION pointAlong(pi_line IN mdsys.sdo_geometry
,pi_pct IN NUMBER
,pi_srid IN NUMBER)
RETURN mdsys.sdo_geometry IS
---
last1 INTEGER;
st_set BOOLEAN := TRUE;
i INTEGER;
l_segment mdsys.sdo_geometry;
x_current NUMBER;
y_current NUMBER;
x_next NUMBER;
y_next NUMBER;
l_geom_type NUMBER;
n_total_length NUMBER;
n_seg_length NUMBER;
n_addr_length NUMBER;
n_run_length NUMBER;
n_dist NUMBER;
n_grad NUMBER;
l_pnt mdsys.sdo_geometry;
v_dist_to_point NUMBER;
inval_geom_type EXCEPTION;
inval_pct EXCEPTION;
--
BEGIN
-- Initialise vars
n_grad := NULL;
n_dist := NULL;
l_geom_type := pi_line.sdo_gtype;
--
IF pi_pct NOT BETWEEN 0 AND 100 THEN
RAISE inval_pct;
END IF;
--
--Check that we have a simple line geometry
--
IF l_geom_type = 2002 THEN
-- Get number of vertices
last1 := pi_line.sdo_ordinates.LAST;
-- Initialise counter
i := 1;
-- Get total length of line
n_total_length := get_line_length(pi_line);
-- Calculate required length
n_addr_length := n_total_length*(pi_pct/100);
WHILE TRUE LOOP
IF st_set THEN
x_current := pi_line.sdo_ordinates(i);
y_current := pi_line.sdo_ordinates(i + 1);
x_next := pi_line.sdo_ordinates(i + 2);
y_next := pi_line.sdo_ordinates(i + 3);
-- Greate segment geometry
l_segment := get_line_from_xy(x_current
,y_current
,x_next
,y_next
,pi_srid);
-- Get segment length
n_seg_length := get_line_length(l_segment);
n_run_length := n_seg_length;
IF n_run_length >= n_addr_length THEN
-- Found segment where point should be on. Calculate exact location and exit loop
n_grad:= get_grad_from_xy(x_current,
y_current,
x_next,
y_next);
n_dist:= get_length_from_xy(x_current,
y_current,
x_next,
y_next);
v_dist_to_point := n_addr_length;
EXIT;
END IF;
st_set := FALSE;
i := i + 2;
ELSE
IF i + 1 >= last1 THEN
-- We are on the last pair
IF n_grad IS NULL OR n_dist IS NULL THEN
-- Something went wrong...
raise_application_error(
-20000
,'Could not loop through geometry vertices. Invalid Geometry?');
END IF;
EXIT;
END IF;
x_current := pi_line.sdo_ordinates(i);
y_current := pi_line.sdo_ordinates(i + 1);
x_next := pi_line.sdo_ordinates(i + 2);
y_next := pi_line.sdo_ordinates(i + 3);
-- Greate segment geometry
l_segment := get_line_from_xy(x_current
,y_current
,x_next
,y_next
,pi_srid);
-- Get segment length
n_seg_length := get_line_length(l_segment);
n_run_length := n_run_length + n_seg_length;
IF n_run_length >= n_addr_length THEN
n_grad := get_grad_from_xy(x_current,
y_current,
x_next,
y_next);
n_dist := get_length_from_xy(x_current,
y_current,
x_next,
y_next);
v_dist_to_point := n_seg_length -
(n_run_length - n_addr_length);
EXIT;
END IF;
i := i + 2;
END IF;
END LOOP;
l_pnt := get_point_at_distance(x_current
,y_current
,v_dist_to_point
,n_grad);
RETURN l_pnt;
ELSE
RAISE inval_geom_type;
END IF;
EXCEPTION
WHEN inval_geom_type THEN
raise_application_error(-20000,'Funtion does not support this spatial type (
' || to_char(l_geom_type) || ')');
WHEN inval_pct THEN
raise_application_error(-20000,'Value must be between 0 and 100');
END pointAlong;

Oracle Spatial for the thrifty developer: Part 1

Oracle Spatial contains a wealth of spatial functions that cater for most needs. Problem is of course they all come with a ridiculously high price tag! No matter how simple your application may be, if you want to take full advantage of spatial functionality you have to pay for the Enterprise Edition of Oracle, plus the so-called ‘Spatial Option’ i.e. 1000’s of euros (or dollars,  or whatever). A price that actually got increased in June 2009 as reported here.

Luckily, Oracle also provides Oracle Locator which is available on every Oracle version and is a cut down version of spatial which allows you to basically store spatial data in your tables along with some (actually 3) spatial functions. These are:

  • SDO_GEOM.SDO_DISTANCE()
  • SDO_GEOM.VALIDATE_GEOMETRY_WITH_CONTEXT()
  • SDO_GEOM.VALIDATE_LAYER_WITH_CONTEXT() 

So what I will try to do here is attempt to extent this poor offering, by providing a set of simple spatial functions that can be run on any Oracle edition including Oracle 10g Express. This obviously by no means intends to replace the full- blown oracle functions but should hopefully be useful for people who require some simple spatial functionality for their applications. This ‘series’ of  posts will start with some simple coordinate geometry functions and will end (at least at the time of writing these lines) with creating some somewhat more complex routines like returning a point along a line at a specified instance and line merging.

Create a Point from a pair of XYs

Dead simple as you may expect, the function in question is:

FUNCTION get_pt_from_xy(pi_x IN NUMBER
,pi_y IN NUMBER
,pi_srid NUMBER DEFAULT NULL) RETURN mdsys.sdo_geometry IS
retval mdsys.sdo_geometry;
BEGIN
--
retval := mdsys.sdo_geometry(2001,
pi_srid,
mdsys.sdo_point_type(pi_x,pi_y,NULL),
NULL,
NULL);
RETURN retval;
END;

Get Length from a pair of XYs

This function returns the straight line/cartesian length given the start and end XY coordinates

FUNCTION get_length_from_xy(pi_x1 NUMBER
,pi_y1 NUMBER
,pi_x2 NUMBER
,pi_y2 NUMBER) RETURN NUMBER IS
retval NUMBER;
BEGIN
retval := sqrt(power((pi_x2 - pi_x1),
2) + power((pi_y2 - pi_y1),
2));
RETURN retval;
END get_length_from_xy;

Get the angle from a pair of XYs

This function returns the the angle between 2 points defined by start and end XY coordinates


FUNCTION get_grad_from_xy(pi_x1 NUMBER
,pi_y1 NUMBER
,pi_x2 NUMBER
,pi_y2 NUMBER) RETURN NUMBER IS
retval NUMBER;
dx NUMBER;
dy NUMBER;
gab NUMBER;
BEGIN
dx := pi_x2 - pi_x1;
dy := pi_y2 - pi_y1;
gab := atan(abs(dx / dy));
------Find quadrant ---------
-- IF DX >0
IF dx > 0 THEN
IF dy > 0 THEN
retval := gab;
END IF;
IF dy = 0 THEN
retval := c_pi / 2;
END IF;
IF dy < 0 THEN
retval := c_pi - gab;
END IF;
END IF;
-- IF DX < 0
IF dx < 0 THEN
IF dy > 0 THEN
retval := (c_pi * 2) - gab;
END IF;
IF dy = 0 THEN
retval := (c_pi * 1.5) - gab;
END IF;
IF dy < 0 THEN
retval := c_pi + gab;
END IF;
END IF;
-- IF DX=0
IF dx = 0 THEN
IF dy > 0 THEN
retval := 0;
END IF;
IF dy = 0 THEN
-- Start and end points are the same!
raise_application_error(-20014,'Could not determine the angle');
END IF;
IF dy < 0 THEN
retval := c_pi;
END IF;
END IF;
RETURN retval;
END get_grad_from_xy;

Return distance at Point
This function returns the point at a specified distance and angle from another given point.

FUNCTION get_point_at_distance(pi_x NUMBER
,pi_y NUMBER
,pi_dist NUMBER
,pi_grad NUMBER
,pi_srid NUMBER DEFAULT NULL)
RETURN mdsys.sdo_geometry IS
----------------------------------------------------------------------------------------
n_x NUMBER;
n_y NUMBER;
retval mdsys.sdo_geometry;
BEGIN
n_x := pi_x + pi_dist * sin(pi_grad);
n_y := pi_y + pi_dist * cos(pi_grad);
retval:=mdsys.sdo_geometry(2001
,pi_srid
,mdsys.sdo_point_type(n_x,n_y,NULL)
,NULL
,NULL);
RETURN retval;
 

Next : Oracle Spatial for the thrifty developer- Part 2: Calculating line lengths and more…