Hennings Dissblog

%% DENSITY SLICE   Scatterplot mit eingefärbter Punktdichte
% bild=densityslice(x,y,xp,yp);
%
% x und y müssen gleich groß sein, dürfen Vektoren oder Matrizen sein.
%
%   densityslice(x,y); Plot mit je 300 Pixeln in x- und y-Richtung
%
%   densityslice(x,y,xp,yp); Plot mit xp und xp Pixeln in x- und y-Richtung
%
% Output 'bild' kann mit image(bild) oder imagesc(bild) dargestellt werden.
% Wenn kein Output angegeben wird, wird der Plot nur dargestellt.

function bild=densityslice(x,y,xp,yp)

if nargin<4, yp = 300; end
if nargin<3, xp = 300; end

dat=[reshape(x,size(x,1)*size(x,2),1),reshape(y,size(y,1)*size(y,2),1)];

xmin = min(dat(:,1));
xmax = max(dat(:,1));
xstep= (xmax-xmin)/xp;
ymin = min(dat(:,2));
ymax = max(dat(:,2));
ystep= (ymax-ymin)/yp;

bild = zeros(xp,yp);

ii=xmin:xstep:xmax;
jj=ymin:ystep:ymax;
for xx=1:xp
    indexz = find(dat(:,1)>=ii(xx) & dat(:,1)<ii(xx)+xstep);
    datz = dat(indexz,:);
    for yy=1:yp;
        index=find(datz(:,2)>=jj(yy) & datz(:,2)<jj(yy)+ystep);
        bild(xx,yy)=length(index);
    end
end

%% Graphik
figure
image(bild); axis xy, colorbar

%% Achsenbeschriftung
XEP= [1; 1/3*xp; 2/3*xp; xp];
XE = [xmin; xmin+1/3*(xmax-xmin); xmin+2/3*(xmax-xmin); xmax];
set(gca,'xtick',XEP, 'xticklabel', num2str(XE,3));
YEP= [1; 1/3*yp; 2/3*yp; yp];
YE = [ymin; ymin+1/3*(ymax-ymin); ymin+2/3*(ymax-ymin); ymax];
set(gca,'ytick',YEP, 'yticklabel', num2str(YE,3));

% set(gca,'fontsize',14,'fontweight','bold')
%% Farbschema mit weißem Hintergrund
% j=jet;
% j(1,:)=[1 1 1];
% colormap (j)

Dieses Thema hat zwar nur wenig mit meiner Diss zu tun, aber ich mag es einfach gern, auch, weil es die schöne These "Science - it works, bitches" stützt. Ich versuche hier, die Eichung des Photometers ohne Formeln zu erklären, die Variante mit Formeln gibt es natürlich auch.

Mit dem Photometer können wir die Helligkeit der Sonne in verschiedenen Spektralkanälen, also verschiedenen Farben des Lichts messen (das weiße Sonnenlicht ist zusammengesetzt aus rotem, grünem, blauem und unsichtbarem infrarotem und ultraviolettem Licht). Was uns eigentlich interessiert, ist allerdings nicht die Ausstrahlung der Sonne, sondern die Trübung der Atmosphäre. Die Ausstrahlung wird als konstant angenommen, was natürlich nur näherungsweise stimmt.

Nun ist es sehr schwer, ein Meßgerät zu bauen, das in absoluten physikalischen Einheiten mißt, und das auch noch über Jahre. Unsere Photometer haben schon einige Jahre auf dem Buckel und geben nicht mehr das gleiche Ergebnis bei gleichen Einstrahlung wie in ihrer Jugend. Das macht aber nichts, denn wir können sie mit ihren eigenen Meßdaten eichen.

Wie gesagt interessiert uns die Trübung der Atmosphäre. Das Licht wird also beim Weg durch die Atmosphäre abgeschwächt. Und zwar bei einem längeren Weg stärker als bei einem kurzen. Die Idee hinter der Eichung ist daher, daß wir an einem sonnigen Tag mehrere Messungen machen: Bei niedrigem Sonnenstand und langem Weg durch die Atmosphäre morgens oder abends und bei hohem Sonnenstand und kurzem Weg mittags. Und zur Sicherheit noch ein paar dazwischen. Wenn man jetzt diese Messwerte auf der einen Achse und den Sonnenstand auf die andere Achse (beide etwas transformiert) aufträgt, erhält man eine Gerade. Und wenn man diese verlängert, kann man in der Graphik (die übrigens Langley-Plot genannt wird) ablesen, welchen Wert das Gerät ohne Atmosphäre gemessen hätte. Und schon hat man die Eichkonstante. Jetzt kann man für jeden Meßwert das Verhältnis des gemessenen Wertes zu dem Wert ohne Atmosphäre bilden und kennt die Transmissivität, also die Trübung der Atmosphäre.

