Extensive spatial functions package

From WickyWiki

This is an extensive package with many PL/SQL and Spatial functions.

define SRID=NULL
define allLEFT=3.0
define allBOTTOM=50
define allRIGHT=8.0
define allTOP=54
define accX=0.00000005
define accY=0.00000005

/* system privilages needed for the following settings

connect system

--disable the recyclebin
alter system set recyclebin=off;

--grant for oracle bug workaround
grant update on mdsys.SDO_INDEX_METADATA_TABLE to geouser;

--explicit grants for use from within procedures
grant create any table to geouser;
grant create any sequence to geouser;

connect geouser
*/

-- create insert statements for spatial metadata:
select 'INSERT INTO USER_SDO_GEOM_METADATA VALUES('''|| tname ||''','''|| cname ||''', MDSYS.SDO_DIM_ARRAY(MDSYS.SDO_DIM_ELEMENT(''X'', &allLEFT, &allRIGHT, &accX),MDSYS.SDO_DIM_ELEMENT(''Y'', &allBOTTOM, &allTOP, &accY) ), &SRID );' as sql
  from col
  where col.coltype like '%SDO_GEOMETRY%'
  and not tname like '%$%'
  and cname = 'SHAPE'
  and NOT (col.tname, col.cname) in (
  select md.table_name, md.column_name
    from user_sdo_geom_metadata md)
  order by tname, cname
  ;
  
-- repair, rebuild and create spatial indexes
call geo.geoCheck_sp_idx(NULL, 1);

/* package preperation
-------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------
*/

-- log table
create sequence s_script_log START WITH 1 INCREMENT BY 1;
create sequence s_script_log_step START WITH 1 INCREMENT BY 1;

-- select s_script_log.nextval, s_script_log.nextval from dual;
-- drop table script_log;
create table script_log(id number
  , logtime DATE default sysdate
  , todo NUMBER
  , step NUMBER
  , lvl number(1) default 1
  , msg VARCHAR2(500));
  
CREATE or replace TRIGGER script_log_TRIG_ID BEFORE INSERT ON script_log
  for each row begin select s_script_log.nextval into :new.ID from dual; end;
/

-- script_log tabel opschonen, laatste 100 logregels bewaren
delete from script_log where id<
  (select min(id) from 
  (select id from script_log order by id desc) where rownum<=100)
  ;

-- geo tables
define geotabel=tmp_geoResult
define geokolom=SHAPE
DROP INDEX &geotabel._SP_IDX;
drop table &geotabel.;
DELETE FROM USER_SDO_GEOM_METADATA WHERE TABLE_NAME=Upper('&geotabel');
create table &geotabel.(rowid1 rowid, rowid2 rowid, datasrcid number(22,11), datasrcid2 number(22,11), shape SDO_geometry);
INSERT INTO USER_SDO_GEOM_METADATA VALUES
  (Upper('&geotabel'),Upper('&geokolom'), MDSYS.SDO_DIM_ARRAY
  (MDSYS.SDO_DIM_ELEMENT('X', &allLEFT, &allRIGHT, &accX),
    MDSYS.SDO_DIM_ELEMENT('Y', &allBOTTOM, &allTOP, &accY) ), &SRID );
CREATE INDEX &geotabel._SP_IDX ON &geotabel. ( &geokolom ) INDEXTYPE IS MDSYS.SPATIAL_INDEX;
CREATE INDEX &geotabel._rowid1_IDX ON &geotabel. ( rowid1 );
CREATE INDEX &geotabel._rowid2_IDX ON &geotabel. ( rowid2 );
CREATE INDEX &geotabel._datasrcid_IDX ON &geotabel. ( datasrcid );
CREATE INDEX &geotabel._datasrcid2_IDX ON &geotabel. ( datasrcid2 );

define geotabel=tmp_geo
define geokolom=SHAPE
DROP INDEX &geotabel._SP_IDX;
drop table &geotabel.;
DELETE FROM USER_SDO_GEOM_METADATA WHERE TABLE_NAME=Upper('&geotabel');
create table &geotabel.(rowid1 rowid, rowid2 rowid, datasrcid number(22,11), datasrcid2 number(22,11), shape SDO_geometry);
INSERT INTO USER_SDO_GEOM_METADATA VALUES
  (Upper('&geotabel'),Upper('&geokolom'), MDSYS.SDO_DIM_ARRAY
  (MDSYS.SDO_DIM_ELEMENT('X', &allLEFT, &allRIGHT, &accX),
    MDSYS.SDO_DIM_ELEMENT('Y', &allBOTTOM, &allTOP, &accY) ), &SRID );
CREATE INDEX &geotabel._SP_IDX ON &geotabel. ( &geokolom ) INDEXTYPE IS MDSYS.SPATIAL_INDEX;
CREATE INDEX &geotabel._rowid1_IDX ON &geotabel. ( rowid1 );
CREATE INDEX &geotabel._rowid2_IDX ON &geotabel. ( rowid2 );
CREATE INDEX &geotabel._datasrcid_IDX ON &geotabel. ( datasrcid );
CREATE INDEX &geotabel._datasrcid2_IDX ON &geotabel. ( datasrcid2 );

-- WORKAROUND oracle Bug No. 4729792 "HEAVY DML ON RTREE INDEX CAN LEAD TO ORA-13236 ERRORS "
-- system grant
update mdsys.SDO_INDEX_METADATA_TABLE 
   set SDO_DML_BATCH_SIZE = 1 
   where sdo_index_owner = Upper('geouser') and SDO_DML_BATCH_SIZE<>1;
commit;

/* package declarations
-------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------
*/

CREATE OR REPLACE PACKAGE GEO AS

  -- Declare externally visible types, cursor, exception
  g_version varchar2(30) := 'Version 4.04, 201112';
  g_timeout number := 20;
  g_SRID number := &SRID.;
  g_allLEFT number := &allLEFT.;
  g_allBOTTOM number := &allBOTTOM.;
  g_allRIGHT number := &allRIGHT.;
  g_allTOP number := &allTOP.;
  g_accX number := &accX.;
  g_accY number := &accY.;  
  -- accT>=accX, when accT=accX then simplify is disabled
  g_accT number := 0.00000005;

  /*  Script_log usage
       - start:
       pid := script_log_start( 1000, 'description' );
       - progress:
       script_log_add( pid, 500 );
       - other logmessages:
       script_log_add(null,null, 'omschrijving' );

    Progress query:
      select LAMIN.Step as Stap
        , TO_CHAR(LAMIN.logtime,'dd mon HH24:MI:ss') as start_tijd_stap
        , TO_CHAR(24*60*(LAMAX.logtime-LAMIN.logtime),'9990.0') as min_bezig
        , Decode(LAMIN.todo-LAMAX.todo,0,NULL
          ,TO_CHAR(24*60*(LAMAX.logtime-LAMIN.logtime) * LAMIN.todo/(LAMIN.todo-LAMAX.todo),'9990.0')) as min_nodig
        , LAMIN.msg as Omschrijving
      from script_log LAMIN, script_log LAMAX
      where LAMIN.ID=(select min(ID) from script_log where step=LAMIN.step group by step)
        and LAMAX.ID=(select max(ID) from script_log where step=LAMIN.step group by step)
      order by LAMIN.step desc;       
  */

  /*  drop_if_exists('objecttype','object');

    drops object (index, table, procedure) if it exists, providing 
     that the user/schema has rights to drop this type of object using a
     pl/sql procedure
     
    examples:
      create table a_test(id number, text varchar2(10));
      create index a_index on a_test(id);
      call drop_if_exists('INDEX','a_index');
      call drop_if_exists('TABLE','a_test');
  */
  procedure drop_if_exists( objecttype1 VARCHAR2, objectname1 VARCHAR2);
    
  -- geo functions:
  
  /* getShapeType(geometry)
  returns:
     0 null geometry
     1 point(s)
     2 lines(s)
     3 polygon(s)
    40 collection null geometry
    41 collection only points
    42 collection only points or lines
    43 collection points, lines or polygons, note: sdo_gtype=2004 allows a mix of geometry types
  */
  function getShapeType(geom SDO_Geometry) return number;

  /* geoDifference(table1,table2) : tmp_geoResult := table1 - table2
  call geoDifference('MDM_METADATA where PROCESNAME LIKE ''%Create ENC%''', 'TOP10_BUFFER');
  call geoDifference(NULL, 'table3');
  -- SDO_UTIL.GETNUMVERTICES(shape)

  geoMode +1  nologging
  */
  PROCEDURE geoDifference(table1 varchar, table2 varchar, geoMode integer default 0);

  /* geoIntersect(table1,table2) : tmp_geoResult := table1 INTERSECT table2
  call geoIntersect('MDM_METADATA where PROCESNAME LIKE ''%Create ENC%''','TOP10_BUFFER_KB');

  geoMode +1 nologging
  geoMode +2 union all geometries that where intersected with 1 geometry of table1
  geoMode +4 do not first truncate result table
  geoMode +32 result in tmp_geo instead of tmp_geoResult
  */
  PROCEDURE geoIntersect(table1 varchar, table2 varchar, geoMode integer default 0);

  /* geoSelect(table1,table2) : tmp_geoResult := select geometries from table2 using table1
  call geoSelect('MDM_METADATA where PROCESNAME LIKE ''%Create ENC%''','TOP10_BUFFER_KB');

  geoMode  +1 nologging
  geoMode  +2 only geometries that are completely contained 
  geoMode  +4 do not first truncate result table
  geoMode  +8 only select the lines of the geometries
  geoMode +16 only geometries that interact with the borders. NB: do not use with +2
  geoMode +32 result in tmp_geo instead of tmp_geoResult
  */
  PROCEDURE geoSelect(table1 varchar, table2 varchar, geoMode integer default 0);

  /* geoPartition(table,width,height)
   tmp_geoResult := opdeling van vlakken in table in rechthoekige vlakken van maximaal width x height
   geoMode +1 nologging
  */
  PROCEDURE geoPartition(
    table1 varchar
    , v_width number
    , v_height number
    , geoMode integer default 0);

  /* map collection (2004) to simpletype geometry by removing other types
     keep_type=polygon(3), polyline(2) or point(1)
  */
  FUNCTION geoCollection_to_multi( Geo IN Mdsys.Sdo_geometry, keep_type number, geoMode integer default 0) RETURN Mdsys.Sdo_geometry;

  /* geoSelectWithCover('sa_geotable','geotable', sdo_cover_geometry);
    geoMode +1 select only when fully within p_geom_cover
  */
  procedure geoSelectWithCover(
     p_from_table varchar2
   , p_to_table varchar2
   , p_shape_type number --1=point, 2=line, 3=polygon
   , p_geom_cover MDSYS.SDO_geometry
   , geoMode integer default 0);

  /* geoCheck_sp_idx(NULL, 1);
  chkMode 1 :  rebuild all indexes

  Note: user needs system privileges:
  - grant update on mdsys.SDO_INDEX_METADATA_TABLE to geouser;
  - grant create any table to geouser;
  */
  procedure geoCheck_sp_idx(pid integer, chkMode integer default 0);
  
