# setup
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
import warnings
warnings.filterwarnings('ignore')
from osgeo import gdal
import os
import numpy as np
import matplotlib.pyplot as plt
Wizualizacja danych rastrowych¶
Wizualizacja danych geoprzestrzenych w Pythonie obejmuje co najmniej dwa zagadnienia:
- wizualizację informacji przestrzennej (w postaci map lub serii map)
- Wizualizację informacji nieprzestrzennej - w postaci diagramów, wykresów itp.
Podstawowy moduł GDAL nie zawiera żadnych narzędzi do wizualizacji danych geoprzestrzennych. Python nie posiada w struktur danych które były by odpowiedzialne bezpośrednio za dane geoprzestrzenne (poza mało przydatnym rasterio). Są one importowane bezpośrednio do formatu numpy array i tym samym cała informacja geoprzestrzenna zostajOznacza to że jeżeli chcemy wizualizować dane kategoryzowane - tabelę kolorów musimy importować osobno do formatu zrozumiałego dla języka Python jest o tyle istotne, że jeżeli chcemy używać tabeli barwnej w prezentacjach nieprzestrzennych musimy korelować kody barwne pomiędzy mapami a wykresami.
W ramach tej części kursu zapoznamy się z następującymi zagadnieniami:
- jak wydobywać wbudowane tabele barwne z plików geoprzestrzennych jak wizualizować pojedyncze mapy i serie map wraz z ich własnymi paletami kolorów (jeżeli dostępne)
- jak kolerować tabele barwne dla map i dla wykresów wykonywanych przy pomocy bibliotek seaborn i matplotlib
Import tabeli barwnej¶
tabela barwna jest własnością warstwy(band) i jest pozyskiwana funkcją band.GetRasterColorTable(). Jeżeli tabela nie jest dostępna - na przkład jest to plik numeryczny funkcja zwraca None. W ten sposób możemy sprawdzić czy tabela barwna istnieje. Tabel jest przechowywana jako krotka trzech lub czterech wartości (RGBA), natomiast dla języka Python lepiej jest przechowywać kod koloru w postaci szesnastkowej #RRGGBBAA. W tym celu można krotkę barwną zmodyfikować poprzez fomratowanie łańcucha tekstowego: "#{0:02x}{1:02x}{2:02x}".format(*ct.GetColorEntry(i), gdzie 02x oznacza prezentowanie wartości jako dwuznakowego kodu szesnastkowego a* przed wywołaniem funkcji to operator rozwinięcia krotki.
e zredukowana do postaci macierzy liczb. Pozostałe elementy takie jak zasięg, typ danych, tabela kolorów są importowane do osobnych struktur.
Systemy informacji geograficznej używają kolorów do prezentacji informacji dla człowieka. Sposób obsługi różni się w zależości od typu danych przestrzennych
- dane kategoryzowane - reprezentowane przez typ unit8 i unit16 (nieujemne), które mogą przechowywać wewnętrznie tabelę kolorów
- dane numeryczne - reprezentowane przez pozostałe typy danych, które nie posiadają wbudowanej tabeli kolorów i są wizualizowane przy pomocy wskazanej tabeli barw.
Oznacza to że jeżeli chcemy wizualizować dane kategoryzowane - tabelę kolorów musimy importować osobno do formatu zrozumiałego dla języka Python jest o tyle istotne, że jeżeli chcemy używać tabeli barwnej w prezentacjach nieprzestrzennych musimy korelować kody barwne pomiędzy mapami a wykresami.
W ramach tej części kursu zapoznamy się z następującymi zagadnieniami:
- jak wydobywać wbudowane tabele barwne z plików geoprzestrzennych jak wizualizować pojedyncze mapy i serie map wraz z ich własnymi paletami kolorów (jeżeli dostępne)
- jak kolerować tabele barwne dla map i dla wykresów wykonywanych przy pomocy bibliotek seaborn i matplotlib
Import tabeli barwnej¶
tabela barwna jest własnością warstwy(band) i jest pozyskiwana funkcją band.GetRasterColorTable(). Jeżeli tabela nie jest dostępna - na przkład jest to plik numeryczny funkcja zwraca None. W ten sposób możemy sprawdzić czy tabela barwna istnieje. Tabel jest przechowywana jako krotka trzech lub czterech wartości (RGBA), natomiast dla języka Python lepiej jest przechowywać kod koloru w postaci szesnastkowej #RRGGBBAA. W tym celu można krotkę barwną zmodyfikować poprzez fomratowanie łańcucha tekstowego: "#{0:02x}{1:02x}{2:02x}".format(*ct.GetColorEntry(i), gdzie 02x oznacza prezentowanie wartości jako dwuznakowego kodu szesnastkowego a* przed wywołaniem funkcji to operator rozwinięcia krotki.
os.chdir('/home/jarekj/Dropbox/dydaktyka/programowanie2/landData')
raster = gdal.Open('10_landuse.tif')
band = raster.GetRasterBand(1)
ct = band.GetRasterColorTable()
catlist = np.unique(band.ReadAsArray())
cols=[]
for i in catlist.tolist():
cols.append("#{0:02x}{1:02x}{2:02x}".format(*ct.GetColorEntry(i)))
cols
Wizualizacja pojedynczej mapy z paletą barwną.¶
Wizualizację pojedynczej mapy bez palety barwnej - takie jak pobór danych, wyznaczenie zasięgu i wizualizacja macierzy została zaprezentowana w poprzedniej części kursu dlatego w tej części ograniczymy się jedynie do omówienia tych zagadnień, które dotyczą palety barwnej. Całość zostanie zamknięta w postaci uproszczonej funkcji, której efektem jest wyświetlenie mapy.
Po nawiązaniu połączenia z plikiem i pobraniu warstwy (band) pobierana jest tabela. Jeżeli istnieje, z danych pobierane są wartości unikalne lub histogram. Następnie dla każdej wartości unikalnej lub pozycji histogramu, która zawiera elementy budowany jest kod barwny i dodawany do listy kolorów. Jeżeli chcemy tę paletę kolorów użyć w bibliotece matplotlib musimy ją skonwertować do ListedColorMap z modułu colors biblioteki matplotlib. Ostatnim krokiem jest wyświetlenie mapy przy pomocy znanego już polecenia plt.imshow() ze wskazaniem krotki zawierającej zasięg oraz paletę kolorów.
from matplotlib import colors
def plot_raster_by_file(file,ax=None,cmap=None):
'''
docstring
'''
ax = ax or plt.gca() # pobranie aktywnego axes
rast = gdal.Open(file)
if not rast:
print("data not found or not raster data")
return None
gt = rast.GetGeoTransform()
east,south = gdal.ApplyGeoTransform(gt,rast.RasterXSize,rast.RasterYSize) # prawy dolny
west,north = gdal.ApplyGeoTransform(gt,0,0) #lewy górny
extent = [west,east,south,north]
band = rast.GetRasterBand(1)
arr = band.ReadAsArray()
ct = band.GetRasterColorTable()
if ct and not cmap:
arr[arr == band.GetNoDataValue()]=0
hist = band.GetHistogram()
cols = ['#FFFFFF'] # simulate NA
for i in range(ct.GetCount()):
if hist[i]:
cols.append("#{0:02x}{1:02x}{2:02x}".format(*ct.GetColorEntry(i)))
cmap = colors.ListedColormap(cols)
ax.imshow(arr,cmap=cmap,extent = extent)
Wyświetlanie kolekcji map¶
Aby wyświetlić kolekcję map musimy skorzystać z poprzednio utworzonej funkcji a następnie zastosować narzędzie suplots() z biblioteki matplotlib. Jak to wykonać?
W pierwszym kroku musimy pobrać listę wszystkich plików, które chcemy wizualizować. Gdy lista jest gotowa (może być wykonana ręcznie), możemy wykonać naszą rycinę (fig), którą definiujemy jako serię podfigur (subplots). Każda figura zawiera wykres (axes), która jest dodawana do pierwszej wolnej podfigury. Funkcja subplots zwraca dwa elementy: figurę oraz listę wykresów (axes) do których przy pomocy pętli dodajemy kolejne mapy. W naszym przykładzie zdefiniowaliśmy macierz wykresów o wymiarach 2 wiersze x 5 kolumn oraz wymiary obrazka 15 x 5 cali. Niestety, nie ma jednoznacznego klucza co do doboru wielkości, parametry te należy dobrać ręcznie.
Często niepożądanym efektem takich wykresów jest nakładanie się legend. Niestety, obrócenie etykiet pod pewnym kątem nie jest proste, a właściwie nie jest łatwe do znalezienia dla kolekcji wykresów. Aby to wykonać, należy po wykonaniu figury ale przed jej uruchoniemiem, dla każdego wykresu (axes) zmodyfikować rotację xticks to prządanej wartości (tu 45 stopni). Następnie można wyświetlić figurę.
Note: funkcje: sca i gca¶
Powyższe funkcje to sktóry od set current axes i get current axes są bardzo pomocne w sytuacji gdy włączymy wykresy na jednej figurze.
Funkcje pochodzą z modułu pyplot biblioteki matplotlib. gca() pobiera aktywny wykres (axes) zgodny z argumentami, jeżeli nie może takiego pobrać tworzy nowy. sca() ustawia aktwne axes jako domyślne w danej figurze.
Kolory: Problem wartości pustych.¶
W plikach geoprzestrzennych wartość pusta jest reprezentowana przez wybraną liczbę, która w sposób jawny jest wskazana we właściwościach pliku geoprzestrzennego. Dla pliku numerycznego może do być liczba symbolizująca NaN (Not a Number) często jest to jednak liczba, której wartość nie występuje w zbiorze danych (często stosowana jest -9999). Dla danych kategoryzowanych może to być wartość 0 (nie ma kategorii zerowej) ale często stosuje się wartość 255. Ponieważ kategorii jest z reguły nie więcej niż kilkanaście w zbiorze danych, wartość 255 rozciąga nam histogram na pełen zakres wartości, co oznacza problemy z dopasowaniem palety kolorów. W takich sytuacjach lepiej jest zastąpić wartości "NULL" w pliku wartością 0 i przypisać dla tej barwy kolor biały. Natępnie dla pozostałych kategorii przypisujemy kolory we właściwej im kolejności.
os.chdir('/home/jarekj/Dropbox/dydaktyka/programowanie4/landData')
files=[f for f in os.listdir(".") if f.endswith('.tif')]
tmpraster = gdal.Open(files[2])
div = 1
xsize,ysize=(int(np.round(tmpraster.RasterXSize/div)),int(np.round(tmpraster.RasterYSize/div)))
files.sort()
files
# Uruchomić całość w jednym kroku
fig, axes = plt.subplots(2, 5, sharex='all', sharey='all',figsize=(15,5))
for ax,file in zip(axes.flat[0:10],files):
plot_raster_by_file(file,ax)
for ax in fig.axes:
plt.sca(ax) # atywowanie axes
plt.xticks(rotation=45)
