Raster Dateien mit einer Bounding Box filtern

29. Juni 2017 | Category blog

Oftmals ist es nicht nötig, alle verfügbaren Raster-Dateien zu prozessieren, sondern nur einen Teil in einem bestimmten geographischen Gebiet. Diese Region wird meist durch eine rechteckige Bounding Box definiert, d.h. „ein Gebiet, dass durch zwei Längen- und Breitengrade bestimmt wird“.

Rasterdateien finden

Zuerst möchten wir einen Überblick über die vorhandenen Dateien bekommen. In der Unix-Shell gibt es die Funktion find, welche praktischerweise alle Dateien mit der Endung „tif“ im Zielordner findet:

find /path/to/file -type f -name "*.tif" > list_of_raster.txt

1	find /path/to/file -type f -name "*.tif" > list_of_raster.txt

Mit gdalbuiltvrt anhand der Bounding Box filtern

In einem zweiten Schritt werden nun die Dateien anhand der Bounding Box gefiltert. Die freie Programmbibliothek GDAL hat für diese Problemstellung das passende Tool. Die Funktion gdalbuiltvrt erzeugt eine virtuelle Datei (VRT), welche eine Liste aller Input-Dateien enthält. Dabei ist gdalbuiltvrt sehr schnell und verbraucht kaum Speicherplatz. Die Option -te setzt die Bounding Box anhand der spezifizierten Breiten- und Längengrade. Der Befehl, um alle Rasterdateien innerhalb der Landesgrenzen von NRW zu filtern, lautet daher:

gdalbuildvrt -te 5.8887,50.3332,9.4702,52.5212 output.vrt -input_file_list list_of_raster.txt

1	gdalbuildvrt -te 5.8887,50.3332,9.4702,52.5212 output.vrt -input_file_list list_of_raster.txt

GDAL Python-Funktionen

Nun wird die GDAL Funktion GetFileList() für Python benutzt, welche die Rasterliste in der virtuellen Datei auslesen kann. Damit nur Rasterdateien ausgelesen werden, muss der erste Eintrag (output.vrt) ignoriert werden. Dies erreicht man indem der erste Wert der Dateiliste durch [1:] übersprungen wird (Python fängt bei 0 an zu zählen):

dataset = gdal.Open("output.vrt", GA_ReadOnly)
filelist = dataset.GetFileList()

#writes it to a new file
with open(dir + "/vrt_list.txt", 'w+') as f:
    for line in filelist[1:]:
        f.write(str(line) +"\n")
f.close()

dataset = gdal.Open("output.vrt", GA_ReadOnly)

filelist = dataset.GetFileList()

#writes it to a new file

with open(dir + "/vrt_list.txt", 'w+') as f:

for line in filelist[1:]:

f.write(str(line) +"\n")

f.close()

Die erstellte Test-Datei kann somit einfach weiter benutzt werden, z.B. indem man eine while-Schleife in Bash ausführt:

while read line; do
    #do something here
done < vrt_list.txt

while read line; do

#do something here

done < vrt_list.txt

Das dazugehörige Skript kann auf GitHub eingesehen werden und steht frei zur Verfügung.

Schlagwörter: bounding box, GDAL, python