Zarządzanie danymi rastrowymi - Python

Dane rastrowe obsługiwane są przez bibliotekę GDAL. Biblioteka dostarcza narzędzi do otwarcia, zapisu i manipulacji ponad 100 formatów rastrowych. Lista formatów znajduje się na stronie projektu GDAL: (http://gdal.org/formats_list.html).

Proces zarządzania danymi rozpoczyna się od ustalenia lokalizacji pliku i przygotowania dla niego ścieżki dostępu. Należy pamietać o użyciu ścieżki dostępu właściwej dla systemu operacyjnego oraz uniwersalnej nazwy użytkownika. Skrypty pythona są znacznie częściej używane poza środowiskiem Windows.

Definicja rastra

Geoprzestrzenny obiekt rastrowy to obraz rozszerzony o parametry geoprzestrzenne pozwalające na lokalizację każdego pixela (komórki rastra) w odniesieniu do przestrzenii geograficznej.

In [2]:
import os #dostosować ścieżkę dostępu
path = os.path.join(os.path.expanduser("~"),"Dropbox","dydaktyka","program_geodata","01_raster") 
rastfile = os.path.join(path,"data","dem.tif")
print rastfile

/home/jarekj/Dropbox/dydaktyka/program_geodata/01_raster/data/dem.tif

Otwarcie pliku

Do zarządzania danymi będziemy potrzebowali biblioteki GDAL, którą trzeba zaimportować (form osgeo). Następnie przy pomocy plecenia open() otwieramy połączenie z plikem. Otwarcie połączenia nie wczytuje danych.

In [3]:
from osgeo import gdal
raster = gdal.Open(rastfile)

Dostępne funkcje można sprawidzć poleceniem dir(raster). Najważniejsze informacje można pozyskać na temat obiketu raster to przy pomocy wybranych funkcji (z () na końcu) lub atrybutów (bez ()):

Nazwa Sterownika

Pozwala określić jakiego typu jest plik z danymi wejściowymi. Często zależy nam aby zapisać wyniki w tym samym formacie

In [4]:
raster.GetDriver().LongName
Out[4]:
'GeoTIFF'

Liczba pasm

Wiele zbiorów geoprzestrzennych posiada więcej niż jedno pasmo. Pliki RGB 3, multispektralne - więcej. Minimalna liczba to 1. Będzie wykorzystywana przy pobieraniu danych.

In [5]:
raster.RasterCount
Out[5]:
1

Ścieżka dostępu

Określa gdzie zapisany jest plik wejściowy

In [6]:
raster.GetDescription()
Out[6]:
'/home/jarekj/Dropbox/dydaktyka/program_geodata/01_raster/data/dem.tif'

Metadane

Podstawowe informacje o pochodzeniu pliku

In [7]:
raster.GetMetadata()
Out[7]:
{'AREA_OR_POINT': 'Area',
 'TIFFTAG_SOFTWARE': 'GRASS GIS 7.2.2 with GDAL 2.2.1'}

Rozmiary pliku:

ilość wierszy (XSize) i kolumn (YSize)

In [8]:
raster.RasterXSize,raster.RasterYSize # rozmiary zwraca tuplet
Out[8]:
(93, 83)

Projekcja

Elipsoida i układ odniesienia

In [9]:
raster.GetProjection() # projekcja
Out[9]:
'PROJCS["Transverse Mercator",GEOGCS["grs80",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",19],PARAMETER["scale_factor",0.9993],PARAMETER["false_easting",500000],PARAMETER["false_northing",-5300000],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]'

Geotransformacje

Parametry pozwalające przeliczyć siatkę rastra (macierz wartości) na punkty w przestrzeni geograficznej:

  • 0 zachodnia krawędź (minimum x)
  • 1 przesunięcie po osi x od 0 (z reguły dodatnie)
  • 2 przesunięcie początku komórki względem krawędzi 0, 0 narożnik, 1/2 * 1 środek komórki
  • 3 północna krawędź (maximum y)
  • 4 przesunięcie początku komórki względem krawędzi 3, 0 narożnik, 1/2 * 5 środek komórki
  • 5 przesunięcie po osi y od 3 (z reguły ujemne)
In [10]:
gt=raster.GetGeoTransform() # geotransformacje
gt
Out[10]:
(350400.0, 200.0, 0.0, 516600.0, 0.0, -200.0)

Parametry geo-transformacji pozwalają przeliczyć dowolne współrzędne siatki (wiersz i kolumna) na współrzędne geograifczne, na przykład dla komórki col,row=(7,2)

In [11]:
col,row=(7,2)
x = gt[0] + col * gt[1] + row * gt[2]
y = gt[3] + col * gt[4] + row * gt[5]
x,y
Out[11]:
(351800.0, 516200.0)

Lub poprzez funkcję ApplyGeoTransform()

In [12]:
x,y=gdal.ApplyGeoTransform(gt,col,row)
x,y
Out[12]:
(351800.0, 516200.0)

W przeciwnym kierunku:

In [13]:
col = int((x - gt[0]) / gt[1]) 
row = int((y - gt[3]) / gt[5])
col,row
Out[13]:
(7, 2)

Pobranie danych

W celu pobrania danych należy wskazać pasmo, z którego będziemy pobierać dane. Wszystkie pasma posiadają ten sam zakres przestrzenny i rozdzielczość (to cecha rastra) ale mogą posiadać inny typ danych i inną reprezentację wartości pustej (NULL). Wybrane dostępne właściwości to:

In [14]:
band = raster.GetRasterBand(1)

Typ danych

Dostępne typy danych:

GDT_Unknown = 0, GDT_Byte = 1, GDT_UInt16 = 2, GDT_Int16 = 3, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6, GDT_Float64 = 7, GDT_CInt16 = 8, GDT_CInt32 = 9, GDT_CFloat32 = 10, GDT_CFloat64 = 11, GDT_TypeCount = 12

In [15]:
band.DataType #zwraca integer,
Out[15]:
6

Podstawowe statystyki

Dla każdego pasma można określić podstawowe statystyki: Min, Max, Avg, Sd

In [16]:
band.GetStatistics(1,0) # zgodność z orginalną funkcją C, pierwszy parametr: akceptacja wartości przybliżonych, drugi wymuszenie przeliczania
Out[16]:
[0.0, 0.0, 0.0, -1.0]

Wartość pusta

Zwraca wartość symbolizującą wartość pustą

In [17]:
band.GetNoDataValue() #Nieokreślona

Histogram

Do pełniejszej charakterystyki warstwy możemy użyć histogramu, który jest jedną z najpopularniejszych metod charakteryzowania rozkładu zmiennej. W trakcie obliczania hitogramu można podać zakres wartości i ilość przedziału historamu w wyniku otrzymamy listę liczb całkowitych określających liczbę komórek w każdym przedziale.

In [18]:
import math
minVal,maxVal,meanVal,sdVal=band.GetStatistics(0,1) # (wyjasnic expand tuple)
nbin=100
hist=band.GetHistogram(math.floor(minVal),math.floor(maxVal),nbin) # min, max , planowana liczba
print hist

[1, 1, 5, 32, 16, 20, 48, 41, 38, 29, 28, 30, 24, 28, 28, 31, 43, 46, 29, 41, 45, 35, 27, 39, 50, 52, 72, 71, 80, 108, 86, 89, 97, 129, 170, 162, 123, 161, 126, 121, 85, 76, 81, 62, 83, 58, 71, 76, 73, 62, 61, 65, 71, 86, 88, 84, 64, 50, 23, 32, 28, 27, 17, 15, 21, 17, 12, 19, 11, 12, 7, 16, 16, 12, 14, 4, 8, 6, 5, 3, 4, 6, 6, 5, 8, 4, 4, 5, 3, 7, 5, 2, 4, 2, 2, 1, 0, 0, 0, 0]

W następnym kroku przygotujemy pomocnicze elementy do wyświetlenia histogramu. Użyjemy bliblotek numpy i matplotlib do przygotowania wykresu:

In [19]:
import numpy as np
x=np.arange(nbin) # zbudowanie listy przedziałów
xtics=np.linspace(math.floor(minVal),math.floor(maxVal),nbin/10) # zbudowanie listy etykiet (co dziesiąty przedział)
xat=np.linspace(0,nbin,nbin/10) # zbudowanie listy lokalizacji etykiet

import matplotlib.pyplot as plt
plt.bar(x, np.asarray(hist),linewidth=0,width=1,color='#A11207') #dodanie histogramu
plt.xticks(xat, xtics.astype(int)) # dodanie etykiet jako liczb całkowitych, brak cudzysłowia
plt.show()

/usr/lib/python2.7/dist-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

Pobranie danych w formie tablicy numpy

Dane do daleszego przetwarzania najprościej jest pobrać w postaci dwuwymiarowej tablicy numpy, nunkcją ReadAsArray(). Należy uważać na wielkość pliku wiejściowego.

In [20]:
arr=band.ReadAsArray()
arr
Out[20]:
array([[  80.21151733,   82.81385803,   82.42944336, ...,  108.05825806,
         105.45590973,  105.29399109],
       [  80.43751526,   83.74990845,   85.00001526, ...,  107.35435486,
         104.8393631 ,  100.41969299],
       [  76.46885681,   82.84009552,   85.81338501, ...,  104.08059692,
         100.50229645,  100.85406494],
       ..., 
       [  78.47695923,   78.66946411,   78.7088623 , ...,   87.54154968,
          87.5704422 ,   84.02039337],
       [  78.56521606,   78.70487976,   78.75489044, ...,   87.21909332,
          87.30551147,   85.19948578],
       [  77.95922852,   78.79443359,   78.81874084, ...,   84.71836853,
          86.27947998,   85.21043396]], dtype=float32)

Wyświetlenie rastra w postaci macierzy

Wizualizacja danych jest jednym z podstawowych kroków analizy. Biblioteka matplotlib daje tu nieograniczone możliwości.

Proste wyświetlenie macierzy

Macierz zostanie wyświetlona w domyślnych kolorach i we wspórzędnycyh macierzy

In [21]:
plt.imshow(arr)
plt.show()

Modyfikacja osi i tabeli kolorów

Osie zostaną zmodyfikowane do ukladu odniesienia oraz dodany zostanie pasek legendy. Dla kolorów zostanie wybrana inna paleta: (https://matplotlib.org/examples/color/colormaps_reference.html) Do wyznaczenia zasięgów osi, zostaną wykorzystane geotransformacje:

In [22]:
east,south = gdal.ApplyGeoTransform(gt,raster.RasterXSize,raster.RasterYSize) # prawy dolny
west,north = gdal.ApplyGeoTransform(gt,0,0) #lewy górny

plt.imshow(arr,cmap='viridis',extent=[west,east,south,north])
plt.colorbar()
plt.show()

Wyświetlenie fragmentu rastra i dodanie wartości jako etykiet

Przydatna opcja jeżeli zależy nam na poznaniu dokładych wartości. Dostosujemy zasięgi poprzez modyfikację geotransformat. Należy wyłączyć domyślną interpolację aby zobaczyć pojedyncze komórki. inne przykłady interpolacji można znaleźć w linku. (https://matplotlib.org/examples/images_contours_and_fields/interpolation_methods.html)

Dodanie wartości (zaokrąglonych) komórek zostanie wykonane przy pomocy funkcji text. W takiej sytuacji nie możemy zastosować współrzędnych geograficznych (przynajmniej w prosty sposób). W celu nadpisania wartości została wykonana tablica wspórzędnych (np.meshgrid), która następnie została połączona we wspólną listę współrzędnych.

In [23]:
col1,row1=(40,10) #lewy górny róg
col2,row2=(50,20) # prawy dolny róg

arrP=arr[col1:col2,row1:row2]
#fig,ax=plt.subplots()
 
plt.imshow(arrP,cmap='viridis',interpolation='none') #Nie można użyć extend

ind_array = np.arange(0, 10, 1)
x, y = np.meshgrid(ind_array, ind_array)
for x_val, y_val in zip(x.flatten(), y.flatten()):
    plt.text(x_val, y_val, '%.0f' % arrP[x_val, y_val], va='center', ha='center')

plt.colorbar()
plt.show()

Tworzenie rastra od podstaw

Aby utworzyć raster (grid) musimy znać docelowy rozmiar, parametry geotransformacji (north,west,resolutution) oraz układ odniesienia. W ramach prezentacji wykonamy tablicę o rozmiarze $100 \times 100$,dla którego następnie zostaną wyliczone pozostałe parametry. Taka tablica nie ma żadnego odniesienia geoprzestrzennego, jest jedynie zbiorem liczb poukładanych w wiersze i kolumny. Istnieje kilka funkcji inicjujących taką macierz, my użyjemy funkcję zeros(), która wypełnia tablicę zerami, jako typ przyjmiemy unit8 czyli Byte. W następnych krokach tablica zostanie wypełniona różnymi wartościami, tworząc mozaikę. Wartości przypisywane są wycinkom tablicy, wybieramym poprzez indeksowanie:

In [24]:
arr=np.zeros((100,100),dtype='uint8')  # tuplet jako argument

arr[10:20,15:38]=1
arr[36:78,21:66]=2
arr[81:100,81:100]=3
arr[0:3,0:3]=4
arr[45:72,11:44]=6
plt.imshow(arr)
plt.show()

Układ odniesienia

Zostanie utworzony przy pomocy funkcji ImportFromEPSG(), która buduje odpowiedni układ na podstawie podanego kodu (użyjemy 2180). Więcej o układach odniesienia w osobnej części kursu.

In [25]:
from osgeo import osr
epsg=2180
spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(epsg) 
spatialRef.ExportToWkt()
Out[25]:
'PROJCS["ETRS89 / Poland CS92",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",19],PARAMETER["scale_factor",0.9993],PARAMETER["false_easting",500000],PARAMETER["false_northing",-5300000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","2180"]]'

Parametry geotransformacji i tworzenie topologii rastra

Po wyborze parametrów geotransformacji tworzymy plik rastrowy. W pierwszej kolejności wybieramy format pliku (driver - tu użyjemy najbardziej uniwersalny GeoTIFF), a następnie:

  • tworzymy plik, jako rozmiar rastra podajemy parametry tablicy (shape).
  • Określamy też typ danych docelowych, użyjemy Byte, tak, jak wartości w tablicy. W przypadku niezgodności typów danych wystąpią problemy.
  • nadajemy plikowi projekcję
  • dodajemy geotranformacje - najważniejszą część aby plik graficzny (Tiff) mógł stać się geoprzestrzennym (GeoTifff)

Taki plik nie zawiera danych.

In [26]:
north=450000
west=520000
res=250

drv=gdal.GetDriverByName('GTiff')
target = drv.Create("mozaika.tif", arr.shape[0],arr.shape[1], 1, gdal.GDT_Byte)
target.SetProjection(spatialRef.ExportToWkt())
target.SetGeoTransform((west, res, 0, north, 0, -res))
Out[26]:
0

Dodanie danych

na tym etapie należy określić jaka wartość będzie reprezentowała wartość pustą. Użyjemy maxymalnej, 255. następnie musimy pobrać pasmo danych (band) i przypisać do niego wartość NODATA, dane WriteArray() oraz wyczyścić chache.

In [27]:
noData=255
tband=target.GetRasterBand(1)
tband.SetNoDataValue(noData)
tband.WriteArray(arr)
tband.FlushCache()
tband.ComputeStatistics(True)
Out[27]:
[0.0, 6.0, 0.9233, 1.7958332634184055]

Zamknięcie pliku

Ostatnim krokiem jest zamknięcie pliku. Realizujemy to poprzez przypisanie obiektowy rastra wartoci pustej (None)

In [28]:
target=None

Kopiowanie rastra - CreateCopy

Aby utworzyć nowy obiekt rastrowy można skorzystać z funkcji CreateCopy z biblioteki GDAL, która pozwala utworzyć kopię otwartego pliku rastrowego, zawierającego wszystkie niezbędne informację, w tym również projekcję, zasięg, rozdzielczość itp. Ponieważ w analizach geoprzestrzennych najczęściej obiekt wynikowy będący efektem przetworzenia danych naczęściej obejmuje ten sam obszar co dane źródłowe jest to szybka forma przygotowania odpowiedniego zbioru źródłowego. Zarówno format źródłowy jak i docelowy są od siebie niezależne, ale nie wszystkie formaty są obsługiwane.

Po utworzeniu kopii rastra można nadpisać dane a następnie przeładować Cache. Ostatnim krokiem jest zamknięcie obu dowiązań do obiektów (rasterCopy = None)

In [29]:
# Orginalne dane
raster = gdal.Open(rastfile)
dane = raster.GetRasterBand(1).ReadAsArray()


# kopia danych 
driver = gdal.GetDriverByName("GTiff")
rasterCopy = driver.CreateCopy(os.path.join(path,"demCopy.tif"),raster,strict=0)

# Przetwarzanie danych
dane = dane / 0.2991  # przeliczenie na stopy

# Zapis danych
rasterCopy.GetRasterBand(1).WriteArray(dane)
rasterCopy.GetRasterBand(1).FlushCache()

# Zamknięcie plików
raster = None
rasterCopy = None
del dane