Die Graphik zeigt Photometermessungen über einen Tag (sorry, ich wollte schon immer mal eine Graphik wie in einer Zeitschrift formatieren). Man sieht, daß die Einstrahlung morgens und abends schwächer ist als mittags. Das liegt am längeren Weg durch die Atmosphäre bei niedrig stehender Sonne. Nebenbei kann man auch noch sehen, warum der Himmel blau ist: Je kurzwelliger das Licht, umso stärker wird es von der Luft gestreut. Deshalb ist der Himmel voller blauem Licht.

Eine Langley-Plot-Graphik folgt noch.

Hier gibt es ein paar Daten.

path='H:\Laser\rohdaten-klassifiziert\holzbach-kacheln\';
file='holzbach_fe_25575_54920.asc';
data=textread([path file]);
index=find(data(:,1)>2557665 &...
data(:,1)<2557920 &...
data(:,2)>5492130 &...
data(:,2)<5492280);

d=data(index,:);
file2='holzbach_ground_25575_54920.asc';
data=textread([path file2]);
index=find(data(:,1)>2557665 &...
data(:,1)<2557920 &...
data(:,2)>5492130 &...
data(:,2)<5492280);
g=data(index,:);
plot3(d(:,1),d(:,2),d(:,3),'.','markersize',1)
hold on
plot3(g(:,1),g(:,2),g(:,3),'.k','markersize',1)
film=rotate(360);
film=orbit(360);
movie2avi(film,'kkh2.avi');

Um eine Convex Hull und damit ein Contiunuum Removal in einer Tabellenkalkulation zu berechnen, gibt es wahrscheinlich viele Möglichkeiten, einige mögen geschickter sein als die hier vorgestellte, aber sie funktioniert und ist recht gut steuerbar und damit flexibel.

Convex Hull ist eine Hüllkurve, die (in unserem Fall) ein Reflexionsspektrum umschließt. Man geht hier davon aus, daß das reine (Boden-)Spektrum eine konvexe (rechtsgekrümmte) Kurve, ein Kontinuum, ist, in das die Absorptionsbanden eingekerbt sind. Im Beispiel ist das Spektrum blau, die Hüllkurve violett. Das Continuum Removal teilt jetzt einfach das Spektrum bei jeder Wellenlänge durch die Hüllkurve, dadurch entsteht die grüne Kurve, die nur noch die normalisierten Absorptionsbanden enthält. Diese können nun unabhängig von einigen Randbedingungen wie Beleuchtung miteinander verglichen oder quantitativ ausgewertet werden (wie Peaks bei der Gaschromatograsphie).

Wie wurde das nun in Excel umgesetzt?

In den ersten beiden Spalten stehen Wellenlänge und Reflexionsgrad des betrachteten Spektrums, in diesem Fall Reaktorkorn aus einer Wasseraufbereitungsanlage, das nicht das erhoffte reine Kalkspektrum war, was aber für diese Betrachtung egal ist. In Spalte C werden Stützpunkte der Hüllkurve markiert. Das Verfahren ist also nicht automatisch, sondern erfordert einiges an Handarbeit, vor allem bei Laborspektren; in Excel mit 2151-zeiligen Datensätzen zu arbeiten ist kein Vergnügen. An den Stützpunkten – in der Graphik mit blauen Kreisen markiert – wird die Hüllkurve an das Spektrum "angeschmiegt". Hier haben beide den gleichen Wert, die kontinuumsbereinigte Kurve damit also den Wert 1. Optimalerweise hätten die Teile des Spektrums, die von alleine schon konvex sind, alle einen Stützpunkt, da hier beide Kurven aufeinander verlaufen sollen. Wenn irgendwo die R/H-Kurve einen Wert größer 1 hat, wurden zu wenige Stützpunkte gesetzt, hier also im Wellenlängenbereich um 700 nm. Macht aber nichts, weil meistens nur einzelne Absorptionsbanden betrachtet werden, und dafür spielt es nur eine Rolle, daß Anfangs- und Endstützpunkt richtig gesetzt sind.

