Rasterio: Numpy for Geospatial Data

Hallo, hier ist Alex von DEV STORIES DOT EU. Rasterio: Numpy für Geodaten, Folge eins von fünf. Du willst Pixelwerte aus einem Satellitenbild extrahieren, aber Standard-Image-Libraries stoßen bei Geodatenformaten an ihre Grenzen, und die offiziellen Bindings fühlen sich an, als würdest du C-Code in Python schreiben. Genau deshalb gibt es Rasterio. Rasterio ist eine Library zum Lesen und Schreiben von rasterbasierten Geodaten wie GeoTIFF-Dateien. Du hast wahrscheinlich schon von GDAL gehört, der branchenüblichen Geodaten-Library. Ein weit verbreiteter Irrtum ist, dass Rasterio GDAL ersetzt. Das tut es nicht. Rasterio nutzt GDAL unter der Haube. Der Unterschied ist das Interface. Wenn du jemals die nativen GDAL Python Bindings verwendet hast, weißt du, dass sie eine umständliche C-API mit manuellem Resource Management und einer für Python völlig ungewohnten Syntax bieten. Rasterio lässt das hinter sich. Es bietet ein modernes, Python-first Interface. Es behandelt Raster-Datasets wie normale Dateien und ihre Pixel wie Standard-Numpy-Arrays. Angenommen, du hast ein GeoTIFF mit Landsat-Bildern auf deiner Festplatte. Du öffnest diese Datei mit der rasterio dot open Funktion. Da Rasterio Standard-Python-Patterns folgt, machst du das mit einem with statement. Das erstellt einen Context Manager. Genau wie beim Öffnen einer einfachen Textdatei stellt der with block sicher, dass das Dataset sauber geschlossen und der Speicher freigegeben wird, sobald dein Code durchgelaufen ist. Du musst Objekte nie manuell zerstören oder dir Sorgen um Memory Leaks durch nicht geschlossene Dateien machen. Innerhalb dieses with blocks gibt dir die open Funktion ein Dataset Reader Objekt. Bevor du die eigentlichen Pixel liest, kannst du die Datei inspizieren. Du kannst das Dataset nach seinem count fragen, was dir die Gesamtzahl der Bands, oder Layer, im Bild verrät. Bei einem Landsat-Bild könnte dieser count zum Beispiel sieben oder elf sein, was verschiedene Lichtwellenlängen repräsentiert. Du kannst die width und height Properties auslesen, um die räumlichen Dimensionen in Pixeln zu bekommen. Du kannst auch die dtypes Property checken. Das gibt dir den Data Type der Pixel für jedes Band, wie zum Beispiel unsigned 16-bit Integers oder 32-bit Floats. All diese Metadaten sind sofort verfügbar, ohne dass die schweren Bilddaten in den Memory geladen werden müssen. Hier ist der entscheidende Punkt. Wenn du bereit bist, dir die eigentlichen Pixel anzusehen, rufst du die read Methode auf dem Dataset Objekt auf. Du kannst einen bestimmten Index übergeben, um nur einen Layer zu laden. Pass hier gut auf: Rasterio Bands sind one-indexed, das heißt, du fragst nach Band eins, nicht Band null. Wenn du read aufrufst und die Zahl eins übergibst, gibt die Methode diese Pixel als standardmäßiges zweidimensionales Numpy-Array zurück. Wenn du read ohne Argumente aufrufst, liest es alles und gibt ein dreidimensionales Numpy-Array zurück, das alle Bands enthält. Es gibt kein spezielles, proprietäres Pixel-Objekt. Du bekommst einfach ein Array aus Zahlen. Sobald du dieses Numpy-Array hast, bist du wieder in bekanntem Terrain. Du kannst es slicen, statistische Operationen ausführen oder es direkt in Machine Learning Models füttern, und zwar mit exakt demselben Code, den du für jede andere Matrix benutzt. Rasterio fungiert als Brücke. Es kümmert sich um die komplizierten File-Formate und Metadaten auf der Disk und übergibt dir ein sauberes mathematisches Objekt im Memory. Wenn dir die Show gefällt und du unsere Arbeit unterstützen willst, kannst du auf Patreon nach DevStoriesEU suchen – das ist eine riesige Hilfe. Die wahre Stärke von Rasterio ist nicht, dass es Geospatial Processing neu erfindet, sondern dass es das Ganze langweilig macht; indem es komplexe Satellitenbilder in Standard-Numpy-Arrays verwandelt, erlaubt es dir, nicht länger mit GIS-Formaten zu kämpfen und stattdessen echte Data Science zu betreiben. Danke fürs Zuhören, happy coding zusammen!

