12.1 Land Cover klasifikace s využitím Sentinel-2 a ST_LUCAS¶

"Template" pro řešení semestrálního projketu.

In [ ]:

import os 
import glob
import rasterio
import rasterio.features
from rasterio.enums import MergeAlg
from rasterio.plot import show
from skimage import exposure
import numpy as np
from shapely.geometry import box
import geopandas as gpd 
from st_lucas import LucasRequest
from st_lucas import LucasIO
import matplotlib.pyplot as plt 
import matplotlib.patches as mpatches

import os import glob import rasterio import rasterio.features from rasterio.enums import MergeAlg from rasterio.plot import show from skimage import exposure import numpy as np from shapely.geometry import box import geopandas as gpd from st_lucas import LucasRequest from st_lucas import LucasIO import matplotlib.pyplot as plt import matplotlib.patches as mpatches

1/ Popis úkolu¶

Rámcový popis problému: ...

Předpoklady úlohy: ...

Očekávaný výsledek: ...

2/ Průzkumná analýza dat¶

Sentinel-2¶

EPSG: 32633

In [ ]:

# Data Sentinel-2 
path = './data'
s2prod = 'S2*.SAFE/GRANULE/*/IMG_DATA/R20m/*B*.jp2'
ips = glob.glob(os.path.join(path, s2prod))
ips.sort()

# Data Sentinel-2 path = './data' s2prod = 'S2*.SAFE/GRANULE/*/IMG_DATA/R20m/*B*.jp2' ips = glob.glob(os.path.join(path, s2prod)) ips.sort()

In [ ]:

# bands = ['B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B07', 'B11', 'B12', 'B8A'] 
bands = []
for ip in ips: 
    print(os.path.basename(ip))
    ds = rasterio.open(ip) 
    array = ds.read(1)
    bands.append(array)
array = None 
print('done!')

# bands = ['B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B07', 'B11', 'B12', 'B8A'] bands = [] for ip in ips: print(os.path.basename(ip)) ds = rasterio.open(ip) array = ds.read(1) bands.append(array) array = None print('done!')

In [ ]:

ds.crs

ds.crs

In [ ]:

# Raster extent 
bounds  = ds.bounds
geom = box(*bounds)
# print(geom.wkt)

# Raster extent bounds = ds.bounds geom = box(*bounds) # print(geom.wkt)

In [ ]:

raster_bounds = gpd.GeoDataFrame({"id":1,"geometry":[box(*bounds)]})
raster_bounds.set_crs(32633, inplace=True)
# raster_bounds.crs

raster_bounds = gpd.GeoDataFrame({"id":1,"geometry":[box(*bounds)]}) raster_bounds.set_crs(32633, inplace=True) # raster_bounds.crs

In [ ]:

s2_array = np.array(bands)

s2_array = np.array(bands)

In [ ]:

p2, p98 = np.percentile(bands[2], (2,98))
image = exposure.rescale_intensity(bands[2], in_range=(p2, p98)) / 100000

show(image, transform=ds.transform, cmap='pink')

p2, p98 = np.percentile(bands[2], (2,98)) image = exposure.rescale_intensity(bands[2], in_range=(p2, p98)) / 100000 show(image, transform=ds.transform, cmap='pink')

ST_LUCAS¶

EPSG: 3035

In [ ]:

# Data - ST_LUCAS
request = LucasRequest()
lucasio = LucasIO()

# Data - ST_LUCAS request = LucasRequest() lucasio = LucasIO()

In [ ]:

# AOI filter
# request.bbox = (4504276, 3020369, 4689608, 3105290) 
request.countries = ['CZ']

# AOI filter # request.bbox = (4504276, 3020369, 4689608, 3105290) request.countries = ['CZ']

In [ ]:

# Temporal filter 
request.years = [2018]

# Temporal filter request.years = [2018]

In [ ]:

# Thematic filter 
request.group = 'LC_LU' 

# Thematic filter request.group = 'LC_LU'

In [ ]:

# Download LUCAS points based on request 
lucasio.download(request)

# Download LUCAS points based on request lucasio.download(request)

In [ ]:

# Geopandas
lucas_df = lucasio.to_geopandas()
# lucas_df.head(3)

# Geopandas lucas_df = lucasio.to_geopandas() # lucas_df.head(3)

In [ ]:

