# setup
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
import warnings
warnings.filterwarnings('ignore')
import sqlite3
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
os.chdir('/home/jarekj/Dropbox/dydaktyka/programowanie4/powiatData')
Czytanie kolekcji danych rastrowych do pandas dataFrame¶
Modelowanie predykcyjne stosuje się również, a może przede wszystkim dla danych rastrowych. Wymaga to jednak przetworzenia kolekcji plików rastrowych do postaci tabelarycznej, gdzie każdy wiersz będzie odpowiadał komórce o określonej lokalizacji a każda kolumna określonej warstwie rastrowej. Tak zorganizowane dane noszą nazwę uporządkowanych (ang. tidy) i podobnie jak atrybuty piku wektorowgo odpowiadają wymogą pierwszej postaci normalnej. Oczywiście zakładamy, że dane spełniają wymóg zgodności geoprzestrzennej, który na tym etapie nie jest sprawdzany (zob. poprzednie części kursu na temat sprawdzania zgodności geoprzestrzennej). Poniważ dane rastrowe składają się z bardo dzużej ilości komórek, przed zamianą do postaci dataFrame należy w pierwszej kolejności ocenić wielkości przetwarzanego pliku aby nie zbudować ramki danych która nie będzie mogła być obsłużona w pamięci. W przypadku, gdy dane są zbyt duże należy zastosować inne rozwiązania (kótre zostaną omówione w innej części kursu). Należy pamiętać, że dane rastrowe zbudowane są tylko z danych liczbowych, które mogą przechowywać zarówno wartości (liczby zmiennoprzecinkowe i "duże" typy całkowitoliczbowe) oraz kategorie (typy całkowite bez znaku: unit8 i unit16) przechowujace tylko wartości nieujemne. Oznacza to że w jednej dataFrame możemy mieć kolumny o różnych typach danych.
Kolejność importu kolekcji warst rastrowych jest następująca:
- Należy przygotować listę plików które zostaną przetworzone do. Listę można przygotować ręcznie lub poprzez skanowanie katalogu
- Pandas wymagają nazw kolumn. W tym celu z nazw plików trzeba usunąć rozszerzenia aby przygotować nazwy kolumn dla dataFrame (nie są dodawane automatycznie)
- Należy przygotować pusty słownik: jest to najprostrza metoda budowania dataFrame
- Następnie w pętli otwieramy po kolei wszystkie pliki i wczytujemy ich zawartość do kolejnych kluczy słownika (A) w formie spłaszczonej (B), gdzie kluczem jest nazwa kolumny a wartością - odczytany numpy array
- Po zakończeniu pętli przekształcamy słownik do dataFrame
Jako przykład zostannie przeanalizowana kolekcja danych zaierająca zarówno dane kategoryzowane jak i numeryczne. Kolekcja zostanie przygotowana przy pomocy pętli skanującej katalog i pobierającej wszystkie pliki z rozszerzeniem *tif.
from osgeo import gdal
os.chdir("/home/jarekj/Dropbox/dydaktyka/programowanie4/landData/")
files=[f for f in os.listdir(".") if f.endswith('.tif')] #1
files.sort()
files
colnames=[os.path.splitext(file)[0] for file in files] #2
dct = {} #3
for file,colname in zip(files, colnames): #4A
raster = gdal.Open(file)
dct[colname] = raster.GetRasterBand(1).ReadAsArray().flatten() #4B
df = pd.DataFrame(dct)
df.sample(8)
df.info()
Note: Komórki które nie zawierają wartości¶
W zbiorach danych liczbowych i kategoryzowanych wartości reprezentowane są przez liczby oznaczające wartości lub kategorie. Często jednak w danym miejscu wartości nie są zdefiniowane. W przypadku danych mamy następujące możliwości:
- wartość w danym miejscu nie jest znana lub nie należy do zbioru
- wartości w danym miejscu nie udało się wyliczyć (najczęściej jako efekt niedozwolonej operacji matematycznej: na przykład pierwiastek z liczby ujemnej) lub jest nieskończona - na przykład w wyniku dzielenia przez 0
Wartości nie będące liczbą nie powinny występować i są reprezentowane pzez wartości puste. W danych geoprzestrzennych wartości puste mogą być one reprezentowane jako "nie liczba" (NaN, w praktyce wartość reprezentująca nie-liczbę) lub może to też być dowolna wartość liczbowa która nie występuje w zakresie danych i symbolizuje wartość pustą. Często jest to wartość -9999 w plikach z wartościami numerycznymi. W przypadku danych całkowitych bez znaku (uint8 i uint16), które są stosowane do danych kategoryzowanych są to wartości 0 lub 255.
W języku Python w numpy nie ma możliwości reprezentowania wartości pustej przez nie-liczbę dla typów całkowitych. W zamian stosuje się liczbę całkowitą - najczęściej 255. W każdym przypadku, niezależnie od typu danych gdy wartość pusta jest reprezentowana przez liczbę po wczytaniu danych do dataFrame będzie ona zwykłą wartością. Aby wskazać ją jako wartość pustą, należy ręcznie zamienić wszystkie wystąpienia liczbowe reprezentujące wartość pustą na symbole not-a-number.
Aby wprowadzić do danych wartości puste należy ponownie przetworzyć listę plików i po otwarciu pliku pobrać informację o NodataValue do zmiennej. Następnie należy sprawdzić czy zmienna jest zdefiniowana (nie jest None) i czy jest liczbą (nie jest NaN). Jeżeli nie jest zdefiniowana ma wartość None oznacza to że albo w zbiorze danych nie ma wartości pustych albo plik został źle przygotowany. jeżeli natomiast jest zdefiniowana należy przeprowadzić test czy wartość jest liczbą. Najprościej wykonać go poprzez porównanie zmiennej z... samą sobą. Otóż jeżeli porównamy dwie nie-liczby, w wyniku otrzymamy zawsze wartość False. Zakładając że zmienna ma nazwę nodata test wygląda następująco: if nodata and nodata == nodata: Jeżeli zmienna jest zdefiniowana i jest reprezentowana przez liczbę należy zamienić ją na np.NaN w dataFrame. Problem polega na tym, że zamiana serii liczbowej na np.nan spowoduje kowersję typu danych do float64 (numerycznego) co nie jest porządane. W związku z tym kolumny co do których uważamy że przechowują dane kategoryzowane należy zamienić na typ kategory. Uznajemy że dane kategoryzowane ograniczajś się do typu uint8 (_GDTByte).
W tym celu dla każdej kolumny, sprawdzamy czy nie dane nie są typu Byte i zamieniamy ją na typ kategoryzowany (astype('category')). Następnie, niezależnie od typu danych każda która posiada wartość pustą reprezentowaną przez liczbę należy wszystkie wystąpienia liczby zamienić na np.NaN. Wykorzystujemy do tego metodę df.replace(): df[colname] = df[colname].replace(nodata,np.NaN)
for file,colname in zip(files, colnames):
raster = gdal.Open(file)
band = raster.GetRasterBand(1)
if band.DataType == gdal.GDT_Byte:
df[colname] = df[colname].astype('category')
nodata = band.GetNoDataValue()
if nodata and nodata == nodata: # test czy jest to liczba
df[colname] = df[colname].replace(nodata,np.NaN)
df.info()
Oczywiście całą procedurę najlepiej zamknąć w funkcji. Funkcja jako argument przyjmuje listę nazw plików w formie rozwijanej krotki. Jeżeli zamierzamy przekazać gotową listę należy ją poprzedzić operatorem rozwijania krotki: *.
def rasterToDataFrame(*files):
dct = {}
colnames=[os.path.splitext(file)[0] for file in files]
for file,colname in zip(files, colnames):
raster = gdal.Open(file)
dct[colname] = raster.GetRasterBand(1).ReadAsArray().flatten()
df = pd.DataFrame(dct)
for file,colname in zip(files, colnames):
raster = gdal.Open(file)
band = raster.GetRasterBand(1)
if band.DataType == gdal.GDT_Byte:
df[colname] = df[colname].astype('category')
nodata = band.GetNoDataValue()
if nodata and nodata == nodata: # test czy jest to liczba
df[colname] = df[colname].replace(nodata,np.NaN)
return df
data = rasterToDataFrame(*files)
data.info()
Zapis kolumny data frame do pliku rastrowego¶
Aby zapisać wybraną kolumnę do pliku rastrowego należy wykonać czynności w odwrotnej kolejności co do jego czytania. Dodatkowo należy określić wymiary rastra i docelowy typ danych. Najprostszym sposobem pozyskania informacji o rozmiarze pliku jest użycie jednego z plików wejściowych jako źródła rozmiaru, projekcji i geotransformacji (parametry te zostały opisane w innej części kursu).
Cała procedura zostanie zamknięta w postaci funkcji. Funkcja przyjmuje cztery argumenty: dwa z nich, seria danych oraz raster referencyjny są obowiązkowe, dwa kolejne: nazwa pliku oraz wartość nodata są opcjonalne. W przypadku braku nazwy zostanie ona pobrana z nazwy serii danych (1), natomiast wartość nodata - dla danych zmiennoprzecinkowych przyjmuje wartość -9999. Jako seria danych może zostać przekazana zarówno kolumna dataFrame jak i jednowymiarowa tablica numpy. W pierwszej kolejności funkcja sprawdza czy dane są typu kategoryzowanego i w takiej sytuacji nastąpi zamiana na typ uint8 (2). Następnie funkcja sprawdzi czy dostarczona seria jest typu seria pandas czy numpy array (3), i w przypadku gdy jest to seria danych zostanie zamieniona na numpy array. Następnie jeżeli w danych występują wartości nan (nie liczby) i zamienione są na wartości nodata (4).
Utworzenie pliku geoprzestrznnego wymaga podanie typu danych, który podawany jest w formie stałej GDT_. W tym celu korzystamy z funkcji NumericTypeCodeToGDALTypeCode z modułu gdal_array, która pozwala na określenie typu danych gdal na podstawie typu danych numerycznych tablicy numpy (5). Po określeniu typu danych pozostała część funkcji to standardowa procedura utworzenia nowego pliku danych, gdzie parametry rastra (za wyjątkiem typu danych) tworzona jest na podstawie rastra referencyjnego.
from osgeo import gdal_array
def seriesToRaster(series,raster,newfile=None,nodata=-9999):
if not newfile:
newfile = series.name + '.tif' #1
if series.dtype.name== 'object':
return False
if series.dtype.name == 'category': #2
nodata=0
series.replace(np.NaN,nodata,inplace=True)
series = series.astype('uint8')
if type(series) == pd.core.series.Series: #3
series = series.values
series[np.where(np.isnan(series))]=nodata #4
gdt_Type = gdal_array.NumericTypeCodeToGDALTypeCode(series.dtype) #5
target = gdal.GetDriverByName('GTiff').Create(newfile, raster.RasterXSize, raster.RasterYSize, 1, gdt_Type) #6
target.SetProjection(raster.GetProjection())
target.SetGeoTransform(raster.GetGeoTransform())
tband=target.GetRasterBand(1)
series.shape=(raster.RasterYSize, raster.RasterXSize)
tband.WriteArray(series)
tband.FlushCache()
tband.SetNoDataValue(nodata)
target=None
return True
Jako przykład dodamy kolumnę zawierającą dane losowe na obszarze zajmowanym przez granice miasta Poznania. Przygotowanie tej kolumny wymaga kilku narzędzi pandas, z którymi zapoznamy się w dalszej części kursu
complete = df['20_dzielnice'].notna()
df['wynik'] = 0
df['wynik'][complete] = np.random.randint(1,10,size=sum(complete))
df['wynik'] = df['wynik'].astype('category')
df['wynik'].replace(0,np.NaN,inplace=True)
# funkcja wykonująca konwersję
seriesToRaster(df['wynik'],raster,'dzielnice.tif')
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/programowanie4/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)
Aby wyświetlić histogram zgodny z kolorami danych rastrowych mapy należy wykorzystać wcześniej utworzoną paletę kolorów cols, która jest po prostu listą kolorów.
import seaborn as sns
sns.countplot(x="10_landuse",data=df.sample(10000),palette=cols)