END GEO;
/

  
/* functions body
-------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------
*/

CREATE OR REPLACE PACKAGE BODY GEO AS

function script_log_start(
  v_todo NUMBER
  , v_msg varchar
  ) RETURN NUMBER IS PRAGMA AUTONOMOUS_TRANSACTION;
BEGIN
  DECLARE
    currstep NUMBER default NULL;
  BEGIN    
    select s_script_log_step.nextval into currstep from dual;
    INSERT INTO script_log(todo, step, lvl, msg) VALUES (v_todo, currstep, 0, v_msg);
    INSERT INTO script_log(todo, step, lvl, msg) VALUES (v_todo, currstep, 1, v_msg);
    COMMIT;
    RETURN currstep;
  END;
EXCEPTION WHEN OTHERS THEN 
  ROLLBACK;
  RETURN NULL;
END;

procedure script_log_add(
  v_pid NUMBER
  , v_todo NUMBER
  , v_msg varchar default NULL
  ) IS PRAGMA AUTONOMOUS_TRANSACTION;
BEGIN
  DECLARE
    currpid NUMBER;
  BEGIN
    IF NOT v_pid IS NULL and v_msg IS NULL THEN
      -- voortgang/einde: select script_log_add(21, 500) from dual;
      update script_log set todo=v_todo, logtime=sysdate where lvl=1 and step=v_pid;
    ELSIF NOT v_pid IS NULL and NOT v_msg IS NULL THEN
      -- voortgang/einde: select script_log_add(21, 500, 'omschrijving') from dual;
      INSERT INTO script_log(todo, step, lvl, msg) VALUES (v_todo, v_pid, 2, v_msg);
    ELSE
      -- overige logmeldingen
      select s_script_log_step.nextval into currpid from dual;
      INSERT INTO script_log(todo, step, lvl, msg) VALUES (0, currpid, 3, v_msg);    
    END IF;
    COMMIT;
  END;