# reprojekce
lucas_df_32633 = lucas_df.to_crs(32633)

# reprojekce lucas_df_32633 = lucas_df.to_crs(32633)

Prostorový výběr LUCAS bodů¶

In [ ]:

f, ax = plt.subplots()
raster_bounds.plot(color='white', edgecolor='black', ax=ax, zorder=1)
lucas_df_32633.plot(marker='o', color='red', markersize=0.5, ax=ax, zorder=2) 

f, ax = plt.subplots() raster_bounds.plot(color='white', edgecolor='black', ax=ax, zorder=1) lucas_df_32633.plot(marker='o', color='red', markersize=0.5, ax=ax, zorder=2)

In [ ]:

# výběr 
lucas_df_32633_sel = lucas_df_32633.overlay(raster_bounds, how='intersection')
lucas_df_32633_sel.plot(alpha=0.5, edgecolor='k', cmap='tab10');

# výběr lucas_df_32633_sel = lucas_df_32633.overlay(raster_bounds, how='intersection') lucas_df_32633_sel.plot(alpha=0.5, edgecolor='k', cmap='tab10');

In [ ]:

# lucas_df_32633_sel.head(3)

# lucas_df_32633_sel.head(3)

In [ ]:

lucas_df_32633_sel.shape

lucas_df_32633_sel.shape

In [ ]:

# Translace LUCAS level-1 na land cover kody 
translace = {
  'A': 1, # umělé povrchy (zástavba) 
  'B': 2, # zemědělství 
  'C': 3, # les
  'D': 4, # křoviny
  'E': 5, # louky 
  'F': 6, # holá půda, skály 
  'G': 7, # voda
  'H': 8, # mokřady 
  'X': 9  # ostatní 
}

# Translace LUCAS level-1 na land cover kody translace = { 'A': 1, # umělé povrchy (zástavba) 'B': 2, # zemědělství 'C': 3, # les 'D': 4, # křoviny 'E': 5, # louky 'F': 6, # holá půda, skály 'G': 7, # voda 'H': 8, # mokřady 'X': 9 # ostatní }

In [ ]:

# level-1 
lucas_df_32633_sel['lc1_l1'] = lucas_df_32633_sel['lc1'].str[0:1]

# level-1 lucas_df_32633_sel['lc1_l1'] = lucas_df_32633_sel['lc1'].str[0:1]

In [ ]:

lucas_df_32633_sel['lc'] = lucas_df_32633_sel['lc1_l1'].map(translace)

lucas_df_32633_sel['lc'] = lucas_df_32633_sel['lc1_l1'].map(translace)

In [ ]:

kody = list(lucas_df_32633_sel['lc'].unique())
kody

kody = list(lucas_df_32633_sel['lc'].unique()) kody

Rasterizace reference¶

In [ ]:

lucas_df_32633_sel['buf'] = lucas_df_32633_sel.buffer(20)

lucas_df_32633_sel['buf'] = lucas_df_32633_sel.buffer(20)

In [ ]:

lucas_df_32633_sel['lc'].unique() 

lucas_df_32633_sel['lc'].unique()

In [ ]:

# geometrie s hodnotami 
geom_value = ((geom,value) for geom, value in zip(lucas_df_32633_sel['buf'], lucas_df_32633_sel['lc']))

# Rasterizace 
lucas_array = rasterio.features.rasterize(geom_value,
                                out_shape = s2_array[0, :, :].shape, 
                                transform = ds.transform,
                                all_touched = True,
                                fill = 0,   
                                merge_alg = MergeAlg.replace,
                                dtype = np.int16)

# geometrie s hodnotami geom_value = ((geom,value) for geom, value in zip(lucas_df_32633_sel['buf'], lucas_df_32633_sel['lc'])) # Rasterizace lucas_array = rasterio.features.rasterize(geom_value, out_shape = s2_array[0, :, :].shape, transform = ds.transform, all_touched = True, fill = 0, merge_alg = MergeAlg.replace, dtype = np.int16)

In [ ]:

# kntrola
np.unique(lucas_array)

# kntrola np.unique(lucas_array)

Výběr referencí¶

In [ ]:

# maska referencí 
lucas_array.shape

# maska referencí lucas_array.shape

In [ ]:

# příprava vstupních dat pod maskou referencí lucas_array
X = pass 

# příprava vstupních dat pod maskou referencí lucas_array X = pass

