Posts mit dem Label GeoRasterViewer werden angezeigt. Alle Posts anzeigen
Posts mit dem Label GeoRasterViewer 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.

Donnerstag, 15. Mai 2014

In wenigen Schritten zu Rasterdaten in der Oracle Datenbank

Der Management von Vektordaten (SDO_GEOMETRY) ist in diesem Blog schon recht ausführlich behandelt worden. Daher möchte ich in meinem Beitrag heute zeigen, wie auf einfache Art und Weise Rasterdaten als SDO_GEORASTER geladen werden können. Voraussetzung dafür ist, dass die Oracle Database Examples von OTN heruntergeladen und installiert werden.
Relevante Links für die Datenbank-Version 12c sind:
Was finden Sie nach der Installation vor?
Im Verzeichnis %ORACLE_HOME%\md wurden der Ordner demo und weitere themenspezifische Unterordner neu angelegt. So auch einer mit dem Namen georaster.

Für die weiteren Ausführungen nutze ich die SQL-Skripte im Ordner plsql als Basis sowie den GeoRasterViewer, der im Ordner java liegt und den ich schon mal starte.

Was wird noch benötigt? Natürlich ein paar Beispiel-Bilder. Ich nutze .tif-Dateien von San Francisco mit den zugehörigen World Files (.wld).

Von der Anzeige eines Rasterbildes aus der Oracle Datenbank heraus bin ich jetzt nur noch 3 Schritte entfernt:
  1. Anlegen einer SDO_GEORASTER-Tabelle und der zugehörigen Raster Data Table(s)
  2. Initialisieren des SDO_GEORASTER Objekts (oder gleich mehrerer)
  3. Laden des Rasterbildes über den GeoRasterViewer
Das sieht dann so aus:
-------------------------------------------------------------------
-- Schritt 1: Anlegen einer SDO_GEORASTER-Tabelle.
--            Anlegen der sogenannten Raster Data Table (RDT).
--
-- Tabelle muss eine Spalte vom Typ SDO_GEORASTER enthalten.
-------------------------------------------------------------------

drop table city_images_rdt_1 cascade constraints purge
/
drop table city_images  cascade constraints purge
/
create table city_images (
  id          number primary key,
  city        varchar2(100),      
  type        varchar2(32), 
  georaster  mdsys.sdo_georaster)
/
create table city_images_rdt_1 of mdsys.sdo_raster (
  primary key (rasterId, pyramidLevel, bandBlockNumber, rowBlockNumber, columnBlockNumber))
  lob(rasterblock) store as securefiles (nocache nologging)
/

-------------------------------------------------------------------
-- Schritt 2: Initialisieren des SDO_GEORASTER Objekts.
-------------------------------------------------------------------

insert into city_images values( 
  1, 
  'San Francisco', 
  'TIFF', 
  mdsys.sdo_geor.init('CITY_IMAGES_RDT_1'))   -- RasterID wird automatisch generiert 
/                                             -- bei fehlendem 2. Parameterwert
commit
/
Für Schritt 3 wird im GeoRasterViewer eine Verbindung auf die Datenbank geöffnet mit dem Nutzer, welcher Eigentümer der zuvor angelegten Tabellen ist.
Zunächst erscheint beim Auswählen des initialisierten SDO_GEORASTER-Objekts (es hat die automatisch generierte ID 29) eine Fehlermeldung:
Das ist richtig so. Habe ich doch das Rasterbild selbst noch gar nicht geladen.
Aber ich hole das jetzt nach mit Tools > Import into DB. Bei den Ladeoptionen entscheide ich mich dafür, kein Blocking anzuwenden.
Die Bestätigung, dass das Rasterbild geladen wurde, erfolgt prompt.
Und das Ergebnis selbst sieht dann so aus:
Ggf. müssen die Daten noch mal gelesen werden über Rasters > Refresh List, damit das Rasterbild wie erwartet angezeigt wird. Aber das war es dann auch schon.
Ausführliche Hinweise finden sich wie immer in der Online Dokumentation im Oracle Technology Network. Für die Rasterdatenverwaltung mit Oracle Spatial gibt es ein eigenes Dokument (hier der Link auf die Dokumentation für Version 12.1).

Zusätzliche Hinweise:
  • Eine Besonderheit ist bei den Beispiel-Dateien ist zu beachten. ESRI WorldFiles enthalten keine Information &uunml;ber das Bezugssystem. Daher muss die SRID beim oder nach dem Laden manuell gesetzt werden. Dazu wird die Prozedur sdo_geor.setModelSRID verwendet.
  • Wie von den Vektordaten mit SDO_GEOMETRY bekannt, werden auch für SDO_GEORASTER Metadaten registriert (also nach dem Anlegen der Tabelle in Schritt 1).
  • Der in früheren Oracle Spatial Versionen (ich habe mit 12.1 getestet) manuell für jede Tabelle mit SDO_GEORASTER anzulegende DML-Trigger (sdo_geor_utl.createDMLTrigger), wird automtisch erzeugt. Dieser trägt den Präfix GRDMLTR.
  • "Tonnenweise" Rasterdaten in Oracle Spatial gibt es übrigens im Geoproxy des Freistaates Thüringen.
-------------------------------------------------------------------
-- Prozedur zum Setzen der SRID
-------------------------------------------------------------------
declare
  geor  sdo_georaster;
begin
  select georaster into geor from city_images where id = 1 for update;
  sdo_geor.setModelSRID(geor, 26943);
  update city_images set georaster = geor where id = 1;
  commit;
end;
/

-------------------------------------------------------------------
-- Metadaten prüfen
-------------------------------------------------------------------
select * from user_sdo_geor_sysdata
/

-------------------------------------------------------------------
-- SDO_GEORASTER Objekt validieren
--
-- Jeweils erwartetes Ergebnis: isvalid = TRUE
-------------------------------------------------------------------
select 
  r.id,
  mdsys.sdo_geor.validategeoraster(r.georaster) isvalid
from 
  city_images r 
order by 
  id
/

select 
  r.id,
  mdsys.sdo_geor.schemavalidate(r.georaster) isvalid
from 
  city_images r 
order by
  id
/