In Spalte D werden nochmal die Reflexionswerte an den Stützpunkten wiederholt [=WENN(C2;B2;"")], das dient nur der graphischen Darstellung der blauen Kreise. Zwischen den Stützpunkten wird linear interpoliert. Hier muß man sich kurz mal an die Schulmathematik erinnern. Bei einer Gerade braucht man Steigung und y-Achsenabschnitt. Die Steigung m ergibt sich aus dem Steigungsdreieck bzw. dem Differenzenquotient Δy/Δx, der Achsenabschnitt an der Stelle x läßt sich aus 2 gegebenen Punkten (x1, y1) und (x2, y2) mit y1+m(x-x1) berechnen.

Dazu habe ich 4 Hilfsspalten eingebaut: E und F enthalten den x- und y-Wert des vorherigen Stützpunkts [=WENN($C2;A2;E1) und =WENN($C2;B2;F1)]. Es wird geschaut, ob die aktuelle C-Spalte einen Wert enthält (ein Stützpunkt ist). Wenn ja, wird der aktuelle x- oder y-Wert genommen, ansonsten derjenige, der letzte Zeile genommen wurde. Zeile 1 muß dann natürlich Stützpunkt sein. Entsprechend stehen in G und H die x- und y-Werte des nächsten Stützpunktes [=WENN(C3;A3;G3)...]. Hierfür habe ich noch fieser getrickst. Es wird geschaut, ob die nächste Zeile ein Stützpunkt ist. Wenn ja, wird der x- oder y-Wert der nächsten Zeile genommen, ansonsten der Wert, der auch letztes mal (von unten gesehen) schon genommen wurde. Vielleicht fällt mir irgendwann ein Weg ein, das verständlich auszudrücken.

Mit Hilfe dieser vier Spalten (die man natürlich auch in einer einzigen großen Formel unterbringen könnte) ist es dann recht leicht, in I die Hüllkurve berechnen zu lassen: =F2+(H2-F2)/(G2-E2)*(A2-E2).

In J wird schließlich mit =B2/I2 die kontinuumsbereinigte Kurve berechnet, und wahlweise in K noch mit =1-J die umgekehrte Kurve, in der man dann einfach Tiefen und Flächen von Absorptionsbanden/Peaks berechnen kann.

sind ja etwas verwirrend, also eine kleine Zusammenfassung durch Fragen (nicht umfassend, aber fernerkundlich):

Wieviel Energie strahlt die Sonne ab?
- Strahlungsenergie (J)

Wieviel Energie strahlt die Sonne pro Sekunde;
wieviel Leistung hat die Sonnenstrahlung?
- Strahlungsfluß (W)

Wieviel Energie strahlt die Sonne in Richtung Erde ab?
- Strahlungsintensität (W/sr)

Wieviel Energie kommt auf einer Fläche an?
- Bestrahlungsstärke (W/m²)

Wieviel Energie strahlt/reflektiert eine Fläche in Richtung eines Sensors?
- Strahldichte (W/sr·m²)

Wieviel Energie kann der Sensor davon in einem Kanal messen?
- Spektrale Strahldichte (W/sr·m²·µm)

3 Ansichten des Frankelbach-Höhenmodells

Zwischen dem Punktdichteverhältnis (First/Only gegen Ground) und dem Gelände-Kronenschlußgrad besteht ein relativ schwacher Zusammenhang. Aber wenn die Betrachtung auf Fichtenbestände eingeschränkt wird, gibt es immerhin eine Korrelation von 0.80.

%% Konfusionsmatrix
anz_klassen = 6;
format compact;

% ordner='F:\Bilddaten\Laser\Morbach\grids\5m-produkte\';
% file1='Klassi_ML_19MNF';
% [datak,datakx,dataky]=envilesengeo(ordner, file1);
% 
% file2='ROI_Vali4_ClassImage';
% ordner2='E:\Publikationen\Diplomarbeit\SAM\';
% [datav,datavx,datavy]=envilesengeo(ordner2, file2);
% 
% subplot(1,2,1),imagesc(datak),subplot(1,2,2),imagesc(datav)