EXCEPTION WHEN OTHERS THEN 
  ROLLBACK;
END;

procedure drop_if_exists( objecttype1 VARCHAR2, objectname1 VARCHAR2) as
BEGIN
  DECLARE
    i number;
  BEGIN
    select count(*) into i from user_objects  
      where object_type=Upper(objecttype1) and object_name=Upper(objectname1);
    IF i=1 THEN
      EXECUTE IMMEDIATE 'drop '||objecttype1||' '||objectname1;
    END IF;
  END;
END;

-- geo Functies
function getShapeType(geom SDO_Geometry) return number is
begin 
  declare  
    geom_dim number default 0;
    i integer;
  begin
    case 
      when geom is NULL then --null
        geom_dim:=0;
      when geom.sdo_gtype=2001 or geom.sdo_gtype=2005 then --point(s)
        geom_dim:=1;
      when geom.sdo_gtype=2002 or geom.sdo_gtype=2006 then --line(s)
        geom_dim:=2;
      when geom.sdo_gtype=2003 or geom.sdo_gtype=2007 then --polygon(s)
        geom_dim:=3;
      when geom.sdo_gtype=2004 then --collection
      begin
        -- look into SDO_ELEM_INFO
        geom_dim:=40;
        i:=2;
        LOOP
          EXIT WHEN i>geom.SDO_ELEM_INFO.COUNT;
          
          if geom_dim<41 and geom.SDO_ELEM_INFO( i ) = 1 then -- point
            geom_dim:=41; 
          elsif geom_dim<42 and geom.SDO_ELEM_INFO( i ) = 2 then -- line
            geom_dim:=42; 
          elsif (geom.SDO_ELEM_INFO( i ) = 1003 or geom.SDO_ELEM_INFO( i ) = 2003) then -- polygon
            geom_dim:=43; 
            i:=geom.SDO_ELEM_INFO.COUNT;
          end if;
          i:=i+3;
        END LOOP;
      end;
    
    end case;

    return geom_dim;
  end;
end;

