Dienstag, 30. November 2010

Lauter kleine Helferlein

Heute soll es darum gehen, Geometrien, die nicht Simple Feature-konform sind, soweit möglich zu korrigieren.
Für eine solche Korrektur steht seit 11gR1 die Funktion RECTIFY_GEOMETRY im Package SDO_UTIL bereit, welche die korrigierte (rektifizierte) Geometrie zurückliefert.
Ein Blick in die Oracle Spatial Online-Dokumentation verrät dabei, welche Fehler korrigierbar sind:
  • doppelte Stützpunkte (bezogen auf die Toleranz)
  • sich selbst schneidende Polygone
  • falsche Orientierung von innerem oder/und äusseren Ring(en) eines Polygons
Dort findet sich ebenso das Format beschrieben.
SDO_UTIL.RECTIFY_GEOMETRY(
  geom      IN SDO_GEOMETRY,
  tolerance IN NUMBER
) RETURN SDO_GEOMETRY;

Wie kann man nun diese Funktion in Anwendung bringen? Schauen wir uns das einfach mal genauer an.
Zunächst wird eine Tabelle angelegt, in welcher alle fehlerhaften Geometrien registriert werden sollen.
drop table geometry_errors purge;

-- Fehlertabelle anlegen
create table geometry_errors (
  table_name      varchar2(30),
  column_name     varchar2(30),
  obj_rowid       rowid,
  geometry        sdo_geometry,
  tolerance       number,
  error_code      char(5),
  error_message   varchar2(256),
  error_context   varchar2(256),
  is_fixed        number(1) default 0
);
Die eigentliche Prozedur zum Überprüfen der Geometrien, kann dann wie nachfolgend aussehen. Die problematischen Geometrien darunter werden in die Fehlertabelle übernommen.
-- Prozedur zum Prüfen auf invalide Geometrien im gesamten Schema
-- Achtung: Kann einige Zeit in Anspruch nehmen!
declare
  DEFAULT_TOLERANCE   number := 0.0000005;
  COMMIT_FREQUENCY    number := 100;
  geom_cursor         sys_refcursor;
  v_diminfo           sdo_dim_array;
  v_srid              number;
  v_tolerance         number;
  v_rowid             rowid;
  v_geometry          sdo_geometry;
  v_num_rows          number;
  v_num_errors        number;
  v_error_code        char(5);
  v_error_message     varchar2(256);
  v_error_context     varchar2(256);
  v_status            varchar2(256);
begin

  -- Alle Tabellen mit SDO_GEOMETRY Spalten prozessieren
  for t in (
    select table_name, column_name
    from   user_tab_columns
    where  data_type = 'SDO_GEOMETRY'
    order by table_name, column_name )
  loop
  
    -- Lies Toleranz aus den Metadaten
    begin
      select diminfo, srid into v_diminfo, v_srid
      from user_sdo_geom_metadata
      where table_name = t.table_name
        and column_name = t.column_name;
    exception
      when no_data_found then
        v_diminfo := null;
        v_srid := null;
    end;
 
    -- Bei fehlenden Metadaten, Default-Toleranz annehmen
    if v_diminfo is null then
      v_tolerance := DEFAULT_TOLERANCE;
    else
      v_tolerance := v_diminfo(1).sdo_tolerance;
    end if;
 
    -- Geometrien prozessieren
    v_num_rows := 0;
    v_num_errors := 0;

    open geom_cursor for
      'select rowid,' || t.column_name || ' from ' || t.table_name;
    loop
      v_status := NULL;
   
      -- Hole Geometrie
      fetch geom_cursor into v_rowid, v_geometry;
        exit when geom_cursor%notfound;
      v_num_rows := v_num_rows + 1;
   
      -- Validiere Geometry
      v_status := sdo_geom.validate_geometry_with_context (v_geometry, v_tolerance);
   
      -- Wenn Fehler, dann protokollieren
      if v_status <> 'TRUE' then
   
        -- Fehler hochzählen
        v_num_errors := v_num_errors + 1;
  
        -- Format the error message
        if length(v_status) >= 5 then
          v_error_code := substr(v_status, 1, 5);
          v_error_message := sqlerrm(-v_error_code);
          v_error_context := substr(v_status,7);
        else
          v_error_code := v_status;
          v_error_message := null;
          v_error_context := null;
        end if;
  
        -- Fehler wegschreiben
        insert into geometry_errors (
          table_name,
          column_name,
          obj_rowid,
          geometry,
          tolerance,
          error_code,
          error_message,
          error_context
        )
        values (
          t.table_name,
          t.column_name,
          v_rowid,
          v_geometry,
          v_tolerance,
          v_error_code,
          v_error_message,
          v_error_context
        );
      end if;
   
      -- Commit entsprechend COMMIT_FREQUENCY
      if mod(v_num_rows,COMMIT_FREQUENCY) = 0 then
        commit;
      end if;
    end loop;
  end loop;
  
  -- Finales Commit
  commit;