matrix = zeros(anz_klassen,anz_klassen);
pmatrix=matrix;

for ii= 1:anz_klassen
    
    for jj= 1:anz_klassen
        
        index=find(datav==jj);
        index2=find(datak(index)==ii);
        matrix(ii,jj)=length(index2);
    end
end
disp ('Konfusionsmatrix nach Pixeln (mit Summen):')      
disp ([matrix, sum(matrix,2); sum(matrix,1) sum(sum(matrix))])


for ii= 1:anz_klassen
    for jj= 1:anz_klassen
        pmatrix(ii,jj)=matrix(ii,jj)./sum(matrix(:,jj));
    end
    
end
disp (' ')
disp ('Prozentuale Konfusionsmatrix:')
disp (100*pmatrix)

%% Prod Acc, User Acc, Error of O, Error of C
acc = zeros(anz_klassen,4);
for ii = 1:6
    acc(ii,1) = matrix(ii,ii)/sum(matrix(:,ii),1);
    acc(ii,2) = matrix(ii,ii)/sum(matrix(ii,:),2);
end
acc(:,3) = 1 - acc(:,1);
acc(:,4) = 1 - acc(:,2);
disp (' ')
disp('Prod Acc, User Acc, Error of O, Error of C')
disp(100*acc)
    
%% Kappa, nach Richards S. 267
P = sum(sum(matrix));
Sxkk = sum(diag(matrix));
Sxkppk = sum(matrix,1)*sum(matrix,2);

disp (' ')
Kappa = (P*Sxkk-Sxkppk)/(P^2-Sxkppk)

%% OAA
OAA = Sxkk / P


clear P Sxkk Sxkppk ii jj

Matlab kriegt es nicht hin, mit

[a b]=meshgrid(pr1(1):1:pr2(1),pr1(2):-1:pr2(2));
H=griddata(data(:,1),data(:,2),data(:,3),a,b,'nearest');

ALS-Daten ordentlich zu gridden, auch nicht mit 'linear', 'cubic' oder 'v4'. Die beiden letzten ergeben überhaupt keine Bilder, die anderen siehe Erweiterten Eintrag.

Datenausschnitt war pr1=[2582500;5520000];pr2=[2583000;5519600];

Daß mit Hilfe von data=fread(fid,[spalten zeilen],'float')'; Envibilder gut einzulesen sind, ist ja altbekannt. Ich hab jetzt mal eine Routine geschrieben, um auch den Header einzulesen, also Zeilen- und Spaltenzahl, Datentyp (noch ungetestet) und Geoinformation zu lesen. Bisher nur für einkanalige Bilder.

%% Read Envi
% [data,datax,datay]=envilesengeo(folder, file)
% Kann noch nicht wirklich mit unterschiedlichen Datentypen umgehen,
% momentan nur float. Bilddatei darf keine Dateiendung haben.
% Wird wahrscheinlich scheitern, wenn die Leerzeichenanzahl im Header
% anders als im Testdatensatz ist.
% Georef steckt auch noch in den Kinderschuhen...
% datax ist Vektor mit X-Werten, datay mit Y-Werten
% Also ist Koordinate von [a,b] [datax(a), datay(b)];
% [punkt(1)-xk yk-punkt(2)]
function [data,datax,datay]=envilesengeo(folder, file)

%file='morbachground_55234';
%folder='F:\Bilddaten\Laser\Morbach\grids\1m\';
dtstring='float';
datax=0;datay=0; baender=1;