PROCEDURE geoDifference(table1 varchar, table2 varchar, geoMode integer default 0) is
BEGIN
  DECLARE
    pid integer;
    p_start date;
    i integer;
    TYPE geoCurTyp  IS REF CURSOR;
    CURSOR geo_cursor1 is  select shape from tmp_geoResult FOR UPDATE OF shape;
    geo_cursor2 geoCurTyp;
    shape1 tmp_geoResult.shape%TYPE;
    shape2 tmp_geoResult.shape%TYPE;
  BEGIN
    IF NOT table1 IS NULL THEN
      execute immediate 'truncate table tmp_geoResult';
      execute immediate 'insert into tmp_geoResult(rowid1,datasrcid,shape) select rowid,datasrcid,shape from '||table1;
    END IF;
    
    select count(*) into i from tmp_geoResult;

    IF bitand(geoMode,1)=0 THEN
      pid:=geo.script_log_start(null,'  geoDifference("'|| nvl(table1,'null') ||'","'|| table2 ||'",'||to_char(nvl(geoMode,0))||')');
      p_start := sysdate;
    END IF;
    
    OPEN geo_cursor1;
    LOOP
      FETCH geo_cursor1 INTO shape1;
      EXIT WHEN geo_cursor1%NOTFOUND;
      
      OPEN geo_cursor2 FOR 'select shape from ' || table2 
        || ' where SDO_ANYINTERACT(shape, :j)= ''TRUE''' USING shape1;
      LOOP
        FETCH geo_cursor2 INTO shape2;
        EXIT WHEN geo_cursor2%NOTFOUND;
        shape1 := SDO_GEOM.SDO_DIFFERENCE(shape1, shape2, g_accX);

        IF bitand(geoMode,1)=0 and 24*60*(sysdate-p_start)>g_timeout THEN
          geo.script_log_add(pid,i,'  geoDifference afgebroken, timeout '||to_char( g_timeout )||' minuten');
          raise_application_error(-20002, 'geoDifference afgebroken, timeout '||to_char( g_timeout )||' minuten');
        END IF;
        EXIT WHEN shape1 is NULL;
      END LOOP;
      CLOSE geo_cursor2;
      
      update tmp_geoResult set SHAPE=shape1 where current of geo_cursor1;

      IF bitand(geoMode,1)=0 THEN
        i:=i-1;
        geo.script_log_add(pid,i);
      END IF;
    END LOOP;

    CLOSE geo_cursor1;

    IF bitand(geoMode,1)=0 THEN
      geo.script_log_add(pid,0);
    END IF;
  END;
END;

PROCEDURE geoIntersect(table1 varchar, table2 varchar, geoMode integer default 0) is
BEGIN
  DECLARE
    pid integer;
    i integer;
    TYPE geoCurTyp IS REF CURSOR;
    geo_cursor1 geoCurTyp;
    geo_cursor2 geoCurTyp;
    geom MDSYS.SDO_geometry;
    shape1 tmp_geoResult.shape%TYPE;
    shape2 tmp_geoResult.shape%TYPE;
    rowid1 tmp_geoResult.rowid1%TYPE;
    rowid2 tmp_geoResult.rowid1%TYPE;
    datasrcid1 tmp_geoResult.datasrcid%TYPE;
    datasrcid2 tmp_geoResult.datasrcid%TYPE;
  BEGIN
    IF bitand(geoMode,4)=0 THEN
      IF bitand(geoMode,32)=0 THEN
        execute immediate 'truncate table tmp_geoResult';
      ELSE
        execute immediate 'truncate table tmp_geo';
      END IF;
    END IF;
    
    OPEN geo_cursor1 FOR 'select rowid,shape,datasrcid from '||table1;    
    
    execute immediate 'select count(*) from '||table1 into i;

    IF bitand(geoMode,1)=0 THEN
      pid:=geo.script_log_start(i,'  geoIntersect("'|| table1 ||'","'|| table2 ||'",'||to_char(nvl(geoMode,0))||')');
    END IF;
    
    geom := NULL;
    LOOP 
      FETCH geo_cursor1 INTO rowid1,shape1,datasrcid1;
      EXIT WHEN geo_cursor1%NOTFOUND;

      OPEN geo_cursor2 FOR 'select rowid,shape,datasrcid from ' || table2 
        || ' where SDO_ANYINTERACT(shape, :j)= ''TRUE''' USING shape1;
      LOOP
        FETCH geo_cursor2 INTO rowid2, shape2, datasrcid2;
        EXIT WHEN geo_cursor2%NOTFOUND;

        IF bitand(geoMode,2)=0 THEN
          -- no union
          IF bitand(geoMode,32)=0 THEN
            insert into tmp_geoResult(rowid1, rowid2, datasrcid, datasrcid2, shape) values ( rowid1, rowid2, datasrcid1, datasrcid2, SDO_GEOM.SDO_INTERSECTION(shape1, shape2, g_accX));        
          ELSE
            insert into tmp_geo(rowid1, rowid2, datasrcid, datasrcid2, shape) values ( rowid1, rowid2, datasrcid1, datasrcid2, SDO_GEOM.SDO_INTERSECTION(shape1, shape2, g_accX));        
          END IF;
        ELSE
          geom:=SDO_GEOM.SDO_UNION(geom, SDO_GEOM.SDO_INTERSECTION(shape1, shape2, g_accX), g_accX);
        END IF;
      
      END LOOP;
      CLOSE geo_cursor2;
      
      IF bitand(geoMode,2)>0 and NOT geom IS NULL THEN
        IF bitand(geoMode,32)=0 THEN
          insert into tmp_geoResult(rowid1, rowid2, datasrcid, datasrcid2, shape) values ( rowid1, rowid2, datasrcid1, datasrcid2, geom);
        ELSE
          insert into tmp_geo(rowid1, rowid2, datasrcid, datasrcid2, shape) values ( rowid1, rowid2, datasrcid1, datasrcid2, geom);
        END IF;
        geom := NULL;
      END IF;

      IF bitand(geoMode,1)=0 THEN
        i:=i-1;
        geo.script_log_add(pid,i);
      END IF;
    END LOOP;
    CLOSE geo_cursor1;
    
    IF bitand(geoMode,1)=0 THEN
      geo.script_log_add(pid,0);
    END IF;
  END;
