Posts mit dem Label gdal werden angezeigt. Alle Posts anzeigen
Posts mit dem Label gdal werden angezeigt. Alle Posts anzeigen

Freitag, 24. Oktober 2014

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.

Freitag, 23. Mai 2014

Ein GeoRaster ETL Werkzeug als grafische Oberfläche für gdal zum Laden nach SDO_GEORASTER

In meinem letzten Blog-Posting, dem 1. zum Thema Rasterdaten, hatte ich den Georaster-Viewer vorgestellt.
Für den Einstieg ins Thema und den schnellen Erfolg, ein Rasterbild in die Datenbank zu laden, war das sicher erst mal hilfreich.
Was das Laden oder Exportieren von Rasterdaten betrifft, kommt man sicher an gdal nur schwer vorbei. Daher wird die Nutzung auch schon länger von Oracle empfohlen. Allerdings braucht es doch ein wenig Übung, auf Kommandozeilenebene damit umzugehen. (Infos u.a. hier)

Eine grafische Oberfläche vom Spatial Development Team schafft Abhilfe. Diese verwendet intern gdal. Die Rede ist vom Georaster Batch ETL Werkzeug (GDAL-based GeoRaster Concurrent Batch Loading and Exporting Wizard). Dessen neueste Version kann von den Oracle Spatial and Graph Sample Code Webseiten im OTN heruntergeladen werden. Es kann sowohl für die Version 12c als auch 11g verwendet werden. Ich habe einen ersten Test mit Version 12cR1 durchgeführt.

Nach dem Download, Auspacken und Anpassen des Shell- bzw. Batch-Skripts (startGeoRasterETL.sh bzw. startGeoRasterETL.bat) kann es auch schon losgehen. Für den Einstieg hilft zudem noch ein kleines Handbuch.
Wie der Name nahelegen soll, können mehrere Rasterdaten-Files gleichzeitig geladen werden. Diese werden ausgeählt, mit den entsprechenden Ladeparameter-Werten versehen (u.a. Blocking oder Anzahl der Pyramid Levels).
Alles zusammen wird als XML-Parameterdatei gespeichert und anschließend als Job ausgeführt.
Die vier in die Tabelle CITY_IMAGES geladenen Rasterdaten-Dateien überprüfe ich dann mit dem schon bekannten GeoRaster-Viewer. Sie wurden mit den (automatisch vergebenen) RasterIDs 49−52 geladen.
Was es mit Blocking und Pyramiden auf sich hat, darüber soll es dann in einem neuen Blog-Posting gehen.