# setup
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
import warnings
warnings.filterwarnings('ignore')
W przypadku dużych zbiorów danych rastrowych stosowanie modeli predykcyjnych, zarówno nadzorowanych jak i nie nadzorowanych posiada wiele ograniczeń. Przede wszystkim dane posiadają strukturę uporządkowanej macierzy z nieregularnie rozłożonymi wartościami pustymi, które trzeba usunąć przed dokonaniem predykcji.
Praca z dużymi zbiorami danych rastrowych niesie za sobą pewien problem: dane są zbyt duże aby je jednorazowo wczytać do pamięci - w takiej sytuacji musimy dane przetwarzać porcjami (chunk), zarówno na etapie próbkowania jak i predykcji.
Do realizacji tego zadania użyjemy biblioteki GDAL oraz metody ReadAsArray()
, która pozwala na czytanie do pamięci wybranych fragmentów zbioru rastrowego. Ćwiczenie zostanie wykonane na małym zbiorze, który zostanie podzielony w sposób sztuczny jedynie po to aby zilustrować cały proces. Zbiory danych o rozmiarze do kilku milionów komórek możemy z powodzeniem przetwarzać w jednym kroku. Zbiór danych liczy 8 warstw rastrowych - zmiennych modelu - ilustrujących wybrane cechy pokrycia terenu miasta Poznania. Do zbudowania procesu przetwarzania w porcjach musimy zdefiniować wielkość kroku (ilość komórek rastra) o jaka jest doczytywana w trakcie predycji na nowych danych.
import os
import sqlite3
import numpy as np
import scipy as sc
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from osgeo import gdal
os.chdir("/home/jarekj/Dropbox/dydaktyka/programowanie4/urbanData")
Dane uczące najczęściej znajdują się w osobnych zbiorach wektorowych, dla których należy wykonać próbkowanie w miejscach gdzie wartość zmiennej zależnej jest znana przy pomocy dowolnego systemu GIS. W poniższym ćwiczeniu użyjemy danych już opróbowanych zapisanych w pliku spatialite. Dane te będą niezbędne na etapie budowania modelu. Dane rastrowe będą stanowiły nowe dane, dla których zostanie wykonana predykcja.
database = 'proby.sqlite' #1A
table = 'proby' #1B
conn=sqlite3.connect(database) #2A
cur = conn.cursor() #2B
query = "SELECT * FROM {0};".format(table)#3
data = pd.read_sql_query(query,conn)#4
data = data.drop("GEOMETRY",axis=1)
data.columns
Kolejnym elementem pracy z danymi jest synchronizacja nazw kolumn danych wektorowych przechowujących zbiór uczący z nazwami plików danych rastrowych. W naszym przypadku nazwy kolumn są skrócone do 10 znaków (wymogi shapefile i niektórych innych baz danych) w związku z tym musimy zmodyfikować nazwy kolumn data Frame na takie jakie pochodzą z nazw plików rastrowych. W tym celu:
files = os.listdir(".")
files = [file for file in files if file.startswith("0") and file.endswith("tif")]
files.sort() #1
files
current_colnames = [file.split(".")[0][0:10] for file in files] #2
colnames = [file.split(".")[0] for file in files]
rename = dict(zip(current_colnames,colnames)) #3
rename
data.rename(columns = rename,inplace=True) #4
Dane uczące (zmienne niezależne) wymagają uprzedniego przeskalowania to wartości standaryzowanych
from sklearn.preprocessing import StandardScaler
X = data.loc[~pd.isnull(data.tempr),colnames]
yg = data.tempr[~pd.isnull(data.tempr)]
scl = StandardScaler().fit(X) # zachować fit
scaled = scl.transform(X)
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(scaled, yg, test_size=0.33)
To oceny jakości modelu możemy użyć prostego wykresu punktowego lub metryk dla modeli regresyjnych.
from sklearn.svm import SVR
model = SVR(C=30,gamma=0.001).fit(X_train,y_train)
values_train = model.predict(X_train)
values_test = model.predict(X_test)
fig, ax = plt.subplots()
g = sns.scatterplot(y_train,values_train,ax=ax)
sns.scatterplot(y_test,values_test,ax=ax)
g.set(ylim=(8,16))
g.set(xlim=(8,16))
from sklearn.metrics import mean_absolute_error
mean_absolute_error(y_train,values_train)
mean_absolute_error(y_test,values_test)
proces predykcji na danych rastrowych wymaga wcześniejszego wczytania danych do środowiska obliczeniowego. Ponieważ wyniki będą zapisywane również jako plik rastrowy rozpoczniemy operację od poboru pliku referencyjnego - pierwszego rastra z listy w celu pobrania rozmiaru, projekcji i geotransforamcji.
ref = gdal.Open(files[0])
x = ref.RasterXSize
y = ref.RasterYSize
proj = ref.GetProjection()
gt = ref.GetGeoTransform()
Następnie wszystkie dane rastrowe pobierzemy do jednej macierzy numpy.
rastdata=np.empty([len(files),x*y])
for i in range(len(files)):
raster = gdal.Open(files[i])
rastdata[i] = raster.GetRasterBand(1).ReadAsArray().flatten()
raster=None
rastdata = rastdata.T
Ponieważ dane zawierają wartości puste musimy je usunąć. W tym celu tworzymy wektor logiczny (bool) o długości równej ilości komórek, gdzie w miejscu gdzie jest komplet wartości jest True, a tam gdzie ich brakuje - False. Służy do tego funkcja isnan()
Any(axis=1)
oznacza że ma ona być zastosowana dla każdego wiersza tablicy.
Gdy mamy już wektor complete używamy go jako selektora aby wybrać z oryginalnej tablicy tylko kompletne dane. Następnie skalujemy dane wykorzystując model zbudowany na próbce, ponieważ dane muszą być skalowane wg tych samych kryteriów.
complete = ~np.isnan(rastdata).any(axis=1)
complete_rastdata = rastdata[complete]
rastscaled = scl.transform(complete_rastdata)
Ponownie budujemy model, tym razem na całości danych treningowych (podział na zbiór testowy i treningowy służy jedynie do oceny wydajności regresora), a predykcję wykonujemy na danych rastrowych z usuniętymi komórkami NULL.
W następnym kroku tworzymy pustą tablicę wynikową, zawierającą wartości "-9999", które potraktujemy jako symbol wartości pustych o długości równej ilości komórek oryginalnego rastra. Do pliku wynikowego tylko dla miejsc gdzie były kompletne dane przypisujemy wyniki predykcji.
model = SVR(C=30,gamma=0.001).fit(scaled,yg)
prediction = model.predict(rastscaled)
results = np.repeat(-9999.,len(complete)) # NA, musi być z przecinkiem lub wymusić typ danych
results[complete] = prediction
Dokonujemy konwersji wektora do oryginalngeo kształu macierzy, gdzie y do ilość wierszy a x kolumn. Macierz możemy podejrzeć, maskując wartości -9999
results.shape = (y,x)
masked_results = np.ma.masked_where(results == -9999, results)
plt.imshow(masked_results)
plt.colorbar()
Ostatnim krokiem jest procedura zapisu do pliku.
temp = gdal.GetDriverByName('GTiff').Create("1000_temperatura.tif", x,y, 1, gdal.GDT_Float32)
temp.SetProjection(proj)
temp.SetGeoTransform(gt)
tband=temp.GetRasterBand(1)
tband.WriteArray(results) # shape nadany przed wizualizacją
tband.FlushCache()
tband.SetNoDataValue(-9999)
del tband
temp = None
Dla tego samego zbioru danych geoprzestrzennych wykonamy klasyfikację, gdzie zmienna zależna posiada kilka klas - w tym wypadku klas pokrycia terenu. Do wykoania klasyfikacji będzie potrzebna LabelEncoder - klasa zamieniająca klasy kodowane jako łańcuchy tekstowe na wartości liczbowe. Klasa działa tak samo jak inne klasy scikit-learn, stosując metody fit()
i transform()
. Zmienne niezależne ponownie muszą zostać przeskalowane do wartości standaryzowanych.
from sklearn.svm import SVC
from sklearn.preprocessing import LabelEncoder
from sklearn.metrics import confusion_matrix
from sklearn.model_selection import GridSearchCV
X = data.loc[:,colnames]
yg = data.landuse
scl = StandardScaler().fit(X) # zachować fit
scaled = scl.transform(X)
encoder = LabelEncoder().fit(yg)
yl = encoder.transform(yg)
X_train, X_test, y_train, y_test = train_test_split(scaled, yl, test_size=0.33)
Klasyfikacje gdy zmienna zależna reprezentuje wiele klas z reguły nie działają dobrze na wartościach domyślnych klasyfikatora. Z tego powodu tuning parametrów jest niezbędnym krokiem w klasyfikacji.
tuned_parameters = {'kernel': ['rbf'], 'gamma': [0.2, 0.12, 0.05, 0.01],'C': [5,8,12,15,]}
clf = GridSearchCV(SVC(),tuned_parameters,cv=5,scoring='accuracy')
clf.fit(scaled,yl)
clf.best_estimator_
model = SVC(gamma=0.11,C=3).fit(X_train,y_train)
classes_train = model.predict(X_train)
confusion_matrix(y_train,classes_train)
classes_test = model.predict(X_test)
confusion_matrix(y_test,classes_test)
Po oszacowaniu wydajności klasyfikatora budujemy nowy klasyfikator na całym zbiorze danych i stosujemy go do nowych danych. Należy pamiętać że wynik klasyfikacji to wektor wartości zaczynających się od 0. Można go zamienić na oryginalne etykiety stosując metodę inverse_transform()
, natomiast w celu zapisu jako raster z kategoriami należy dodać wartość 1.
model = SVC(gamma=0.11,C=3).fit(scaled,yl)
rastscaled = scl.transform(complete_rastdata) # complete już istnieje ale mamy inną transformację
classes = model.predict(rastscaled)
labels = encoder.inverse_transform(classes)
np.unique(labels)
classes+=1
results = np.repeat(0,len(complete))
results[complete] = classes
results.shape = (y,x)
masked_results = np.ma.masked_where(results == 0, results)
plt.imshow(masked_results)
Ostatnim krokiem jest zapisanie wyników jako zbioru rastrowego, pamiętając, że plik wynikowy powinien być typu Byte.
# raster referencyjny już jest
clasf = gdal.GetDriverByName('GTiff').Create("1000_landcover.tif", x,y, 1, gdal.GDT_Byte)
clasf.SetProjection(proj)
clasf.SetGeoTransform(gt)
cband=clasf.GetRasterBand(1)
cband.WriteArray(results)
cband.FlushCache()
cband.SetNoDataValue(0)
del cband
clasf = None
Zmienne niezależne mogą być w różnej relacji do zmiennych zależnych. Wiele może po prostu być nieistotnych - tj nie mieć żadnego związku z wartościami/klasami zmiennych zależnych. W przypadku wielu klasyfikatorów warto rozważyć wcześniejsze wyeliminowanie takich zmiennych. Pakiet scikit-learn dostarcza kilka narzędzi pozwalających na ocenę znaczenia zmiennych. Mogą to być proste testy analizujące związki każdej ze zmiennych osobno ze zmienną zależną lub nieco bardziej skomplikowane metody (nie wszystkie działają) trenujące klasyfikatory/regresrory poprzez usuwanie kolejnych zmiennych.
files
Na nowo zostaną wprzeskalowane zmienne.
X = data.loc[:,colnames]
yg = data.landuse
scl = StandardScaler().fit(X)
scaled = scl.transform(X)
encoder = LabelEncoder().fit(yg)
yl = encoder.transform(yg)
W poniższym przykładzie wykorzystamy dwie metryki: informację wzajemną oraz test f, zarówno dla regresji jak i klasyfikacji. Dla mutual_info_
ważność zmiennych określa udział (od 0 do 1) wspólnej informacji która jest zawarta zarówno w zmiennej zależnej jak i każdej ze zmiennych niezależnych. W przykładzie widać, że zmienne: 3 i 4 (odległość od wody i różnice wysokości) mają najmniejszy wpływ zarówno na klasę terenu natomiast najistotniejsze są nieprzepuszczalność podłoża udział zieleni oraz zabudowa. Te same zmienne wskazuje statystyka F.
from sklearn.feature_selection import mutual_info_classif, f_classif, SelectPercentile
m = mutual_info_classif(scaled,yl)
m
f,p = f_classif(scaled,yl)
f
-np.log10(p)
X = data.loc[~pd.isnull(data.tempr),colnames]
yg = data.tempr[~pd.isnull(data.tempr)]
scl = StandardScaler().fit(X) # zachować fit
scaled = scl.transform(X)
Dla regresora - wyznaczającego temperaturę powierzchni nieistotne są wszystkie zmienne za wyjątkiem nieprzepuszczalności i zabudowy, z tym że statystyka F pokazuje jeszcze pewną istotność gęstości dróg
from sklearn.feature_selection import mutual_info_regression, f_regression
m = mutual_info_regression(scaled,yg)
m
f,p = f_regression(scaled,yg)
f
-np.log10(p)
W przypadku klasyfikacji nienadzorowanych danych rastrowych pojawia się dodatkowy problem: dane są zbyt duże aby je sklasyfikować - ten problem zwykle rozwiązujemy poprzez próbkowanie zbioru i trening klasyfikatora na próbie losowej, a następnie dokonujemy predykcji dla całego zbioru danych. Z tego powodu dla klasyfikacji nienadzorowanych - oprócz wielkości kroku rastra musimy zdefiniować wielkość próby losowej (ułamek) jaki chcemy użyć do zbudowania modelu klasyfikacyjnego.
import os
import numpy as np
import scipy as sc
from osgeo import gdal
os.chdir("/home/jarekj/Dropbox/dydaktyka/programowanie4/urbanData")
# pobieramy wielkość rastra wzorcowego
# lista rastrów
files = os.listdir(".")
files = [file for file in files if file.startswith("0") and file.endswith("tif")]
files.sort()
# tworzymy krok wielkości submacierzy
step = 100
sample_size = 3/100.
Po zdefiniowaniu parametru i wskazaniu zbioru danych potrzebujemy raster referencyjny, z którego pobierzemy wymiary danych. Zakładamy, że wszystkie zbiory są tego samego rozmiaru i że zostało to sprawdzone na etapie ich przygotowania. Z rastra referencyjnego pobieramy wymiary, projekcję i geotransformaty a następnie wyznaczamy 4 listy: listę przesunięć (offsets) dla rastra są to po prostu lista kolejnych wartości oddzielonych od siebie o krok, oraz lista wielkości okna, która jest po prostu powtórzeniem wartości kroku tyle razy ile mamy kroków w danym wymiarze. Ostatni krok z reguły jest mniejszy i z tego powodu musimy go dodać osobno. Ilość powtórzeń wyznaczymy jako wynik dzielenia całkowitego ilości komórek przez wielkość kroków a ostatni krok to wynik operacji modulo (reszta z dzielenia). Jeżeli jest ona różna od zera dodajemy ją na końcu listy.
ref = gdal.Open(files[0])
x = ref.RasterXSize
y = ref.RasterYSize
proj = ref.GetProjection()
gt = ref.GetGeoTransform()
xoffs = list(range(0,x,100))
yoffs = list(range(0,y,100))
# wyliczenie zakresów
ysizes = [step]*(y//step)
if y%step: ysizes.append(y%step)
xsizes = [step]*(x//step)
if x%step: xsizes.append(x%step)
xoffs
xsizes
yoffs
ysizes
Następnie wykonamy procedurę próbkowania w porcjach. Zakładamy, że nie jesteśmy w stanie wczytać całego zbioru do pamięci i w związku z tym wczytujemy porcję określoną przesunięciem i wielkością okna i pobieramy z niej losową pod-próbkę. W tym celu:
xsize=100
ysize=100
xoff = 0
yoff = 0
data=np.empty([len(files),xsize*ysize]) #1
for i in range(len(files)): #2
raster = gdal.Open(files[i])
data[i] = raster.GetRasterBand(1).ReadAsArray(xoff,yoff,xsize,ysize).flatten()
raster=None
data = data.T #3
sample_num = int(xsize*ysize*sample_size) #4
sample = np.random.randint(0,xsize*ysize,sample_num) #5
subsample = data[sample,] #6
Jeżeli procedura działa dla jednej porcji, możemy ją powtórzyć dla wszystkich zbioru. Każdy krok będzie generował pod-próbę, którą będziemy dodawać do zbioru wynikowego.
row_stack()
sample_data = np.empty((0,8),float) #1
for xoff,xsize in zip(xoffs,xsizes): #2A
for yoff,ysize in zip(yoffs,ysizes): #2B
data=np.empty([len(files),xsize*ysize]) #3
for i in range(len(files)):
raster = gdal.Open(files[i])
data[i] = raster.GetRasterBand(1).ReadAsArray(xoff,yoff,xsize,ysize).flatten()
raster=None
data = data.T
sample_num = int(xsize*ysize*sample_size)
sample = np.random.randint(0,xsize*ysize,sample_num)
subsample = data[sample,]
sample_data = np.row_stack((sample_data,subsample)) #4
Nasze dane zawierają wartości puste typowe dla danych rastrowych. Należy je usunąć:
True
jeżeli wiersz zawiera wszystkie dane i False
jeżeli choćby jedna dana była pusta (any()
).complete = ~np.isnan(sample_data).any(1)
complete_sample_data = sample_data[complete,:]
Ponieważ dla każdej porcji będziemy powtarzali tę samą procedurę, warto skorzystać z poznanego wcześniej przetwarzania potokowego. Użyjemy tego samego zestawu funkcji co poprzednio co poprzednio. Po zbudowaniu potoku możemy dokonać predykcji na próbie danych i ocenić jakość klasyfikacji poprzez predict_proba()
from sklearn.pipeline import Pipeline
from sklearn import decomposition
from sklearn import mixture
from sklearn import preprocessing
scale = preprocessing.StandardScaler()
gmm = mixture.GaussianMixture(n_components=5, covariance_type='full')
pca = decomposition.PCA(n_components=3)
pipe = Pipeline(steps=[('scale',scale),('pca',pca),('gmm',gmm)])
model = pipe.fit(complete_sample_data)
base_classes = model.predict(complete_sample_data)+1
base_proba = model.predict_proba(complete_sample_data)
Ostatnim krokiem jest klasyfikacja i zapis danych, które musimy wykonać równolegle. W pierwszym kroku utworzymy nowy raster typu Byte() według poznanej wcześniej procedury,
klasy = gdal.GetDriverByName('GTiff').Create("1000_klasyfikacja.tif", x,y, 1, gdal.GDT_Byte)
klasy.SetProjection(proj)
klasy.SetGeoTransform(gt)
kband=klasy.GetRasterBand(1)
a następnie powtarzamy procedurę zastosowaną przy czytaniu danych, z tą różnicą że dla każdej porcji dokonujemy klasyfikacji używając wcześniej zdefiniowany potok (w tym miejscu widzimy jakie to ułatwienie). Ponieważ jednak nasze dane zawierają wartości puste:
complete
for xoff,xsize in zip(xoffs,xsizes):
for yoff,ysize in zip(yoffs,ysizes):
data=np.empty([len(files),xsize*ysize])
for i in range(len(files)):
raster = gdal.Open(files[i])
data[i] = raster.GetRasterBand(1).ReadAsArray(xoff,yoff,xsize,ysize).flatten()
raster=None
data = data.T
#
complete = ~np.isnan(data).any(1)
classes = np.repeat(0,len(complete))
classes[complete] = pipe.predict(data[complete,])+1
classes.shape = (ysize,xsize)
kband.WriteArray(classes,xoff,yoff)
Ostanie kroki to zamknięcie rastra wynikowego i zdefiniowanie wartości pustej.
kband.FlushCache()
kband.SetNoDataValue(0)
del kband
klasy = None