END;

PROCEDURE geoSelect(table1 varchar, table2 varchar, geoMode integer default 0) is
BEGIN
  DECLARE
    pid integer;
    i integer;
    sql2 varchar(200);
    TYPE geoCurTyp  IS REF CURSOR;
    geo_cursor1 geoCurTyp;
    geo_cursor2 geoCurTyp;
    shape1 tmp_geoResult.shape%TYPE;
    shape2 tmp_geoResult.shape%TYPE;
    rowid1 tmp_geoResult.rowid1%TYPE;
    rowid2 tmp_geoResult.rowid2%TYPE;
    datasrcid1 tmp_geoResult.datasrcid%TYPE;
    datasrcid2 tmp_geoResult.datasrcid%TYPE;
  BEGIN
    IF bitand(geoMode,4)=0 THEN
      IF bitand(geoMode,32)=0 THEN
        execute immediate 'truncate table tmp_geoResult';
      ELSE
        execute immediate 'truncate table tmp_geo';
      END IF;
    END IF;
    
    IF bitand(geoMode,1)=0 THEN
      execute immediate 'select count(*) from '||table1 into i;
      pid:=geo.script_log_start(i,'  geoSelect("'|| table1 ||'","'|| table2 ||'",'||to_char(nvl(geoMode,0))||')');
    END IF;
    
    IF instr(lower(table2),' where ')=0 THEN
      sql2:='select rowid, datasrcid, shape from ' || table2 || ' where ';
    ELSE
      sql2:='select rowid, datasrcid, shape from ' || table2 || ' and ';
    END IF;
    
    IF bitand(geoMode,2)=2 THEN
      sql2:=sql2 || 'SDO_RELATE(shape, :i,''mask=inside+coveredby+overlapbdyintersect+equal'')= ''TRUE''';
      -- mask=inside+coveredby+overlapbdyintersect+equal+overlapbdydisjoint+on
    ELSIF bitand(geoMode,16)=16 THEN
      sql2:=sql2 || 'SDO_RELATE(shape, :i,''mask=overlapbdyintersect'')= ''TRUE''';
      -- mask=overlapbdyintersect+overlapbdydisjoint+on
    ELSE
      sql2:=sql2 || 'SDO_ANYINTERACT(shape, :i)= ''TRUE''';
    END IF;
    
    OPEN geo_cursor1 FOR 'select rowid,datasrcid,shape from '||table1;    
    LOOP 
      FETCH geo_cursor1 INTO rowid1,datasrcid1,shape1;
      EXIT WHEN geo_cursor1%NOTFOUND;

      OPEN geo_cursor2 FOR sql2 USING shape1;    
      LOOP
        FETCH geo_cursor2 INTO rowid2, datasrcid2, shape2;
        EXIT WHEN geo_cursor2%NOTFOUND;
        
        IF bitand(geoMode,32)=0 THEN
          IF bitand(geoMode,8)=0 THEN
            insert into tmp_geoResult(rowid1, rowid2, datasrcid, datasrcid2, shape) 
              values ( rowid1, rowid2, datasrcid1, datasrcid2, shape2 );
          ELSE
            insert into tmp_geoResult(rowid1, rowid2, datasrcid, datasrcid2, shape) 
              values ( rowid1, rowid2, datasrcid1, datasrcid2, SDO_UTIL.POLYGONTOLINE( shape2 ) );
          END IF;
        ELSE
          IF bitand(geoMode,8)=0 THEN
            insert into tmp_geo(rowid1, rowid2, datasrcid, datasrcid2, shape) 
              values ( rowid1, rowid2, datasrcid1, datasrcid2, shape2 );
          ELSE
            insert into tmp_geo(rowid1, rowid2, datasrcid, datasrcid2, shape) 
              values ( rowid1, rowid2, datasrcid1, datasrcid2, SDO_UTIL.POLYGONTOLINE( shape2 ) );
          END IF;
        END IF;

      END LOOP;
      CLOSE geo_cursor2;
      
      IF bitand(geoMode,1)=0  THEN
        i:=i-1;
        geo.script_log_add(pid,i);
      END IF;
    END LOOP;
    CLOSE geo_cursor1;

    IF bitand(geoMode,1)=0 THEN
      geo.script_log_add(pid,0);
    END IF;
  END;
END;

PROCEDURE geoPartition(
  table1 varchar
  , v_width number
  , v_height number
  , geoMode integer default 0) is
