www.martinwerner.de
Motivation
Die Übungsaufgabe, in der eine Menge aus Polygonen A mit Labels versehen werden soll, indem eine andere Menge gelabelter Polygone B mit diesen Polygonen geschnitten wird und als Label des Polygons das Label genommen wird, mit dem die Polygone die meiste Fläche haben, ist in Geometrie sehr aufwändig zu lösen, zum Beispiel
Für jedes Polygon in A
Finde zu einem Polygon aus A alle scheidenden Polygone aus B
www.martinwerner.de
A litte inaccuracy saves a lot of calculation…
Eine einfachere Lösung des Problems wird geliefert, indem man den komplexen Schnitt zweier Polygone durch eine geeignete Menge an Punkt-in-Polygon-Tests ersetzt.
Dazu benötigen wir eine Menge an Punkten, die Flächeninhalt mit einer gewählten Genauigkeit repräsentieren kann.
Eine Möglichkeit besteht in einem regulären Gitter: Setze Punkte mit festen Abständen (z.B. je 1m in X und Y-Richtung), sodass jeder Punkt einen Quadratmeter repräsentiert. Dann kann man zuerst dieses Gitter aus Datensatz B labeln (Punkt in Polygon-Test) und
www.martinwerner.de
Gitter
Ein n-dimensionales reguläres Gitter ist eine Menge an Punkten,
deren Abstände in jeder Koordinate fest sind. Diese Abstände heißen Gitterabstände.
Ein georeferenziertes Gitter ist ein Gitter mit Koordinaten in einem Referenzsystem (Koordinatensystem)
www.martinwerner.de
Gitter Labeln
Werte, die einem regulären Gitter zugeordnet werden, sind in natürlicher Weise als Matrix aufzufassen:
j
i ⋯
⋮ ⋱ ⋮
⋯
www.martinwerner.de
Gitter Labeln
Werden bei festem Gitter jedem Gitterpunkt mehrere Werte zugeordnet, spricht man von Bändern.
www.martinwerner.de
Beispiel: Wellenlängen in optischen Bildern
+ +
www.martinwerner.de
Zur Übungsaufgabe
1. Wähle Gitterabstand.
2. Berechne Gitter / Raster
3. Fülle Raster-Layer mit Klassen aus Datensatz B (z.B. int8) 4. Fülle Raster-Layer mit ID aus Datensatz B (z.B. int64)
5. Reduziere auf ID und zähle
Bem.: Das ist durch ein Ablaufen des Rasters einfach möglich.
www.martinwerner.de
Raster- und Vektorkarten
Vektorkarten
Liste von geometrischen Primitiven (Line, …)
+ beliebig genau
+ verlustfrei skalierbar
+ Datengröße abhängig von Komplexität - Komplexe Flächen- und
Kollisionsanfragen - Spatiale Organisation
- Numerische Probleme (Artefakte)
Rasterkarten
• Information pro Flächenelement / Gitterpunkt
+ konstante Genauigkeit
+ Flächen- und Kollisionsanfragen einfach
+ Systeme können für Bildbearbeitung optimiert sein (MMX, SSE, CUDA etc.)
~ Datengröße abhängig von Fläche - keine verlustfreie Skalierbarkeit - weniger Artefakte (Anti-Aliasing)
- sehr viele Auflösungen für „schöne“ Anzeige benötigt
www.martinwerner.de
Struktur einer Vektorkarte
www.martinwerner.de
Struktur einer Rasterkarte
www.martinwerner.de
Geometrie in Rasterkarte (horizontal, vertikal)
www.martinwerner.de
Geometrie in Rasterkarte (Diagonal)
www.martinwerner.de
Schneiden sich diese zwei Linien?
www.martinwerner.de
Pixel werden
entsprechend ihres
Flächenanteils gefüllt oder gefärbt.
Für Farben einfach (weil vom menschlichen Auge korrekt wahrgenommen).
Für z.B. Raster mit
diskreten Werten (unsere Klassen von vorhin) ist das mitunter etwas schwierig Lösung: Sub-Pixeling / Anti-Aliasing
www.martinwerner.de
Raster-Darstellung in der echten Welt
gdalinfo germany-cloudfree-simple.tiff Driver: GTiff/GeoTIFF
Files: germany-cloudfree-simple.tiff Size is 67235, 94080
Coordinate System is:
LOCAL_CS["WGS 84 / Pseudo-Mercator", GEOGCS["WGS 84",
DATUM["unknown",
SPHEROID["unretrievable - using WGS84",6378137,298.257223563]], PRIMEM["Greenwich",0],
Angaben zu - Format
- Anzahl Pixel
- Koordinatensystem (WKT)
www.martinwerner.de
Raster-Darstellung in der echten Welt
Origin = (641594.910993173718452,7370990.488857678137720) Pixel Size = (15.545613327413848,-15.545613327413848)
Metadata:
AREA_OR_POINT=Area Image Structure Metadata:
INTERLEAVE=PIXEL Corner Coordinates:
Upper Left ( 641594.911, 7370990.489) Lower Left ( 641594.911, 5908459.187) Upper Right ( 1686804.223, 7370990.489) Lower Right ( 1686804.223, 5908459.187) Center ( 1164199.567, 6639724.838)
Band 1 Block=67235x1 Type=Byte, ColorInterp=Red Band 2 Block=67235x1 Type=Byte, ColorInterp=Green
- Ursprung in Geokoordinaten - Größe der Pixel - Fläche oder Punkt - Strukturinformation - Geokoordinaten
(berechnet) - Bänder
www.martinwerner.de
Bemerkungen
Rastergeometrie im Detail:
Ein GeoTIFF-Raster besteht aus einer rechteckigen Matrix mit Werten pro Band. Eine Geotransformation besagt, wie sich mit einem Schritt nach rechts (bzw. nach unten) in dieser Matrix die Weltkoordinaten ändern. Diese Veränderung ist konstant.
Elt. Bedeutung
GT[0] X-Koordinate in CRS oben links (für Pixel 0,0)
GT[1] X-Offset in X-Richtung für jeden Pixel-Schritt nach Rechts
www.martinwerner.de
Berechne Weltkoordinaten für beliebigen Pixel
Mit der Geotransformation bestimmen die folgenden Gleichungen die Weltkoordinaten eines Pixels mit (ganzen) Koordinaten i,j:
= + ⋅ + ⋅
= + ⋅ + ⋅
www.martinwerner.de
Einfache Rasteroperationen
gdal_info:
Zeige Raster-Informationen an gdal_warp:
Zusammensetzen von Teilrastern
Ändern des Gitters / der Projektion (Warping) gdal_translate
Ändere Format
Subsetting / Resampling / Interpolation gdal_rasterize
Übersetze (OGR-)Geometrie in Raster
www.martinwerner.de
Rasteralgorithmen
www.martinwerner.de
Rasterisierung
Übersetzen von
Geometrieobjekten in Raster
Beispiel: Punkte (mit Radius)
Berechne Pixelkoordinaten in Fließkomma
Berechne anteilige Fläche pro Pixel
Modifiziere Pixel (Blending
www.martinwerner.de
Rasterisierung: Linien
Einfachste Möglichkeit: Interpoliere und Raster Punkte
= + ⋅
Problem: Welche Schrittweite? Tradeoff zwischen Aufwand und Wahrscheinlichkeit, alle Pixel auch zu treffen.
www.martinwerner.de
Bresenham-Algorithmus
Beobachtung: Nach Setzen einen Pixels kann der nächste Pixel nur ein Nachbarpixel sein
Einschränkung auf Oktant: Linien, die nach unten-Rechts gehen
Nur Linien in
www.martinwerner.de
Welcher Y-Schritt
Allgemeine Liniengleichung
−
− = −
−
Pro ganzzahligem Schritt in X-Richtung berechnet man dann
= −
− ⋅ − +
/
/ /
www.martinwerner.de
Optimierung
Der Bresenham-Algorithmus berechnet jetzt nicht die echte
Y-Koordinate, sondern verfolgt die ganzzahlige Y-Koordinate und den Fehlerterm (also die Differenz des echten Y-Wertes und des aktuellen ganzzahligen Y-Wertes. Sobald dieser Fehler größer wird als 1, wird er um 1 verringert und die Y-Koordinate angepasst.
Explizit also Pseudo-Code also:
error=0, int y = y0, float deltaerr = abs( (y1-y0)/(x1-x0)) for x from x0 to x1
putpixel(x,y)
www.martinwerner.de
Beispiel: Bresenham-Algorithmus für Linien (2)
Eine Fallunterscheidung unter Ausnutzung der Symmetrie wird jetzt verwendet, um den Basisalgorithmus auf alle 8 Oktanten abzubilden.
Merke: Der Bresenham-Algorithmus hat immer eine schnelle Richtung und eine langsame Richtung. In die schnelle Richtung (X oder Y)
macht man immer einen Schritt, in die langsame Richtung nur ab und zu.
Übungsaufgabe: Implementieren Sie Bresenham‘s Algorithmus für eine Linie in Python. Verwenden Sie zur Eingabe Fließkommazahlen im
Bereich zwischen 0 und 100 und als Raster eine Numpy-Matrix der