In [1]:
# setup
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
import warnings
warnings.filterwarnings('ignore')

Regresja i klasyfikacja dużych zbiorów danych rastrowych

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.

In [2]:
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.

In [3]:
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
Out[3]:
Index(['OGC_FID', '04_terange', '06_highbui', '07_buildup', '08_road_de',
       '05_imperm', '02_green', '01_roughne', '03_distwat', 'landuse',
       'tempr'],
      dtype='object')

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:

  1. sortujemy alfabetycznie nazwy plików,
  2. skracamy nazwy do 10 znaków
  3. tworzymy słownik zmian nazw o nazwie rename - jako klucz stara nazwa jako value nowa
  4. Wykonujemy metodę rename dla data frame
In [4]:
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
Out[4]:
['01_roughness.tif',
 '02_green.tif',
 '03_distwaters.tif',
 '04_terange.tif',
 '05_imperm.tif',
 '06_highbuildup.tif',
 '07_buildup.tif',
 '08_road_dens.tif']
Out[4]:
{'01_roughne': '01_roughness',
 '02_green': '02_green',
 '03_distwat': '03_distwaters',
 '04_terange': '04_terange',
 '05_imperm': '05_imperm',
 '06_highbui': '06_highbuildup',
 '07_buildup': '07_buildup',
 '08_road_de': '08_road_dens'}

Regresja i predykcja do pliku geoprzestrzennego

Dane uczące (zmienne niezależne) wymagają uprzedniego przeskalowania to wartości standaryzowanych

In [5]:
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)
In [6]:
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.

In [7]:
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))
Out[7]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f96b9768588>
Out[7]:
[(8, 16)]
Out[7]:
[(8, 16)]
In [8]:
from sklearn.metrics import mean_absolute_error
mean_absolute_error(y_train,values_train)
mean_absolute_error(y_test,values_test)
Out[8]:
0.94896230951067051
Out[8]:
1.0532537561816631

Pobór danych rastrowych i zapis do pliku rastrowego

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.

In [9]:
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.

  1. W pierwszej kolejności tworzymy pustą macierz o wymiarach: ilość pliów $ \times $ ilość komórek
  2. do każdego wiersza wczytujemy spłaszczoną tablicę z wartościami rastra
  3. Dokonujemy transpozycji macierzy aby obecne wiersze były kolumnami, tak jak w przypadku danych uczących.
In [10]:
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

Operacje na niekompletnych danych (raster zawierający NULL)

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.

In [11]:
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.

In [12]:
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

In [13]:
results.shape = (y,x)
masked_results = np.ma.masked_where(results == -9999, results)
plt.imshow(masked_results)
plt.colorbar()
Out[13]:
<matplotlib.image.AxesImage at 0x7f96b8a33048>
Out[13]:
<matplotlib.colorbar.Colorbar at 0x7f96b89ec828>

Ostatnim krokiem jest procedura zapisu do pliku.

  1. Tworzymy nowy geotiff o wymiarach takich jak plik referencyjny i typie danych Float32, ustawiamy projekcję i geotransformaty
  2. Pobieramy warstwę rastra
  3. Zapisujemy do niej macierz results
  4. zamykamy połączenie do pliku i ustawiamy -9999 jako wartość pustą dla pliku
In [14]:
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
Out[14]:
0
Out[14]:
0
Out[14]:
0
Out[14]:
0

Klasyfikacja - wiele klas (multiclass classification)

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.

In [15]:
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.

In [16]:
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)
Out[16]:
GridSearchCV(cv=5, error_score='raise',
       estimator=SVC(C=1.0, cache_size=200, class_weight=None, coef0=0.0,
  decision_function_shape='ovr', degree=3, gamma='auto', kernel='rbf',
  max_iter=-1, probability=False, random_state=None, shrinking=True,
  tol=0.001, verbose=False),
       fit_params=None, iid=True, n_jobs=1,
       param_grid={'kernel': ['rbf'], 'gamma': [0.2, 0.12, 0.05, 0.01], 'C': [5, 8, 12, 15]},
       pre_dispatch='2*n_jobs', refit=True, return_train_score='warn',
       scoring='accuracy', verbose=0)
Out[16]:
SVC(C=8, cache_size=200, class_weight=None, coef0=0.0,
  decision_function_shape='ovr', degree=3, gamma=0.12, kernel='rbf',
  max_iter=-1, probability=False, random_state=None, shrinking=True,
  tol=0.001, verbose=False)
Out[16]:
array([[120,  13,   0,   0,   1,   3],
       [ 13,  56,   0,   0,   0,   2],
       [  4,   4,   7,   0,   0,   2],
       [  0,   0,   0,   3,   0,   0],
       [  8,   2,   1,   0,  12,   3],
       [  8,   2,   0,   0,   0,  52]])
Out[16]:
array([[60, 13,  1,  0,  0,  2],
       [ 3, 30,  0,  0,  0,  0],
       [ 5,  3,  1,  0,  0,  2],
       [ 0,  0,  0,  0,  0,  1],
       [ 5,  0,  0,  0,  2,  4],
       [ 7,  0,  0,  0,  1, 16]])

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.

In [17]:
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
Out[17]:
array(['AGRI', 'FOREST', 'GREEN', 'HIGH', 'IND', 'URBAN'], dtype=object)
In [18]:
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)
Out[18]:
<matplotlib.image.AxesImage at 0x7f96b8985748>