In [ ]:

# příprava LUCAS referencí pod maskou lucas_array
pass
y = pass 

# příprava LUCAS referencí pod maskou lucas_array pass y = pass

In [ ]:

print(f'Dim. vstupu (S2): {X.shape}')
print(f'Dim. referenci (LUCAS): {y.shape}')

print(f'Dim. vstupu (S2): {X.shape}') print(f'Dim. referenci (LUCAS): {y.shape}')

Vizualizace dat¶

In [ ]:

# vizualizace
# graf rozptylu 
# příznaky 
# t-SNE 

# vizualizace # graf rozptylu # příznaky # t-SNE

In [ ]:

# Definice tříd 
nomenklatura ={
        1: "zástavba",
        2: "zemědělství",
        3: "les",
        4: "křoviny",
        5: "louky",
        6: "holá půda", 
        7: "voda", 
        8: "mokřady", 
        9: "ostatní", 
}

# Definice tříd nomenklatura ={ 1: "zástavba", 2: "zemědělství", 3: "les", 4: "křoviny", 5: "louky", 6: "holá půda", 7: "voda", 8: "mokřady", 9: "ostatní", }

In [ ]:

barvy = {
    1: (255./255,   0./255,   0./255, 0.5),   # red - zástavba
    2: (255./255, 255./255,   0./255, 0.5),   # yeallow - zemědělství  
    3: ( 51./255, 160./255,  44./255, 1  ),   # dark green - les
    4: (223./255, 223./255, 138./255, 1  ),   # light greyish green - křoviny 
    5: (  0./255, 223./255,   0./255, 0.5),   # green - louky
    6: (211./255, 211./255, 211./255, 1  ),   # grey - holá půda
    7: (  0./255,   0./255, 255./255, 0.5),   # blue - voda 
    8: (238./255, 130./255, 238./255, 0.5),   # violet - mokřady 
    9: (255./255, 255./255, 255./255, 0.5)    # white - ostatní
}

# (255./255, 192./255, 180./255, 0.5) # pink 