BEGIN
  DECLARE
    i integer;
    pid integer;
    TYPE geoCurTyp  IS REF CURSOR;
    geo_cursor1 geoCurTyp;
    shape1 tmp_geoResult.shape%TYPE;
    shape2 tmp_geoResult.shape%TYPE;
    geom tmp_geoResult.shape%TYPE;
    rowid1 tmp_geoResult.rowid1%TYPE;
    datasrcid1 tmp_geoResult.datasrcid%TYPE;
    xi integer; 
    yi integer;
    xcount integer; 
    ycount integer;
    x2 number; 
    y2 number;
    v_left number; 
    v_top number; 
    v_right number; 
    v_bottom number;
    width2 number; 
    height2 number;
    shape1_enccover_toporel varchar2(30);
  BEGIN
    -- '||nvl(v_width,'null')||','||nvl(v_height,'null')||')'); 

    IF bitand(geoMode,1)=0 THEN
      execute immediate 'select count(*) from '||table1 into i;
      pid:=geo.script_log_start(i,'  geoPartition("'||table1||'", '||to_char(v_width)||', '||to_char(v_height)||')');
    END IF;

    -- raster maken op basis van vlakken in table1 -> tmp_geoResult
    OPEN geo_cursor1 FOR 'select rowid,datasrcid,shape from '||table1;
    execute immediate 'truncate table tmp_geoResult';
    LOOP
      FETCH geo_cursor1 INTO rowid1, datasrcid1, shape1;
      EXIT WHEN geo_cursor1%NOTFOUND;

      -- v_left < v_right, v_bottom < v_top, v_width, v_height
      -- v_width := 0.0415
      -- v_height := 0.0277
      -- ORA-06533: Subscript beyond count wanneer area(shape1)=0
      v_right := SDO_GEOM.SDO_MAX_MBR_ORDINATE(shape1, 1);
      v_left := SDO_GEOM.SDO_MIN_MBR_ORDINATE(shape1, 1);
      v_top := SDO_GEOM.SDO_MAX_MBR_ORDINATE(shape1, 2);
      v_bottom := SDO_GEOM.SDO_MIN_MBR_ORDINATE(shape1, 2);
      
      xcount := ROUND( ( v_right - v_left ) / v_width );
      ycount := ROUND( ( v_top - v_bottom ) / v_height );      
      IF xcount<1 THEN 
        xcount := 1; 
      END IF;
      IF ycount<1 THEN 
        ycount := 1; 
      END IF;

      width2 := ( v_right - v_left ) / xcount;
      height2 := ( v_top - v_bottom ) / ycount;
      
      for yi in 0 .. ycount-1 loop
        for xi in 0 .. xcount-1 loop

          IF xi < xcount-1 THEN
            x2 := v_left+(xi+1)*width2;
          ELSE
            x2 := v_right;
          END IF;
          
          IF yi < ycount-1 THEN
            y2 := v_bottom+(yi+1)*height2;
          ELSE
            y2 := v_top;          
          END IF;

          shape2 := MDSYS.SDO_GEOMETRY(2003, g_SRID, NULL,
              MDSYS.SDO_ELEM_INFO_ARRAY(1,1003,3),
              MDSYS.SDO_ORDINATE_ARRAY( v_left+xi*width2, v_bottom+yi*height2, x2, y2));

          BEGIN
            shape1_enccover_toporel:=SDO_GEOM.RELATE(shape2, 'DETERMINE', shape1, g_accX);
            IF shape1_enccover_toporel in ('OVERLAPBDYINTERSECT', 'OVERLAPBDYDISJOINT','INSIDE','COVEREDBY','EQUAL') THEN 
              geom := SDO_GEOM.SDO_INTERSECTION(shape1, shape2, g_accX);
              IF NOT geom is NULL THEN
                insert into tmp_geoResult(rowid1, datasrcid, shape) values ( rowid1, datasrcid1, geom);
              END IF;
            END IF;
          EXCEPTION WHEN OTHERS THEN
            -- skip 'impossible' geometries
            geo.script_log_add(pid,i, '  geoPartition() "'||DBMS_UTILITY.FORMAT_ERROR_BACKTRACE || SQLERRM||'"');
            geom := NULL;
          END;
        end loop;
      end loop;

      IF bitand(geoMode,1)=0 THEN
        i:=i-1;
        geo.script_log_add(pid,i);
      END IF;
    END LOOP;
    
    IF bitand(geoMode,1)=0 THEN
      geo.script_log_add(pid, 0);
    END IF;
  END;
END;

FUNCTION geoCollection_to_multi( Geo IN Mdsys.Sdo_geometry, keep_type number, geoMode integer default 0) RETURN Mdsys.Sdo_geometry
AS  v_geo Mdsys.Sdo_geometry;
  v_gtype PLS_INTEGER;
  v_numdims PLS_INTEGER;
  v_offset PLS_INTEGER;
  v_elem_type PLS_INTEGER;
  v_interpretation PLS_INTEGER;
  v_numsubelements PLS_INTEGER;
  v_numgeometry PLS_INTEGER;
  v_nextoffset PLS_INTEGER;
  i PLS_INTEGER;
  j PLS_INTEGER;
