Vegetationsindex berechnen mit Oracle Spatial und Raster Algebra in der Datenbank

Seit Oracle 12c können direkt in der Datenbank Raster Algebra Operationen ausgeführt werden.
Damit ist es nun möglich, z.B. Vegetationsindizes wie NDVI für Rasterdaten zu berechnen, die als SDO_GEORASTER in der Oracle Datenbank gespeichert sind.

Hierzu sind im Wesentlichen nur 3 Schritte notwendig:

1. Laden der Rasterdaten in die Datenbank
2. Ausführen der entsprechenden Raster Algebra Operation
3. Anzeige des Ergebnisses

Die 3 Schritte beschreibe ich nachfolgend ein wenig genauer.

Wer das Ganze Schritt für Schritt ausprobieren möchte, ohne die Oracle Datenbank selber installieren zu müssen, die/der kann einen Workshop mit vorkonfigurierter virtueller Maschine und Testdaten dafür nutzen.

Schritt 1: Laden der Rasterdaten

Das Laden habe ich bereits in diesem Blog Posting beschrieben. Alternativ kann man GDAL nutzen und mit Hilfe eines Kommandozeilenbefehls ein Rasterbild direkt in die Datenbank schreiben. Beispielhaft sieht der Aufruf so aus:

gdal_translate -of georaster LT50440332011261PAC01.tif \
	georaster:student/student@$ORACLE_SID,napa_landsat,image \
	-co description="(name varchar2(32), description varchar2(64), image sdo_georaster, primary key (name))" \
	-co insert="('original', 'Landsat 7 bands 2011 Ago 28', sdo_geor.init('napa_landsat_rdt$'))" \
	-co genpyramid=NN \
	-co blockbsize=1

GDAL_TRANSLATE kann optional die Tabelle NAPA_LANDSAT mit IMAGE als SDO_GEORASTER-Spalte im DB Schema STUDENT gleich miterzeugen, falls sie vorher noch nicht angelegt war.

Mit Hilfe des Werkzeuges GeoRasterViewer kann das Ergebnis des Ladevorgangs gleich überprüft werden.

Schritt 2: NDVI direkt in der Datenbank berechnen

Die dem NDVI zugrunde liegende Formel lautet:

NDVI = (NIR-RED)/(NIR+RED)

NIR entspricht dabei dem zurückgestrahlten Licht im Nahen Infrarotbereich, RED dem roten Anteil des sichtbaren Lichts. Bei Landsat-Bildern sind das die Kanäle 3 und 4.

Grundlage für die Berechnung ist das Unterprogramm RASTERMATHOP im PL/SQL Package SDO_GEOR_RA. Dieses wird genutzt, um den NDVI über das zuvor geladene Rasterbild zu rechnen.

set echo on
set serverout on

declare
  gr1 sdo_georaster;
  gr2 sdo_georaster;
  opr sdo_string2_array;
  sts sdo_number_array;
begin
  select image into gr1 from napa_landsat where name = 'original';

-- Tabelle NAPA_IMAGES wurde zuvor angelegt 

  delete from napa_images where name = 'ndvi';
  insert into napa_images 
    values ('ndvi', 'Full scene NDVI', 
            sdo_geor.init('napa_images_rdt$'))
    return image into gr2;

--
-- Raster Algebra Operation - NDVI, Normalized Difference Vegetation Index
--
-- String Array mit Formel für NDVI
    
  opr := sdo_string2_array('(({3}-{2})/({3}+{2}+0.000001))'); -- Nummerierung der Kanäle startet bei 0
  
  sdo_geor_ra.rasterMathOp(inGeoRaster  =&gr; gr1, 
                           operation    =&gr; opr, 
                           storageParam =&gr; 'celldepth=32bit_real', 
                           outGeoRaster =&gr; gr2, 
                           nodata       =&gr; 'TRUE');

-- Zusätzlich Statistiken berechnen
  sts := sdo_geor.generateStatistics(georaster      =&gr; gr2,
                                     pyramidLevel   =&gr; 0,
                                     samplingFactor =&gr; 'samplingFactor=1',
                                     samplingWindow =&gr; null,
                                     bandNumbers    =&gr; null,
                                     nodata         =&gr; 'TRUE',
                                     polygonClip    =&gr; 'FALSE');

-- Statistiken ausgeben
  dbms_output.put_line('Min    = '||sts(1));
  dbms_output.put_line('Max    = '||sts(2));
  dbms_output.put_line('Mean   = '||sts(3));
  dbms_output.put_line('Median = '||sts(4));
  dbms_output.put_line('Mode   = '||sts(5));
  dbms_output.put_line('StdDev = '||sts(6));

-- Bildpyramiden erzeugen
  sdo_geor.generatePyramid(gr2, 'resampling=NN');
  
  sdo_geor.setID(gr2, 'ndvi');
  
-- Rasterbild mit NDVI in die Datenbank schreiben
  update napa_images set image = gr2 where name = 'ndvi';
  commit;
end;
/

exit

Die mit SDO_GEOR.GENERATESTATISTICS berechneten Statistiken haben folgende Werte ergeben:

Min    = -.936254978179932
Max    = .982608675956726
Mean   = .297637717413117
Median = -.936240315437317
Mode   = -.936240315437317
StdDev = .206702946761902

Schritt 3: Ergebnis visualisieren

Der letzte Schritt ist jetzt einfach und nur eine Wiederholung dessen, was schon gemacht wurde. Diesmal nutze ich aber uDig (User-Friendly Desktop GIS) anstelle des GeoRasterViewer für die Anzeige.

uDig unterstützt den direkten Zugriff auf die Oracle Datenbank und kann sowohl mit Vektor- (SDO_GEOMETRY) als auch Rasterdaten (SDO_GEORASTER) arbeiten.

Probiert es mal aus. Nutzt den Workshop "Developing Geospatial Applications With uDig and Oracle Spatial". Einfacher geht es nicht.