end;
/
Zu beachten ist dabei, dass das Ausführen der Prozedur eine Weile in Anspruch nehmen kann, je nach Anzahl der zu überprüfenden Geometrien.
Bevor es an die eigentliche Korrektur geht, schauen wir uns zunächst die Fehlertabelle genauer an.
-- Fehler pro Tabelle, Error Code
select table_name, error_code, error_message, count(*)
from geometry_errors
group by table_name, error_code, error_message 
order by table_name, error_code;

-- Fehler pro Error Code, Tabelle
select error_code, error_message, table_name, count(*)
from geometry_errors
group by error_code, error_message, table_name
order by error_code, table_name;

-- Fehlerdetails
select table_name, obj_rowid, is_fixed, error_message, error_context
from geometry_errors
where is_fixed = 0
order by table_name, obj_rowid;
Für die Korrektur wird nun die Funktion SDO_UTIL.RECTIFY_GEOMETRY eingesetzt.
-- Prozedur zum Korrigieren der Geometrien
declare
  -- Fehlerbehandlung für nicht korrigierbare Geometrien
  -- "ORA-13199: the given geometry cannot be rectified"
  cannot_rectify exception;
  pragma exception_init(cannot_rectify, -13199);
  v_geometry_fixed sdo_geometry;

begin
  -- Invalide Geometrien prozessieren
  for e in (
    select rowid, table_name, column_name, obj_rowid, tolerance, geometry
    from geometry_errors
    where is_fixed = 0
    order by table_name, column_name
  )
  loop
    -- Wenn möglich, Geometrie korrigieren mit SDO_UTIL.RECTIFY_GEOMETRY. 
    -- Andernfalls:
    --   In 11.1.0.6: Funktion gibt NULL zurück
    --   In 11.1.0.7 und 11.2: Funktion wirft ORA-13199 Fehler
    begin
      v_geometry_fixed := sdo_util.rectify_geometry (e.geometry, e.tolerance);
    exception
      when cannot_rectify then
        v_geometry_fixed := null;
    end;

    if v_geometry_fixed is not null then
      -- Korrigiere Geometrie
      execute immediate 'update ' || e.table_name || ' set '|| e.column_name || ' = :g where rowid = :r'
        using v_geometry_fixed, e.obj_rowid;

      -- Kennzeichne korrigierte Geometrien in Fehlertabelle mit is_fixed = 1
      update geometry_errors set is_fixed = 1 where rowid = e.rowid;
      -- Alternativ ist auch DELETE möglich
      -- delete from geometry_error where rowid = e.rowid;
    end if;

  end loop;
end;
/
Was wurde korrgiert? Was nicht? Das Ergebnis kann am Ende in der Fehlertabelle überprüft werden.
-- Fehler nach Tabelle, Korrekturstatus und Error Code
select table_name, is_fixed, error_code, error_message, count(*)
from geometry_errors
group by table_name, is_fixed, error_code, error_message 
order by table_name, is_fixed, error_code;

-- Fehler nach Korrekturstatus und Error Code
select is_fixed, error_code, error_message, count(*)
from geometry_errors
group by is_fixed, error_code, error_message
order by is_fixed, error_code;

-- Details
select table_name, obj_rowid, is_fixed, error_message, error_context
from geometry_errors
order by is_fixed, table_name, obj_rowid;
Möglicherweise konnten nicht alle fehlerhaften Geometrien auf diese Weise korrigiert werden, aber ein erster Schritt in Richtung Qualitätsverbesserung ist getan.
Wer die Fehlertabelle nach erfolgter Korrektur nicht mehr benötigt, sollte übrigens nicht vergessen, diese wieder zu entfernen.
drop table geometry_errors purge;

Montag, 22. November 2010

GPX-Dateien nach SDO_GEOMETRY umwandeln

Jüngst habe ich für einen Kollegen einen GPS-Track in die Datenbank geladen. Ich hatte eine .gpx-Datei, die in etwa so aussieht ...
<gpx xmlns="http://www.topografix.com/GPX/1/1" tc2=" ...
   <trk>
       <name>2010-11-16T03:46:55Z</name>
       <trkseg>
           <trkpt lat="51.4799802" lon="7.9678522">
               <ele>270.7484131</ele>
               <time>2010-11-16T03:46:55Z</time>
           </trkpt>
           <trkpt lat="51.4799436" lon="7.9678679">
               <ele>264.0192871</ele>
               <time>2010-11-16T03:47:06Z</time>
           </trkpt>
           <trkpt lat="51.4797696" lon="7.9681059">
               <ele>264.4998779</ele>
               <time>2010-11-16T03:47:12Z</time>
           </trkpt>
           <trkpt lat="51.4797497" lon="7.9681919">
               <ele>264.9805908</ele>
               <time>2010-11-16T03:47:14Z</time>
           </trkpt>