BEGIN
  IF Geo is NULL THEN
    RETURN NULL;
  END IF;  
  
  IF NOT keep_type in (1,2,3) THEN
    raise_application_error(-20003, 'Value of KEEP_TYPE should be 1,2 or 3');
  END IF;

  v_geo := NULL;
  
  -- new geometry
  v_geo := Mdsys.Sdo_geometry( Geo.Sdo_gtype, Geo.Sdo_srid, NULL, NULL, NULL );
  
  -- geometry type
  -- 1-POINT 5-MULTIPOINT 2-LINE or CURVE 6-MULTILINE or MULTICURVE 3-POLYGON 7-MULTIPOLYGON 4-COLLECTION
  v_gtype := mod(Geo.Sdo_gtype, 100);
  
  IF NOT v_gtype=4 THEN
    RETURN Geo;
  END IF;
  
  -- number of dimensions
  v_numdims := Trunc( v_gtype / 1000 );

  -- elements
  IF NOT Geo.Sdo_elem_info IS NULL THEN
    
    -- new elem_info
    v_geo.Sdo_elem_info := Mdsys.Sdo_elem_info_array( );
    
    -- new ordinates
    v_geo.Sdo_ordinates := Mdsys.Sdo_ordinate_array( );
    
    -- init
    v_numsubelements:=0;
    v_numgeometry:=0;
    
    -- values
    i:=1;
    WHILE i <= Geo.Sdo_elem_info.COUNT-2 LOOP
      v_offset:=Geo.Sdo_elem_info( i );
      i:=i+1;
      
      -- element type, interpretation. Consult oracle documentation for detailed information.
      v_elem_type:=Geo.Sdo_elem_info( i );
      i:=i+1;
      
      v_interpretation:=Geo.Sdo_elem_info( i );
      i:=i+1;
      
      -- subelements (not supported)
      IF (v_numsubelements=0) AND (v_interpretation>1) AND
          v_elem_type in (4 -- linestring
          , 1005, 2005 ) -- polygon exterior/interior
      THEN
        --raise_application_error(-20004, 'subelements not supported');
        RETURN NULL;
      END IF;
            
      IF (keep_type=1 AND v_elem_type=1)
      OR (keep_type=2 AND v_elem_type in (2,4))
      OR (keep_type=3 AND v_elem_type in (1003,2003)) THEN   
      
        IF v_elem_type <> 2003 THEN
          v_numgeometry := v_numgeometry+1;
        END IF;

        v_geo.Sdo_elem_info.EXTEND( 1 );
        v_geo.Sdo_elem_info( v_geo.Sdo_elem_info.COUNT ) := v_geo.Sdo_ordinates.COUNT+1;
        v_geo.Sdo_elem_info.EXTEND( 1 );
        v_geo.Sdo_elem_info( v_geo.Sdo_elem_info.COUNT ) := v_elem_type;
        v_geo.Sdo_elem_info.EXTEND( 1 );
        v_geo.Sdo_elem_info( v_geo.Sdo_elem_info.COUNT ) := v_interpretation;
  
        -- next offset (if any)
        IF i <= Geo.Sdo_elem_info.COUNT-2 THEN
          v_nextoffset:=Geo.Sdo_elem_info( i );
        ELSE
          -- read remaining coordinates
          v_nextoffset:=Geo.Sdo_ordinates.COUNT+1;
        END IF;
        
        -- read element coordinates
        j:=v_offset;
        WHILE j < v_nextoffset and j <= Geo.Sdo_ordinates.COUNT LOOP
          v_geo.Sdo_ordinates.EXTEND( 1 );
          v_geo.Sdo_ordinates( v_geo.Sdo_ordinates.COUNT ) := Geo.Sdo_ordinates( j );        
          j:=j+1;
        END LOOP;
      --ELSE element
      END IF;
      
    END LOOP;
  END IF;
  
  --gtype
  IF NOT v_geo.sdo_elem_info is NULL THEN
    v_geo.sdo_gtype := 2000 + keep_type;
    IF v_numgeometry>1 THEN
     v_geo.sdo_gtype :=  v_geo.sdo_gtype + 4;
    ELSIF v_numgeometry=0 THEN
    v_geo := NULL;
    END IF;
  END IF;
  
  --check
  IF NOT v_geo.sdo_gtype in (2001,2002,2003,2005,2006,2007) THEN
    --raise_application_error(-20005, 'type not supported');
    RETURN NULL;
  END IF;

  RETURN v_geo;
END;

procedure geoSelectWithCover(
  p_from_table varchar2
  , p_to_table varchar2
  , p_shape_type number
  , p_geom_cover MDSYS.SDO_geometry
  , geoMode integer default 0) is