Hallo, hier ist Alex von DEV STORIES DOT EU. Rasterio: Numpy für Geodaten, Folge 2 von 5. Ein einfaches Numpy-Array ist einfach ein flaches Grid aus Zahlen. Wenn du es visualisierst, bekommst du einfach ein farbiges Rechteck auf deinem Monitor. Was verwandelt diese abstrakten Zahlen eigentlich in eine präzise Map des Colorado-Plateaus? Die Antwort ist die Affine Transform. Oft denkt man, Geodaten sind einfach nur ein Image, bei dem die Ecken grob an Map-Coordinates angeheftet sind. Das ist aber nicht der Fall. Die Transform ist ein strikter mathematischer Link, der die genaue Real-World-Position von jedem einzelnen Pixel in deinem Array berechnet. Bevor du Pixel auf die Erde mappen kannst, musst du wissen, wie du die Erde überhaupt misst. Das ist das Coordinate Reference System. Du liest das aus dem Dataset aus, indem du auf dataset dot crs zugreifst. Diese Property sagt Rasterio, ob deine räumlichen Koordinaten in Längen- und Breitengraden gemessen werden, oder in Metern relativ zu einem bestimmten Äquator und Nullmeridian. Sobald das Reference System steht, brauchst du die Mapping-Logik. Die findest du, indem du dataset dot transform checkst. Das liefert dir eine affine Transformationsmatrix. Das klingt nach harter linearer Algebra, aber es ist im Grunde einfach ein Set aus sechs Zahlen, die den physischen Footprint von deinem Grid beschreiben. Die Matrix definiert die exakten räumlichen x- und y-Koordinaten der oberen linken Ecke des allerersten Pixels oben links. Sie definiert die physische Breite von einem Pixel, zum Beispiel exakt dreißig Meter. Sie definiert die physische Höhe von einem Pixel. Die ist meistens negativ, weil Image-Rows nach unten zählen, während Map-Coordinates nach oben zählen. Die Matrix enthält außerdem Rotation-Values, für den Fall, dass der Satellit bei der Bildaufnahme schräg stand. Hier ist die wichtigste Erkenntnis: Du musst diese Matrix-Werte nie von Hand multiplizieren. Rasterio gibt dir direkte Methods, um zwischen Pixel-Space und Coordinate-Space hin und her zu springen. Angenommen, du analysierst ein Image und findest eine Anomalie an einer bestimmten Location, sagen wir Row fünfhundert und Column eintausend. Um herauszufinden, wo das auf dem echten Planeten ist, übergibst du diese beiden Integers an die dataset dot xy Method. Rasterio pusht die Row und Column durch die affine Matrix und gibt dir die exakten räumlichen x- und y-Koordinaten für den Mittelpunkt von genau diesem Pixel zurück. Du kannst diese Maschine auch rückwärts laufen lassen. Sagen wir mal, du hast die exakten räumlichen Koordinaten für eine physische Location. Vielleicht ist es ein Canyon, der genau einhundert Kilometer östlich und fünfzig Kilometer südlich von deinem Image-Origin liegt. Du musst die zugrundeliegenden Pixel-Daten für diesen Canyon extrahieren. Du nimmst deine räumlichen x- und y-Werte und übergibst sie an die dataset dot index Method. Rasterio invertiert die affine Matrix, verarbeitet deine Koordinaten und gibt dir die exakten Row- und Column-Integers zurück. Du nutzt diese Integers dann, um in dein Numpy-Array zu slicen und den exakten Messwert auszulesen, den der Satellit aufgezeichnet hat. Die Affine Transform isoliert dein Raw-Data-Array von deiner räumlichen Geometrie und garantiert so, dass du deine genaue Position auf dem Globus niemals verlierst, egal wie du deine Daten croppst oder liest. Das war's für heute. Danke fürs Zuhören – geh und bau was Cooles.