Zunächst hatte ich ein wenig nach Converter-Software gegoogelt, wurde jedoch nicht so recht fündig (komische Lizenzbedingungen, keine Konvertierung nach GML oder KML, sondern nur in andere Formate). Da es aber um eine einfach zu verstehende XML-Datei geht, beschloß ich, mir einfach selbst eine Funktion zu schreiben, welche die Arbeit erledigt. Und hier ist sie ...
create or replace function convert_gpx_to_sdo(
 p_gpxfile in xmltype,
 p_lrs     in number default 0
) return sdo_geometry is
 v_ordinates sdo_ordinate_array := sdo_ordinate_array();
 v_gtype     number;

 v_cnt       pls_integer := 1;

 v_startts   timestamp;
 v_interval  interval day(9) to second;

 c_gpxns     varchar2(200) := 'xmlns="http://www.topografix.com/GPX/1/1"';
begin
 if p_lrs = 1 then
   v_gtype := 3302;
 else
   v_gtype := 2002;
 end if;

 for tp in (
   select
     extractvalue(value(tp), '/trkpt/@lon', c_gpxns) x,
     extractvalue(value(tp), '/trkpt/@lat', c_gpxns) y,
     extractvalue(value(tp), '/trkpt/time', c_gpxns) m
   from table(xmlsequence(extract(p_gpxfile, '//trkpt', c_gpxns))) tp
 ) loop
   if v_cnt = 1 then
     v_startts := to_timestamp(tp.m, 'YYYY-MM-DD"T"HH24:MI:SS"Z"');
   end if;
   if p_lrs = 1 then
     v_ordinates.extend(3);
   else
     v_ordinates.extend(2);
   end if;
   v_ordinates(v_cnt) := tp.x;
   v_cnt := v_cnt + 1;
   v_ordinates(v_cnt) := tp.y;
   v_cnt := v_cnt + 1;
   if p_lrs = 1 then
     v_interval := to_timestamp(tp.m,  'YYYY-MM-DD"T"HH24:MI:SS"Z"') - v_startts;
     v_ordinates(v_cnt) :=
      extract(MINUTE from v_interval) +
      extract(HOUR from v_interval) * 60 +
      extract(SECOND from v_interval) / 60;
     v_cnt := v_cnt + 1;
    end if;
 end loop;

 return sdo_geometry(v_gtype, 8307, null, sdo_elem_info_array(1,2,1), v_ordinates);
end convert_gpx_to_sdo;
/
sho err
Die Funktion nimmt das GPX-File (XMLTYPE) als ersten Parameter entgegen. Der zweite Parameter gibt an, ob eine Linear Referencing-Geometrie zurückgegeben werden soll. Wenn ja, werden die Measures in Minuten gespeichert. Gebt hier entweder "1" oder "0" an ... Man könnte hier sicherlich noch mehr machen: die GPX-Dateien enthalten eine Höhenangabe - man könnte also auch 3D-Geometrien bauen ...
Vielleicht jemand diese erste Version brauchen - viel Spaß damit!

Freitag, 12. November 2010

Views mit SDO_GEOMETRY

Auch für Geodaten-Tabellen werden recht häufig Views generiert. Daten aus verschiedenen Tabellen stehen damit quasi als "eine Tabelle" bereit. Und man kann mit den Views auch sehr gut arbeiten. Nur lassen sich die Views von vielen GIS-Werkzeugen und auch vom Oracle MapBuilder (das Werkzeug zum Einrichten der Kartendefinitionen für Oracle MAPS) nicht verwenden - die Auswahldialoge zeigen sie einfach nicht an.
Nun befindet sich auf der View zwar kein Spatial-Index, aber das ist ja auch richtig so; auf eine View kann kein Index gelegt werden. Der Grund ist die Spatial-Metadaten-View USER_SDO_GEOM_METADATA, die von all diesen Werkzeugen ausgelesen wird und meist Grundlage der Auswahllisten ist.
Die Lösung ist also einfach: Die View wird, wie die zugrundeliegende Spatial-Tabelle in USER_SDO_GEOM_METADATA eingetragen. Nehmt einfach dieses SQL hier als Vorlage.
insert into user_sdo_geom_metadata values (
  '{TABLE/VIEW NAME}',
  '{GEOM_COL_NAME}',
  sdo_dim_array(
   sdo_dim_element('X', {xmin}, {xmax}, {tolerance}),
   sdo_dim_element('Y', {ymin}, {ymax}, {tolerance})
  ),
  {SRID}
)
/
Einen Index könnt Ihr danach immer noch nicht anlegen - aber das braucht Ihr auch nicht - der View sollte nun in eurem GIS-Werkzeug oder im Oracle MapBuilder sicht- und auswählbar sein.