BEGIN
  DECLARE
    TYPE geoCurTyp IS REF CURSOR;
    geo_cursor1 geoCurTyp;
    r_id number(11);
    r_datasrcid number(22,11); 
    r_shape MDSYS.SDO_geometry; 
    v_shape MDSYS.SDO_geometry; 
    shape1_enccover_toporel varchar2(30);
    v_column varchar2(30);
    v_columns varchar2(500);  
  BEGIN
    -- build query columns, select columns that appear in both tables, exclude the SHAPE column
    v_columns := NULL;
    OPEN geo_cursor1 FOR 'select from_table.cname as cname 
      from col from_table, col to_table 
        where from_table.tname=Upper('''||p_from_table||''')
        and to_table.tname=Upper('''||p_to_table||''')
        and from_table.cname=to_table.cname
        and from_table.cname<>''SHAPE''
        and from_table.cname<>''ID''
        order by from_table.colno';
    LOOP
      FETCH geo_cursor1 INTO v_column;
      EXIT WHEN geo_cursor1%NOTFOUND;    
      v_columns := v_columns || ',' || v_column;
    END LOOP;
    CLOSE geo_cursor1;

    -- select data
    OPEN geo_cursor1 FOR 'select id,datasrcid,shape from '||p_from_table
      || ' where not shape is null and SDO_RELATE(shape, :j, ''mask=inside+coveredby+overlapbdyintersect+overlapbdydisjoint+equal'')= ''TRUE''' USING p_geom_cover;
    LOOP
      FETCH geo_cursor1 INTO r_id, r_datasrcid, r_shape;
      EXIT WHEN geo_cursor1%NOTFOUND;

      IF mod(geo.getShapeType(r_shape),10)=p_shape_type THEN

         -- simplify geometry (Douglas-Peucker algorithm)
        IF g_accT > g_accX THEN
          v_shape:=SDO_UTIL.SIMPLIFY(r_shape, g_accT, g_accX);
          -- validate
          IF SDO_GEOM.VALIDATE_GEOMETRY_WITH_CONTEXT(v_shape, g_accX)<>'TRUE' THEN
            -- simplified shape is not valid, use the original shape
            v_shape := r_shape;
          END IF;
        ELSE
          v_shape := r_shape;
        END IF;

        IF mod(geo.getShapeType(v_shape),10)=p_shape_type THEN
          shape1_enccover_toporel:=SDO_GEOM.RELATE(p_geom_cover, 'DETERMINE', v_shape, g_accX);
          IF shape1_enccover_toporel in ('OVERLAPBDYINTERSECT', 'OVERLAPBDYDISJOINT') THEN 
            IF  bitand(geoMode,1)=0 THEN
              -- intersects cover, use partly if the part fits the geometrytype
              v_shape:=SDO_GEOM.SDO_INTERSECTION(v_shape, p_geom_cover, g_accX);
              IF mod(geo.getShapeType(v_shape),10)=p_shape_type THEN
                execute immediate 'insert into '||p_to_table||'(objectid,id,shape,table_name,enc'||v_columns||') select s_objectid.nextval,s_objectid.currval,:s,:t,:e'||v_columns||' from '||p_from_table||' where ID=:f' using v_shape, p_from_table, r_id;
              END IF;
            END IF;
          ELSE -- use whole
            execute immediate 'insert into '||p_to_table||'(objectid,id,shape,table_name,enc'||v_columns||') select s_objectid.nextval,s_objectid.currval,shape,:t,:e'||v_columns||' from '||p_from_table||' where ID=:f' using p_from_table, r_id;
          END IF;
        END IF; -- shapetype
      END IF; -- shapetype      
    END LOOP;
    CLOSE geo_cursor1;
  END;
END;

procedure geoCheck_sp_idx(pid integer, chkMode integer default 0) is
BEGIN
  DECLARE
    pid_null integer default NULL;
  BEGIN
    IF pid is NULL THEN
      pid_null:=geo.script_log_start(1, 'import_check_sp_idx('||TO_CHAR(chkMode)||')');
    END IF;
    
    -- check, create, re-create, rebuild spatial indexes
    FOR rec IN ( select * from (
      select md.table_name as tname, md.column_name as cname, idx.index_name as iname
        , decode( idx.STATUS || ',' || idx.DOMIDX_STATUS || ',' || idx.DOMIDX_OPSTATUS, 'VALID,VALID,VALID', 'VALID', ',,', 'ONLYMD', 'INVALID') as istate
        from user_sdo_geom_metadata md, user_ind_columns cidx, user_indexes idx
        where md.table_name=cidx.table_name(+)
        and md.column_name=cidx.column_name(+)
        and cidx.index_name=idx.index_name(+)
        and idx.ityp_name(+)='SPATIAL_INDEX'
        and not md.table_name like '%$%'
      union
        select col.tname as tname, col.cname as cname, NULL as iname, 'NONE' as istate 
        from col
        where col.coltype like '%SDO_GEOMETRY%'
        and not tname like '%$%'
        and NOT (col.tname, col.cname) in (
          select md.table_name, md.column_name
          from user_sdo_geom_metadata md )
      ) order by tname,cname,iname )
    LOOP
      IF rec.istate='ONLYMD' THEN
        geo.script_log_add(pid, 1, 'index "'||substr(rec.tname,1,30-7)||'_SP_IDX" create');
        -- create index
        execute immediate 'CREATE INDEX '||substr(rec.tname,1,30-7)||'_SP_IDX ON '||rec.tname||' ( '||rec.cname||' ) INDEXTYPE IS MDSYS.SPATIAL_INDEX  PARAMETERS(''tablespace=geouser_INDEX'')';
      ELSIF rec.istate='INVALID' THEN
        geo.script_log_add(pid, 1, 'index "'||rec.iname||'" re-create');
        -- drop force, re-create
        execute immediate 'DROP INDEX '||rec.iname||' FORCE';
        execute immediate 'CREATE INDEX '||rec.tname||'_SP_IDX ON '||rec.tname||' ( '||rec.cname||' ) INDEXTYPE IS MDSYS.SPATIAL_INDEX  PARAMETERS(''tablespace=geouser_INDEX'')';
      ELSIF bitand(chkMode,1)=1 and not rec.iname is NULL THEN
        geo.script_log_add(pid, 1, 'index "'||rec.iname||'" rebuild');
        -- otherwise rebuild the index
        execute immediate 'ALTER INDEX '||rec.iname||' rebuild';
      END IF;
    END LOOP ;
    
    -- WORKAROUND oracle Bug No. 4729792 "HEAVY DML ON RTREE INDEX CAN LEAD TO ORA-13236 ERRORS "
    -- (system) grant update on mdsys.SDO_INDEX_METADATA_TABLE to geouser;
    update mdsys.SDO_INDEX_METADATA_TABLE set SDO_DML_BATCH_SIZE = 1 
      where sdo_index_owner = 'GEOUSER' and SDO_DML_BATCH_SIZE<>1;
    commit;
    
    IF not pid_null is NULL THEN
      geo.script_log_add(pid_null,0);
    END IF;
    
  END;
END;

END;
/