Hallo, hier ist Alex von DEV STORIES DOT EU. Rasterio: Numpy für Geospatial Data, Folge 3 von 5. Wenn du versuchst, ein zehn Gigabyte großes GeoTIFF direkt in den Arbeitsspeicher deines Laptops zu laden, stürzt dein Script mit ziemlicher Sicherheit ab. Du brauchst vielleicht nur einen winzigen Pixelbereich, um einen Algorithmus zu testen, aber das System versucht, den ganzen Kontinent auf einmal zu schlucken. Die Lösung für dieses Memory-Limit ist Windowed Processing. Mit Windowed Processing kannst du bestimmte rechteckige Subsets eines Raster-Files lesen oder schreiben, während der Rest des Files sicher auf deiner Festplatte bleibt. Das ist komplett indexbasiert. Wir beschäftigen uns heute nicht mit geografischen Koordinaten. Wir arbeiten ausschließlich mit Pixel-Rows und Pixel-Columns. Um ein Subset zu lesen, brauchst du einen Weg, Rasterio mitzuteilen, welche Pixel es sich greifen soll. Das machst du, indem du das Window-Objekt aus dem rasterio windows-Modul importierst. Der Standardweg, ein Window zu definieren, ist die Übergabe von vier Zahlen: Column-Offset, Row-Offset, Width und Height. Pass an dieser Stelle gut auf. Die Reihenfolge der Argumente ist wichtig. Es verlangt zuerst den Column-Offset, dann den Row-Offset. Der Column-Offset ist deine horizontale Startposition, gemessen vom linken Rand. Der Row-Offset ist deine vertikale Startposition, gemessen vom oberen Rand nach unten. Nach den Offsets gibst du an, wie viele Pixel breit und wie viele Pixel hoch deine Auswahl sein soll. Stell dir ein riesiges Test-File im Gigabyte-Bereich vor. Du möchtest eine kleine Kachel extrahieren, die 256 Pixel breit und 512 Pixel hoch ist. Diese Kachel soll eintausend Pixel vom linken Rand und zweitausend Pixel vom oberen Rand entfernt beginnen. Du instanziierst dein Window, indem du eintausend für den Column-Offset, zweitausend für den Row-Offset, 256 für die Width und 512 für die Height übergibst. Als Nächstes öffnest du dein Dataset mit der Standard-open-Funktion. Wenn du die read-Methode aufrufst, lässt du die Argumente nicht leer. Leere Argumente sagen Rasterio, dass es das gesamte File in den Memory laden soll. Stattdessen weist du dein neues Window-Objekt dem window-Keyword-Argument innerhalb der read-Methode zu. Rasterio übersetzt dieses Window in Byte-Offsets auf der Disk, holt sich nur diesen spezifischen 256 mal 512 Chunk und gibt ihn als Numpy-Array zurück. Dein Memory-Verbrauch bewegt sich kaum. Wenn du den ganzen Tag mit Numpy arbeitest, fühlt sich die Übergabe von Offsets und Dimensions vielleicht etwas unnatürlich an. Du bist es wahrscheinlich gewohnt, Slices mit Start- und Stop-Indizes für deine Rows und Columns zu definieren. Rasterio unterstützt diesen Workflow mit einer alternativen Methode namens Window dot from slices. Diese Methode nimmt zwei Inputs. Der erste Input ist ein Tuple, das den Row-Start-Index und Row-Stop-Index definiert. Der zweite Input ist ein Tuple, das den Column-Start-Index und Column-Stop-Index definiert. Beachte, dass die Reihenfolge hier zuerst Rows, dann Columns ist. Das spiegelt genau wider, wie du ein zweidimensionales Numpy-Array slicest. Unter der Haube berechnet diese Methode die Offsets und Dimensions für dich und gibt ein Standard-Window-Objekt zurück, das du genau wie vorher an die read-Methode übergeben kannst. Indem du Windowed Reads in deinen Data-Pipelines erzwingst, entkoppelst du die physische Größe deiner Input-Files dauerhaft von den Hardware-Limits deiner Maschine. Das war's für diese Folge. Danke fürs Zuhören, und keep building!