barvy = { 1: (255./255, 0./255, 0./255, 0.5), # red - zástavba 2: (255./255, 255./255, 0./255, 0.5), # yeallow - zemědělství 3: ( 51./255, 160./255, 44./255, 1 ), # dark green - les 4: (223./255, 223./255, 138./255, 1 ), # light greyish green - křoviny 5: ( 0./255, 223./255, 0./255, 0.5), # green - louky 6: (211./255, 211./255, 211./255, 1 ), # grey - holá půda 7: ( 0./255, 0./255, 255./255, 0.5), # blue - voda 8: (238./255, 130./255, 238./255, 0.5), # violet - mokřady 9: (255./255, 255./255, 255./255, 0.5) # white - ostatní } # (255./255, 192./255, 180./255, 0.5) # pink

In [ ]:

# Viz. 
plt.subplot(131)
plt.imshow(s2_array[3, :, :], cmap=plt.cm.Greys_r)
plt.title('RED')

plt.subplot(132)
plt.imshow(s2_array[4, :, :], cmap=plt.cm.Greys_r)
plt.title('NIR')

plt.show()

# Viz. plt.subplot(131) plt.imshow(s2_array[3, :, :], cmap=plt.cm.Greys_r) plt.title('RED') plt.subplot(132) plt.imshow(s2_array[4, :, :], cmap=plt.cm.Greys_r) plt.title('NIR') plt.show()

In [ ]:

# no 0 
lucas_array2 = lucas_array[~np.all(lucas_array == 0, axis=1)]
idx = np.argwhere(np.all(lucas_array2[:, ...] == 0, axis=0))
zmensene_y = np.delete(lucas_array2, idx, axis=1)

# no 0 lucas_array2 = lucas_array[~np.all(lucas_array == 0, axis=1)] idx = np.argwhere(np.all(lucas_array2[:, ...] == 0, axis=0)) zmensene_y = np.delete(lucas_array2, idx, axis=1)

In [ ]:

# barvy v matici 
barevne_pole = np.zeros((zmensene_y.shape[0], zmensene_y.shape[1], 4)) #rgba
tot_size = dict()
for label, color in barvy.items():
    mask = class_labels == label
    tot_size[label] = np.sum(mask==1)
    barevne_pole[mask] = color

# barvy v matici barevne_pole = np.zeros((zmensene_y.shape[0], zmensene_y.shape[1], 4)) #rgba tot_size = dict() for label, color in barvy.items(): mask = class_labels == label tot_size[label] = np.sum(mask==1) barevne_pole[mask] = color

In [ ]:

# Zmensene pole LUCAS bodu s barvou 
plt.figure(figsize = (15,15))
plt.imshow(barevne_pole, interpolation='none')
patches = [mpatches.Patch(color=barvy[key], label= nomenklatura[key]) for key in barvy.keys() ]
plt.legend(handles=patches, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

# Zmensene pole LUCAS bodu s barvou plt.figure(figsize = (15,15)) plt.imshow(barevne_pole, interpolation='none') patches = [mpatches.Patch(color=barvy[key], label= nomenklatura[key]) for key in barvy.keys() ] plt.legend(handles=patches, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

In [ ]:

# Spojita barevna tabulka 
import matplotlib as mpl
fig, ax = plt.subplots(figsize=(6, 1))
cmap = (mpl.colors.ListedColormap(['red', 'yellow', 'darkgreen', 'olive', 'lightgreen', 'grey', 'blue', 'violet', 'white'])
        .with_extremes(under='black', over='black'))
bounds = [1, 2, 3, 4, 5, 6, 7, 8, 9]
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

fig.colorbar(
    mpl.cm.ScalarMappable(cmap=cmap, norm=norm),
    cax=ax, orientation='horizontal',
    extend='both',
    spacing='proportional',
    label='Tridy',
)

# Spojita barevna tabulka import matplotlib as mpl fig, ax = plt.subplots(figsize=(6, 1)) cmap = (mpl.colors.ListedColormap(['red', 'yellow', 'darkgreen', 'olive', 'lightgreen', 'grey', 'blue', 'violet', 'white']) .with_extremes(under='black', over='black')) bounds = [1, 2, 3, 4, 5, 6, 7, 8, 9] norm = mpl.colors.BoundaryNorm(bounds, cmap.N) fig.colorbar( mpl.cm.ScalarMappable(cmap=cmap, norm=norm), cax=ax, orientation='horizontal', extend='both', spacing='proportional', label='Tridy', )

3/ Příprava dat¶

In [ ]:

# Příprava příznaků (feature engineering) 
# Chybějící data
# Harmonizace dat 
# Rozdělení na trénovací a testovací 

# Příprava příznaků (feature engineering) # Chybějící data # Harmonizace dat # Rozdělení na trénovací a testovací

In [ ]:

from sklearn.model_selection import train_test_split

from sklearn.model_selection import train_test_split

In [ ]:

X_train, X_test, y_train, y_test = pass 

X_train, X_test, y_train, y_test = pass

4/ Trénování vybraného modelu¶

In [ ]:

# import modelu a metrik evaluace F1/accuracy
pass 

# import modelu a metrik evaluace F1/accuracy pass

In [ ]:

# učení modelu 
pass 

# učení modelu pass

In [ ]:

# testovací predikce 
pass 

# testovací predikce pass

5/ Ladění modelu¶

In [ ]:

# Ladění hyperparametrů 
# Analýza nejlešího modelu 
# Evaluace na testovacím datasetu 

# Ladění hyperparametrů # Analýza nejlešího modelu # Evaluace na testovacím datasetu

In [ ]:

# Grid search CV 
# https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html

pass 

# hyperparametry 
hyperparameter_space = pass 

# grid search 
gs = pass

# Grid search CV # https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html pass # hyperparametry hyperparameter_space = pass # grid search gs = pass

In [ ]:

# Optimální parametry
print("Optimální kombinace: ", gs.best_params_)
print("CV přesnost : ", round(gs.best_score_ * 100, 2))

# Optimální parametry print("Optimální kombinace: ", gs.best_params_) print("CV přesnost : ", round(gs.best_score_ * 100, 2))

In [ ]:

# Test přesnosti 
pass 

# Test přesnosti pass

6/ Interpretace výsledků¶

In [ ]:

# Predikce a interpretace indikátorů důležitosti 
pass 

# Predikce a interpretace indikátorů důležitosti pass

In [ ]:

# Vizualizace 
pass 

# Vizualizace pass

In [ ]:

# Komentáře kvality modelu 
# ...

# Komentáře kvality modelu # ...