Zaletą plików tekstowych jest ich jawna i otwarta struktura oraz niezależność od systemów operacyjnych i oprogramowania, co oznacza że mogą być wykorzystywane wszędzie tam, gdzie nie ma możliwości zarządzania danymi binarnymi bezpośrednio (na przykład w przypadku braku dostępu do sterowników GDAL). Każdy język programowania poradzi sobie bez problemów z każdym rodzajem pliku tekstowego, tzn. w każdym języku programowania można napisać kod który otworzy taki plik, wczyta i zinterpretuje jego zawartość jako dane, pozwoli te dane przeliczyć a następnie ponownie zapisać wynik obliczeń w nowym pliku tekstowym. Większość (wszystkie?) popularne formaty tekstowe są obsługiwane przez bibliotekę GDAL bezpośrednio, co oznacza, że nie ma potrzeby ręcznego przetwarzania danych jeżeli mamy dostęp do biblioteki GDAL. Wadą formatów tekstowych jest ich duża objętość i wolny czas odczytu/zapisu.
W poniższym przykładzie wczytamy plik "mozaika.asc" wraz ze stowarzyszonym z nim plikiem "mozaika.prj", który zawiera informację o projekcji w formacie WKT (well-known text). Format *.asc to ESRI ASCII grid, ("AAIGrid") wykorzystywany powszechnie jako przenośny format tekstowy dla danych rastrowych. Więcej na temat formatów na stronie projektu GDAL: (http://gdal.org/formats_list.html).
#Zwrócić uwagę na ścieżki dostępu.
import os
path = os.path.join(os.path.expanduser("~"),"Dropbox","dydaktyka","program_geodata","01_raster")
rastfile = os.path.join(path,"data","mozaika.asc")
print rastfile
Otwieranie fromatu tekstowego odbywa się w dokładnie taki sam sposób jak formatu binarnego. Typ pliku biblioteka GDAL rozpoznaje automatycznie. W przypadku formatu ASCII stowarzyszony plik *.prj zawiera informację na temat układu odniesienia danych w formacie WKT. Nie jest możliwy natomiast zapis/aktualizacja pliku tekstowego.
from osgeo import gdal
raster = gdal.Open(rastfile)
band = raster.GetRasterBand(1)
projection = raster.GetProjection()
raster = None # zamknięcie pliku poprzez przypisanie wartości NULL
Zaletą formatów tektowych jest możliwość samodzielnej obsługi, jeżeli poznamy jego wewnętrzną strukturę.
Wybrany przez nas format ASCII Grid posiada bardzo prostą strukturę. Informacje o danych i geotransformacji znajdują się w pliku *.asc, Jeżeli jest informacja o projekcji to znajduje się w pliku *.prj.
Wczytanie danych o projekcji o po prostu wczytanie całej zawartości pliku *.prj do zmiennej i zamknięcie pliku.
f = open(os.path.join(path,"data","mozaika.prj"), 'r')
projekcja = f.read()
f.close()
projekcja
Pierwsze sześć linii to dane geotransformacji oraz informacja na temat wartości sympolizującej brak danych. Format ASCII Grid nie pozwala ustawić osobnej rozdzielczości (cellsize) dla kierunków n-s i w-e.
W celu wczytania danych nagłówka nalezy odworzyć plik do odczytu i ręcznie wczytać (readline()) pierwsze 6 linii pliku, począwszy od 14 (13 index) znaku. Wartości należy w trakcie czytania konwertować do odpowiednich typów: całkowitego dla liczby kolumn i wierszy oraz zmiennoprzecinkowch dla pozostałych. Typ dla NODATA zależy od oczekiwanego typu danych. Pow wczytaniu wybranych linii dostęp do pliku należy zamknąć.
# header w postaci zmiennych
f = open(os.path.join(path,"data","mozaika.asc"), 'r')
ncols=int(f.readline()[13:])
nrows=int(f.readline()[13:])
xllcorner=float(f.readline()[13:])
yllcorner=float(f.readline()[13:])
cellsize=float(f.readline()[13:])
NODATA_value=int(f.readline()[13:]) #zależy od formatu danych wejściowych
f.close()
Na podstawie danych z nagłówka możemy zbudować tuplet zawierający geotransformacje
geotransforms=(xllcorner,cellsize,0,yllcorner,0,-cellsize)
geotransforms
Po nagłówku znajduje się macierz danych, gdzie kolumny oddzielone są spacjami a wiersze znakiami końca linii właściwymi dla systemu operacyjnego. Dane najlepiej wczytać funkcją loadtxt() z biblioteki numpy, pomijając 6 pierwszych wierszy nagłówka. Należy również ręcznie określić typ danych docelowych (tu Byte)
import numpy as np
arr = np.loadtxt(os.path.join(path,"data","mozaika.asc"), skiprows=6,dtype='uint8')
arr
Dane możemy przetworzyć w dowolny sposób. Tu dodamy do danych 1.
arr +=1
arr
Aby zbudować nagłówek (zmienna header) należy przygotować łańcuch tekstowy (string) składający się z:
header = 'ncols'.ljust(12)+str(ncols)+os.linesep + \
'nrows'.ljust(12)+str(nrows)+os.linesep + \
'xllcorner'.ljust(12)+str(xllcorner)+os.linesep + \
'yllcorner'.ljust(12)+str(yllcorner)+os.linesep + \
'cellsize'.ljust(12)+str(cellsize)+os.linesep +\
'NODATA_value'.ljust(12)+str(NODATA_value)
header # nie będzie widać przerw między liniami
Jeżeli chcemy aby nasze dane miały układ odniesienia musimy albo użyć zbioru z danych wejściowych, albo przygotować własny na podstawie biblioteki osr (więcej o projekcjach w innych cześciach kursu). Do zbudowania pliku projekcji potrzebujemy kodu EPSG (użyjemy 2180). W pierwszym kroku należy zbudować Referencję przestrzenną a następnie zaimportować do niej parametry na podstawie kodu EPSG. Następnie zapisać wynik w pliku *.prj w formacie WKT (ExportToWkt())
from osgeo import osr
proj = osr.SpatialReference()
proj.ImportFromEPSG(2180)
with open(os.path.join(path,'mozaikaMan.prj'), 'w') as file: # ...Man od Manual
file.write(proj.ExportToWkt()) # używając with.... nie musimy martwić się o zamykanie pliku (pythonic way)
Ostatnim krokiem jest zapis pliku z danymi. W tym celu użyjemy analogicznej funkcji do loadtxt: savetxt(). Obok nazwy pliku i źródła danych, funkcja posiada kilka istotnych parametrów:
np.savetxt(os.path.join(path,'mozaikaMan.asc'),arr,fmt='%i',header=header,comments='')