Convert Oracle Spatial data to geojson

So at some point I needed to convert Oracle Spatial data (SDO) into geojson format. I quickly came across this excellent post from SpatialDB Advisor. Which basically does exactly that. Converts SDO data into geojson. However -as he already states in his post- the function does not support attributes. So based on this original function, I created a package that will create a FULL geojson including attributes using a bit of dynamic SQL (DBMS_SQL package). And in fact, you can select which attributes by providing a select statement as a parameter. I also *think* I fixed a couple of minor bugs in the original sdo2geojson function.

You can find the package (ora2geojson.rar) here. Unzip the file and run the package header  (.pkh) and body (.pkb) files e.g.

SQL>@ora2geojson.pkh

and

SQL>@ora2geojson.pkb

Example:

SELECT ora2geojson.sdo2geojson('select * from v_segments_all_sdo',ROWID, shape) FROM v_segments_all_sdo;

Where v_segments_all is the spatial table or view and shape is the name of the spatial column. The ROWID is needed to collect the attributes. I am sure there are more sufficient ways to do this but this one seemed to work fine for me.

A sample usage and result is shown below:

sdo2geojson

 

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…