Hallo, hier ist Alex von DEV STORIES DOT EU. Rasterio: NumPy für Geodaten, Folge 4 von 5. Du hast gerade ein riesiges Satellitenbild heruntergeladen, das ein ganzes Bundesland abdeckt, aber du interessierst dich nur für die Pixel innerhalb einer bestimmten Farmgrenze. Das ganze Grid zu laden, verschwendet Arbeitsspeicher, und das Lesen über Pixel-Offsets ist nutzlos, wenn deine Grenze eine unregelmäßige Polygonform hat. Die Lösung ist Masking mit Vektoren. In einer vorherigen Folge haben wir uns Windowed Reading angesehen. Dieser Ansatz erforderte exakte Pixel-Indizes, so als würdest du nach Zeile fünfzig und Spalte einhundert fragen. Masking funktioniert völlig anders. Masking basiert auf echten geografischen Koordinaten. Anstatt Pixel zu zählen, übergibst du eine Vektor-Geometrie, wie ein Polygon, das durch echte Längen- und Breitengrade definiert ist. Rasterio findet heraus, welche Pixel in diese Form fallen, und übernimmt die komplexe Translation zwischen dem kontinuierlichen Vektor-Koordinatenraum und dem diskreten Raster-Grid. Das Kernwerkzeug für diesen Workflow ist die mask-Funktion, die sich im Modul rasterio dot mask befindet. Du übergibst ihr zwei Haupt-Argumente. Erstens, ein offenes Raster-Dataset. Zweitens, ein Iterable von geometrischen Formen. Diese Formen müssen als Dictionaries formatiert sein, die der GeoJSON-Spezifikation entsprechen. Typischerweise liest du ein Shapefile ein, das dein Farm-Polygon enthält, extrahierst das Geometry-Dictionary für dieses spezifische Feature, packst es in eine Standard-Python-List und übergibst es der mask-Funktion. Hier ist die wichtigste Erkenntnis. Die mask-Funktion hat zwei unterschiedliche Verhaltensweisen, die durch einen einzigen Boolean-Parameter namens crop gesteuert werden. Du musst den Unterschied verstehen, ob du ein Dataset einfach nur maskierst oder es tatsächlich cropst. Wenn du die Funktion mit crop auf false ausführst, schaut sich Rasterio dein Farm-Polygon an und setzt jeden Pixel außerhalb dieser Grenze auf einen nodata-Wert. Es liest den Default-nodata-Wert direkt aus den ursprünglichen Raster-Metadaten. Das Output-Array, das du zurückbekommst, hat exakt dieselben Dimensionen wie dein riesiges, ursprüngliches, landesweites Bild. Die Pixel innerhalb der Farm behalten ihre tatsächlichen Werte, aber alles außerhalb wird ausgeblendet. Die gesamte räumliche Ausdehnung der Datei ändert sich nicht. Wenn du crop auf true setzt, ändert sich das Verhalten. Rasterio setzt die Pixel außerhalb deines Polygons weiterhin auf nodata. Allerdings verwirft es danach die riesigen leeren Bereiche. Es berechnet die rechteckige Bounding Box deines unregelmäßigen Farm-Polygons und schrumpft das zurückgegebene Array so, dass es exakt in diese spezifische Box passt. Dein riesiges, landesweites Array wird physisch auf ein winziges, handliches Grid reduziert, das nur die Farm und ein kleines bisschen nodata-Padding an den unregelmäßigen Rändern enthält. Weil das Cropping die Dimensionen des Grids physisch verändert, ändert sich auch der geografische Startpunkt. Die obere linke Ecke ist nicht mehr die obere linke Ecke des Bundeslandes. Sie ist jetzt die obere linke Ecke der Farm-Bounding-Box. Deshalb gibt die mask-Funktion ein Tuple zurück, das zwei Elemente enthält. Das erste Element ist dein neu gecropptes NumPy-Array. Das zweite Element ist ein brandneues affines Transform-Objekt. Dieser aktualisierte Transform enthält die neuen Ursprungskoordinaten. Wenn du dieses neu gecroppte Array in einer neuen Datei speichern willst, musst du diesen aktualisierten Transform in dein Write-Profile übergeben. Wenn du aus Versehen den ursprünglichen landesweiten Transform wiederverwendest, wird dein winziges Farm-Bild wieder über das gesamte Bundesland gestreckt. Der wahre Wert der mask-Funktion ist, dass sie diese ganze räumliche Neuberechnung automatisch für dich übernimmt. Ich hoffe, das war nützlich. Danke fürs Zuhören und genieß den Rest deines Tages.