%% Header lesen
fid=fopen([folder file '.hdr']);
while (~feof(fid))
    zeile=fgetl(fid);
    if (length(zeile)>6)
        % Spalten
        %if (zeile(1:7)=='samples'), spalten=str2num(zeile(10:end)), end
        if (zeile(1:7)=='samples')
            teile=explode(zeile,'=');
            spalten=str2num(teile{2})
        end
        % Zeilen
        %if (zeile(1:5)=='lines'), zeilen=str2num(zeile(8:end)), end
        if (zeile(1:5)=='lines')
            teile=explode(zeile,'=');
            zeilen=str2num(teile{2})
        end
        % Datentyp, siehe Envihandbuch 3.4 Seite 863
        if (zeile(1:7)=='data ty')
            dt=str2num(zeile(12:end));
            if dt==1, dtstring='int8', end
            if dt==2, dtstring='int16', end
            if dt==3, dtstring='int32', end
            if dt==4, dtstring='float', end
            if dt==5, dtstring='double', end
            if dt==12, dtstring='uint16', end
        end
        % Bänder
        if (zeile(1:5)=='bands')
            teile=explode(zeile,'=');
            baender=str2num(teile{2})
        end        
        % Georef
        if (zeile(1:8)=='map info')
            teile=explode(zeile,','); % Aufteilung in Teilstrings,Funktion
            xs=str2num(teile{2}); % Pixelsize x
            ys=str2num(teile{3}); % Pixelsize y
            xk=str2num(teile{4}) % X-koord Links oben
            yk=str2num(teile{5}) % Y-Koord links oben
            datax=xk:xs:xk+(xs*(spalten-1)); % Vektor mit X-Werten
            datay=yk:-ys:yk-(ys*(zeilen-1)); % Vektor mit Y-Werten
        end
    end
end
fclose(fid);

%% Daten lesen
% falls Dateiendung hier eintragen
ext=[]; %ext='.dat';

fid=fopen([folder file ext]);
data=fread(fid,[spalten zeilen],dtstring)';
fclose(fid);

Kaum sind die Daten da ist auch schon die erste Rechnung fertig.

Rohdaten gekachelt; Kacheln gerastert; Raster von Firstonly subsettet (Spalten 21 bis 9148) und Differenz sieht unschön aus. Mit "0>(b1-b2)<50" sieht das schon besser aus, auch wenns etwas geschummelt ist. So sieht das Histogramm dazu aus, wenn die ersten 3 Einträge, also die Werte um Null rausgelöscht sind:

Über das lange Wochenende den Rechner mal wieder laufenzulassen hat sich gelohnt; der Matlab hat drei ganz nette Lidar-Produkte ausgespuckt. Die Laufzeit hab ich dummerweise nicht gemessen, aber leider war es saulangsam: Ca. 20 h für einen von 40 Streifen.

Was ich getan habe (schließlich hat die Hansa die Ground-Daten noch nicht geliefert):

• Einen Rohdaten-(FirstOnly)-Streifen eingelesen
• Ein Moving Window laufen lassen. Allerdings waren es ja keine Rasterdaten, deshalb auch kein normales MW, sondern mit einer Schrittweite von 5 m ein neues Raster durchlaufen. Für jeden Mittelpunkt quadratisch (Henman wollte ein rundes Fenster probieren, auch nicht schlecht) +/- 5 m Umgebungspunkte gesucht.
• Wenn mehr als 3 Punkte existieren, 3 Parameter abgeleitet:
• Standardabweichung (daher der Name "Homogenität")
• Minimum (Primitive Ground-Suche)
• Höhen (Max. des Fensters - Min des Fensters)
• Weitere, die mir jetzt spontan einfallen
• Anzahl Punkte
• Range (achnee, hab ich ja schon)
• Eine bessere Bodenklassifikation

Die Ergebnismatrizen hab ich mit save maxdiff55168 hoehen -v6 gespeichert und in Envi eingelesen. Headerlänge war 184 Bytes. Die Bilder mußten noch um die y-Achse gespiegelt werden (Flip 270° und transpose).

Die Unterschiede zwischen den (angeblichen) First- und Last-Datensätzen sind tatsächlich sehr gering. Auch bei der Intensität muß man ganz genau hinschauen um was zu entdecken. Das Bild zeigt First-Last-First als RGB, Mosaik in 1 m Auflösung (Mosaikieren geht übrigens sehr schnell wenn man erst alle Kacheln öffnet).

Wenn man nur einen klein genugen Ausschnitt nimmt, kann man sogar ohne Überhöhung gute 3D-Ansichten in Envi bauen. Was mir zuerst nicht ganz klar war: Es wird immer ein Bild drübergelegt, und zwar das, welches gerade im Viewer offen ist. Wenn das Höhenmodell selber gewählt ist, sieht das ganze meist nach schneebedeckten Gipfeln aus. Animationen lassen sich auch recht einfach erstellen, aber für MPEGF-Export gibts hier keine Lizenz. Vielleicht versuch ich mal, das sonstwie zu exportieren.