Ostatnim krokiem jest zapisanie wyników jako zbioru rastrowego, pamiętając, że plik wynikowy powinien być typu Byte.

In [19]:
# 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
Out[19]:
0
Out[19]:
0
Out[19]:
0
Out[19]:
0

Identyfikacja i selekcja istotnych cech

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.

In [20]:
files
Out[20]:
['01_roughness.tif',
 '02_green.tif',
 '03_distwaters.tif',
 '04_terange.tif',
 '05_imperm.tif',
 '06_highbuildup.tif',
 '07_buildup.tif',
 '08_road_dens.tif']

Na nowo zostaną wprzeskalowane zmienne.

In [21]:
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.

In [22]:
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)
Out[22]:
array([ 0.17124148,  0.30713359,  0.01170216,  0.05853935,  0.35431617,
        0.20846526,  0.32283076,  0.20385351])
Out[22]:
array([  35.2492858 ,   61.95509705,    3.30279401,    8.71650884,
        130.53064217,   47.3626166 ,  100.87975363,   66.88850714])
Out[22]:
array([ 29.86472695,  48.73755128,   2.2155101 ,   7.18862933,
        85.53027743,  38.85206564,  71.21829884,  51.87576569])
In [23]:
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

In [24]:
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)
Out[24]:
array([ 0.07074411,  0.10881192,  0.01259474,  0.02713051,  0.29421483,
        0.14886997,  0.2444255 ,  0.052917  ])
Out[24]:
array([  2.34343956e+01,   3.99276230e+01,   1.85985538e-02,
         7.40291975e+00,   1.79560968e+02,   7.98162259e+01,
         1.59815740e+02,   1.08844027e+02])
Out[24]:
array([  5.69345159,   9.04800198,   0.04982484,   2.16288445,
        32.00495736,  16.45922383,  29.19876939,  21.36611803])

Klasyfikacje nienadzorowane w dużych zbiorach danych geoprzestrzennych (raster)

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.

In [25]:
import os
import numpy as np
import scipy as sc
from osgeo import gdal
os.chdir("/home/jarekj/Dropbox/dydaktyka/programowanie4/urbanData")
In [26]:
# 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.

In [27]:
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
Out[27]:
[0, 100, 200]
Out[27]:
[100, 100, 100]
Out[27]:
[0, 100, 200]
Out[27]:
[100, 100, 70]

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:

  1. Tworzymy pustą macierz numpy o liczbie kolumn równej ilości komórek w próbie i liczbie wierszy równej liczbie rastrów
  2. Dla każdego rastra w liście otwieramy, wczytujemy jako array a następnie spłaszczamy
  3. Dokonujemy transpozycji danych (są poukładane w wierszach)
  4. Określamy wielkość próbki: wartość całkowitą iloczyn liczby komórek i frakcji
  5. Generujemy zestaw liczb losowych z przedziału od 0 do liczby komórek o wyliczonej poprzednio wielkości
  6. Pobieramy wybrane wiersze
In [28]:
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.

  1. Zbiór taki trzeba zainicjować, tworzymy pustą tablicę numpy array o wymiarach 0 wierszy i 8 kolumn
  2. Następnie w podwójnej pętli dla każdego wiersza i kolumny w wierszu...
  3. ...pobieramy podpróbkę wg poznanej procedury
  4. Na ostatnim etapie dodajemy pod-próbę do tworzonej próby losowej przy pomocy funkcji row_stack()
In [29]:
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ąć:

  1. Wyznaczyć wektor logiczny complete, zawierający True jeżeli wiersz zawiera wszystkie dane i False jeżeli choćby jedna dana była pusta (any()).
  2. Wybrać tylko te wiersze które zawierają kompletne dane poprzez selekcję logiczną
In [30]:
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()

In [31]:
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,

In [32]:
klasy = gdal.GetDriverByName('GTiff').Create("1000_klasyfikacja.tif", x,y, 1, gdal.GDT_Byte)
klasy.SetProjection(proj)
klasy.SetGeoTransform(gt)
kband=klasy.GetRasterBand(1)
Out[32]:
0
Out[32]:
0

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:

  1. Dla każdego kroku tworzymy wektor complete
  2. Tworzymy pusty wektor wynikowy (classes). Nie zawiera on wartości nan ale zera dla miejsc, gdzie nie określimy klasy
  3. Wykonujemy klasyfikację przy pomocy potoku tylko na danych kompletnych a wynik przypisujemy tylko do tych miejsc w waktorze wynikowym gdzie dane są kompletne. W ten sposób omijamy wartości puste.
  4. Zmieniamy wymiary macierzy. Uwaga wymiary macierzy podawane są odwrotnie: liczba kolumn - liczba wierszy
  5. Zapisujemy macierz z odpowiednim przesunięciem
In [33]:
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)
Out[33]:
0
Out[33]:
0
Out[33]:
0
Out[33]:
0
Out[33]:
0
Out[33]:
0
Out[33]:
0
Out[33]:
0
Out[33]:
0

Ostanie kroki to zamknięcie rastra wynikowego i zdefiniowanie wartości pustej.

In [34]:
kband.FlushCache()
kband.SetNoDataValue(0)
del kband
klasy = None
Out[34]:
0