Hallo, hier ist Alex von DEV STORIES DOT EU. Rasterio: Numpy für Geodaten, Folge 5 von 5. Manchmal hast du wunderschöne, hochauflösende Satellitenbilder, aber dein Klimamodell akzeptiert nur grobe, niedrig aufgelöste Eingaben. Du kannst das Array nicht einfach zerschneiden, da die resultierenden Pixel immer noch genau dasselbe physische Gebiet darstellen müssen, nur in einem Raster mit geringerer Dichte. Um das korrekt zu machen, brauchst du Resampling und Writing. Hier geht es nicht um Reprojektion. Das Koordinatenreferenzsystem bleibt exakt gleich. Wir ändern lediglich die Auflösung. Nehmen wir an, du musst ein Raster um die Hälfte verkleinern. In Rasterio führst du das Resampling genau in dem Moment durch, in dem du die Daten in den Speicher lädst. Wenn du die read Methode für dein Dataset aufrufst, übergibst du einen Parameter namens out shape. Anstatt das gesamte Array einzulesen, liest Rasterio es ein und skaliert es sofort auf die gewünschte Form. Bei einem Einzelband-Raster wäre dein out shape ein Tuple, das ein Band, die halbe ursprüngliche Höhe und die halbe ursprüngliche Breite angibt. Du musst Rasterio außerdem mitteilen, wie die Werte dieser neuen, größeren Pixel berechnet werden sollen. Das machst du, indem du einen resampling Parameter übergibst. Rasterio stellt hierfür ein Resampling Enum bereit. Wenn du mit kontinuierlichen Daten wie Temperaturdaten arbeitest, übergibst du vielleicht Resampling dot bilinear, was einen gewichteten Durchschnitt der umliegenden Pixel berechnet. Wenn du mit kategorialen Daten wie Landbedeckung arbeitest, würdest du Resampling dot nearest verwenden, um die exakten ursprünglichen Klassenwerte beizubehalten. Hier ist die entscheidende Erkenntnis. Du hast nun ein neues, kleineres Numpy Array, aber es repräsentiert weiterhin exakt denselben geografischen Bereich auf der Erde. Da das Array weniger Pixel enthält, deckt jedes einzelne Pixel einen größeren physischen Raum ab. Das bedeutet, dass dein ursprünglicher Affine Transform nun völlig falsch ist. Wenn du das neue Array mit dem alten Transform auf die Festplatte schreibst, schrumpft dein Bild auf der Karte und deckt nur noch ein Viertel seiner ursprünglichen Fläche ab. Du musst den Transform skalieren. Das machst du, indem du den ursprünglichen Dataset Transform nimmst und ihn mit einer Affine Scaling Matrix multiplizierst. Die Skalierungsfaktoren findest du, indem du die ursprüngliche Breite durch deine neue Breite und die ursprüngliche Höhe durch deine neue Höhe teilst. Da wir die Dimensionen halbiert haben, sind beide Faktoren genau zwei. Wenn du den ursprünglichen Transform mit diesen Skalierungsfaktoren multiplizierst, bekommst du einen neuen Transform, bei dem die Pixeldimensionen verdoppelt sind, der Ursprungspunkt oben links jedoch perfekt verankert bleibt. Jetzt hast du dein resampeltes Array und deinen skalierten Transform. Der letzte Schritt ist das Schreiben eines gültigen GeoTIFFs auf die Festplatte. Dafür brauchst du Metadaten. Der sicherste Ansatz ist, die profile Property aus deinem Source Dataset zu kopieren. Dadurch bekommst du ein Dictionary, das den File Driver, das Koordinatenreferenzsystem, den Datentyp und den nodata Value enthält. Anschließend aktualisierst du dieses profile Dictionary mit deiner neuen Höhe, deiner neuen Breite und deinem neu berechneten Transform. Wenn das aktualisierte profile bereit ist, rufst du rasterio dot open auf. Du gibst den neuen Dateinamen an, setzt den mode auf write und entpackst dein profile Dictionary als Keyword Arguments. Dadurch wird ein leeres Dataset auf der Festplatte erstellt, das exakt deinen Vorgaben entspricht. Zum Schluss rufst du die write Methode auf diesem neuen Dataset Object auf und übergibst dein herunterskaliertes Numpy Array. Das Array wird in die Datei geschrieben, und dein neues GeoTIFF mit niedrigerer Auflösung ist bereit für das Klimamodell. Denk daran, das Ändern von Rasterdaten ist immer ein zweiteiliger Vertrag: Wenn du die Form des Numpy Arrays änderst, musst du den Affine Transform explizit anpassen, sonst ist dein räumlicher Fußabdruck kaputt. Damit sind wir am Ende unserer Serie angelangt. Ich ermutige dich, die offizielle Rasterio Dokumentation zu lesen, um andere Resampling Algorithmen zu erkunden und zu versuchen, diese Pipelines hands-on zu bauen. Wenn du Themen für eine zukünftige Serie vorschlagen möchtest, besuche devstories dot eu. Ich möchte mir kurz einen Moment Zeit nehmen, um dir fürs Zuhören zu danken – das hilft uns sehr. Mach's gut!