Mist, da teilt man die Daten schon in 40 Schnipsel auf, und trotzdem.

Ich wollte von meinem IDL-Programm aus Morbach-Firstonly-Daten ein 0.5-m-Intensity-Raster rechnen lassen, und er hat auch recht vielversprechend angefangen. Aber bei Streifen 55177 war dann Schluß. Merken: 250 MB Inputdaten gehen, 263 MB nicht mehr.

Die Dimension des 1-m-Rasters 55177 ist 5526x300, so wären es also 11052*600, bei float also gute 24 MB. Das Speicherproblem kommt aber wohl nicht vom Raster sondern vom TIN (Triangulate in IDL uses double precision for all computations).

Delaunay in der Wikipedia

So sieht das Envi-Mosaik-Tool aus, wenn es mit 40 Scheiben gefüttert wird:

Wär ja zu einfach, wenn man die gelieferten Daten einfach so gebrauchen könnte (wäre es tatsächlich, wo soll dann der Stoff für eine Diss herkommen?).

Die Rasterung und Mosaikierung mit dem IDL-Skript hat zwar hervorragend geklappt (Laufzeit kann ich nicht genau sagen, weil ich es in Teilabschnitten hab rechnen lassen, aber wohl unter 24 Stunden). Aber leider muß ja eigentlich der First-Datensatz zuerst mit dem Only-Datensatz verknüpft werden, bevor ich anfange zu rastern. (Note to self: Illustration dazu hier einfügen).

Die roten Punkte stellen First-, die gelben Only-Pulses dar. Die alleinige Verwendung der First Pulses würde ein Oberflächenmodell wie durch die rote Linie dargestellt ergeben. Nur durch Hinzunahme der Only Pulses wird garantiert daß auch Lücken zwischen Bäumen im Oberflächenmodell erfaßt werden.

Dafür hab ich also mein altbewährtes Matlabskript modifiziert, das sich momentan der 48. Rechenstunde nähert. Input sind die für Morbach in 21, für Otterbach in 14 Segmente geschnittenen Rohdaten. Der Output (40 Kacheln/Streifen) ist 7,24 GB groß. Wahrscheinlich ist es sogar ratsam, das nochmal zu wiederholen, und dabei kleine Überlappungsbereiche einzubauen, damit das Mosaik keine schwarzen Streifen zeigt.

Funktionsweise ist ganz einfach: Aus Min und Max-y-Werten (x-Teiler ist 1, kann aber auch gesetzt werden) und der Anzahl der Teilstücke werden Dateien erzeugt, die jeweils die Datenzeilen aus First und Only enthalten, die zwischen ihrer Nord- und Südgrenze liegen. Zeile für Zeile wird eingelesen, zerteilt und in die entsprechende Outputdatei gesteckt. Erst First dann Only.

Heute kamen endlich die Lidar-Daten (wann war doch gleich die Befliegung? September?)

Sie sind eingeteilt in first, last und only, und die jeweils wieder in 21 Teile. Und für Otterbach nochmal das gleiche (ich weiß noch gar nicht ob das auch 21 sind).

Die alle in Envi zu konvertieren wäre natürlich sehr nervig, also hab ich mal geschaut obs auch in IDL geht. Es scheinen nur 2 Befehle notwendig zu sein: TRIANGULATE und TRIGRID. Auch das Einlesen ging recht gut über read_ascii. Mit ASCII_TEMPLATE() kann in einem GUI das Format des Ascii-Texts eingestellt werden, dies wird mit save und restore gesichert.

Dann blieb (fast) nur noch, das Bild auch zu speichern, und einen Header dazuzuschreiben. Speichern geht über WRITEU, geklammert von OPENW und FREE_LUN. Den Header hab ich mit PRINTF geschrieben. Eine kleine Klippe war noch, daß die erzeugten Bilder auf dem Kopf standen (Nullpunkt unten links statt oben links), jetzt werden sie vor dem Speichern mit rotate(raster,7) (7 steht für Drehung um 270° und transpose) umgedreht.