Припрема из претходних свезака¶

Све свеске се међусобно надограђују, тако да нам је неопходно да поновимо неке делове из претходних свески, што без додатног објашњења чинимо у овој секцији.

Овде нам је ради анализе резултата неопходан списак класа. Такође, да би смо могли да користимо модел поновићемо поступак из ранијих свезака - дакле, увеземо неопходне библиотеке и учитамо модел. Поред тога морамо поновити цео поступак учитавања сателитског снима Србије из претходне свеске. Ово се наводи без објашњења, с обзиром да је поступак идентичан као у претходним свескама

In [1]:

classes = [ 'AnnualCrop',
            'Forest',
            'HerbaceousVegetation',
            'Highway',
            'Industrial',
            'Pasture',
            'PermanentCrop',
            'Residential',
            'River',
            'SeaLake']

In [2]:

from os import environ
environ["OPENCV_IO_ENABLE_JASPER"] = "true"
import torch
import torch.nn as nn
import torch.nn.functional as F
import torchvision

import cv2

import numpy as np

from skimage import exposure
from sklearn import metrics

from matplotlib import pyplot as plt

In [3]:

with torch.no_grad():
    model = torchvision.models.resnet50(pretrained=True)
    num_features = model.fc.in_features
    model.fc = nn.Linear(num_features, 10)

device = "cpu"
model.to(device)

model.load_state_dict(torch.load(r"mldata\resnet_50_land_use.pt", map_location=torch.device(device)))
model.eval()
print("Модел учитан")

Модел учитан

In [4]:

img_path_b02 = r"mldata\S2A_MSIL2A_20210630T093041_N0301_R136_T34TCQ_20210630T122807.SAFE\GRANULE\L2A_T34TCQ_A031450_20210630T093134\IMG_DATA\R10m\T34TCQ_20210630T093041_B02_10m.jp2"
img_path_b03 = r"mldata\S2A_MSIL2A_20210630T093041_N0301_R136_T34TCQ_20210630T122807.SAFE\GRANULE\L2A_T34TCQ_A031450_20210630T093134\IMG_DATA\R10m\T34TCQ_20210630T093041_B03_10m.jp2"
img_path_b04 = r"mldata\S2A_MSIL2A_20210630T093041_N0301_R136_T34TCQ_20210630T122807.SAFE\GRANULE\L2A_T34TCQ_A031450_20210630T093134\IMG_DATA\R10m\T34TCQ_20210630T093041_B04_10m.jp2"

def load_sentinel_2_image_rgb(img_path_b02, img_path_b03, img_path_b04):
    img_b02 = cv2.imread(img_path_b02, cv2.IMREAD_ANYCOLOR | cv2.IMREAD_ANYDEPTH)
    img_b03 = cv2.imread(img_path_b03, cv2.IMREAD_ANYCOLOR | cv2.IMREAD_ANYDEPTH)
    img_b04 = cv2.imread(img_path_b04, cv2.IMREAD_ANYCOLOR | cv2.IMREAD_ANYDEPTH)

    rgb = np.zeros((img_b02.shape[0], img_b02.shape[1], 3))
    rgb[:, :, 0] = img_b04
    rgb[:, :, 1] = img_b03
    rgb[:, :, 2] = img_b02

    rgb = rgb.astype(np.uint16)

    return rgb

image = load_sentinel_2_image_rgb(img_path_b02, img_path_b03, img_path_b04)

def prepare_rgb(image):
    p_low, p_high = np.percentile(image, (1, 99))
    image = exposure.rescale_intensity(image, in_range=(p_low, p_high))
    image = image / image.max()
    
    image = exposure.adjust_gamma(image, 0.9)
    image += 0.05
    image = np.clip(image, 0.0, 1.0)

    return image

# Шабац - севернији део и река Сава
CROP_START_COL = 9584
CROP_START_ROW = 4060
CROP_SIZE = 512

# Шабац - јужнији део
# CROP_START_COL = 9423
# CROP_START_ROW = 4329
# CROP_SIZE = 512

# Рума - јужни део
# CROP_START_COL = 10383
# CROP_START_ROW = 1626
# CROP_SIZE = 512

# Јарак - цело место и меандар Саве
# CROP_START_COL = 9929
# CROP_START_ROW = 2348
# CROP_SIZE = 512

# Вештачко језеро Ровни
# CROP_START_COL = 9721
# CROP_START_ROW = 9908
# CROP_SIZE = 512

sample_image_1 = image[CROP_START_ROW:CROP_START_ROW+CROP_SIZE, CROP_START_COL:CROP_START_COL+CROP_SIZE, :]

Извршавање модела над подацима за Србију¶

Припрема података за потребе модела¶

Као што знамо, наш модел машинског учења као улаз узима слике величине 64x64, док је наш исечак величине 512x512. Да бисмо обрадили цео исечак морамо га поделити у мање исечке величине 64x64. Таквих исечака ће бити 8 по дужини и 8 по ширини, што значи укупно 64. Можда већ и назирете да ће тих 64 мањих исечака чинити хрпу коју ћемо као тавку дати мрежи и од ње добити одговарајући број предикција.

Функције испод се баве управо овом припремом. Конкретно, функција get_images дати исечак дели у скуп малих плочица величине 64x64, а функција prepare_images_for_dnn пакује ту хрпу плочица у структуру падатака адекватне величине коју вештачка неурална мрежа може да прихвати.

In [5]:

def get_images(image):
    image = prepare_rgb(image)

    rows = image.shape[0] // 64
    cols = image.shape[1] // 64
    all_tiles = np.zeros((rows * cols, 64, 64, 3), dtype=np.float32)

    for r in range(rows):
        for c in range(cols):
            all_tiles[r + rows * c, ...] = image[r*64:r*64+64, c*64:c*64+64, :]
    
    return all_tiles

def prepare_images_for_dnn(images):
    images_prepared = torch.zeros((images.shape[0], 3, 224, 224))

    for i in range(images.shape[0]):

            image_resized = cv2.resize(images[i, ...], (224, 224))
            img_normalized = image_resized.transpose(2, 0, 1).astype(np.float32)

            images_prepared[i, ...] = torch.from_numpy(img_normalized)

    return images_prepared

In [6]:

all_tiles = get_images(image[CROP_START_ROW:CROP_START_ROW+CROP_SIZE, CROP_START_COL:CROP_START_COL+CROP_SIZE, :])
tiles_prepared_for_dnn = prepare_images_for_dnn(all_tiles).to(device)

Извршавање модела¶

Модел над овако припремљеним подацима из Србије извршавамо на исти начин као и у ранијим секцијама.

In [7]:

with torch.no_grad():
    pred = model(tiles_prepared_for_dnn)

In [8]:

_, argmax = torch.max(pred, 1)

In [9]:

argmax

Out[9]:

tensor([4, 7, 7, 7, 7, 7, 7, 4, 9, 8, 6, 7, 7, 7, 7, 7, 3, 9, 8, 8, 7, 7, 7, 7,
        7, 0, 6, 8, 3, 4, 7, 7, 6, 0, 6, 3, 8, 4, 4, 6, 6, 6, 6, 6, 0, 8, 4, 4,
        0, 0, 6, 0, 0, 8, 8, 4, 1, 1, 0, 0, 0, 8, 8, 4])

In [10]:

with torch.no_grad():
    softmaxes = F.softmax(pred, dim=1).cpu()

Анализа сигурности модела¶

Анализу сигурности модела такође спроводимо као и у ранијим поглављима, али ће овде бити јако интересантно упоредити са резултатима од раније.

In [11]:

all_scores = np.take_along_axis(softmaxes.numpy().T, np.expand_dims(argmax.numpy(), 0), 0)
print(all_scores)

[[0.28714052 0.99255097 0.99988437 0.99999654 0.99999726 0.84850115
  0.9999993  0.5522834  0.7030461  0.5255811  0.89717245 0.9999397
  0.9995536  0.9999969  0.9999994  0.98643523 0.45484185 0.9644964
  0.9998373  0.59635746 0.7072009  0.7328156  0.99999774 0.98855084
  0.5047223  0.97150755 0.24233337 0.7680217  0.6728504  0.9999906
  0.95893246 0.99998343 0.91004765 0.9227216  0.52570987 0.84776026
  0.82162917 0.56186706 0.999752   0.99741226 0.9826607  0.93476874
  0.6480117  0.8610074  0.85560554 0.9984549  0.99504954 0.94746214
  0.6534794  0.8140811  0.7918098  0.7742024  0.9847025  0.9995018
  0.99851817 0.85900277 0.89284736 0.9786476  0.7984398  0.93689704
  0.98762786 0.9812507  0.9982487  0.98680216]]

In [12]:

all_scores.mean()

Out[12]:

0.8531328

Самостално анализирајте просечну сигурност модела и његове предикције и упоредите је са раније добијеном сигурношћу када се модел извршавао на валидационом скупу. Дакле, у питању је идентичан модел, примењен на различитим скупом података.

Даје анализу и интерпретацију зашто постоје ове разлике? Узмите у обзир шта би се десило да нисмо применили никакву доменску адаптацију (чак можете и пробати да је искључите и видите шта би се десило).

Визуелизација резултата¶

С обзиром да за податке из Србије немамо лабеле, није могуће урадити квантитативну анализу тачности модела. Зато ћемо се задржати на квалитативној анализи посматрајући резултате модела и видети да ли они одговарају реалности. За овакву анализу је потребно дати што бољу визуелизацију резултата, што ћемо и учинити у овој секцији.

Сваку од класа које модела предвиђа ћемо представити једном бојом да би се јасније видела на растерској подлози. Испод су за сваку од постојећих класа придружене боје у виду триплета који означавају црвену, плаву и зелену боју. Намерно смо ставили да поједине јако сличне класе имају исту боју - но, осећајте се потпуно слободно да вредности испод измените по жељи. Тренутно су зелене боје за разне врсте зелених поврина и ораница, црвена боја за насеља, магента за индустрију, наранџаста за аутопут, а плаве за реке/језера.

In [13]:

classes_to_rgb = { 
                'Highway': (1.0, 0.5, 0), 
                'Industrial': (1.0, 0, 0.5),
                'Residential': (1.0, 0, 0), 

                'AnnualCrop': (0, 1.0, 0.5),
                'PermanentCrop': (0, 1.0, 0.5),
                'Forest': (0, 1.0, 0),
                'HerbaceousVegetation': (0, 1.0, 0),
                'Pasture': (0, 1.0, 0),
                
                'River': (0.5, 0, 1.0),
                'SeaLake': (0.5, 0, 1.0)}

Функције испод служе за лепу визуелизацију резултата. Уколико добро познајете Пајтон, требало би углавном да можете самостално да их разумете, мада то и није неопходно.

In [14]:

def show_results_tiles(all_tiles, predictions_argmax, classes):
    for idx, prediction in enumerate(predictions_argmax):
        plt.figure(figsize=(2, 2))
        plt.imshow(all_tiles[idx])
        plt.show()
        print(classes[prediction])

def show_results_tiles_semantic(all_tiles, predictions_argmax, classes, softmaxes, num_rows):
    # I am packing them row majro
    img = np.zeros((num_rows * 64, (all_tiles.shape[0] // num_rows) * 64, 3))
    img_sem = np.zeros((num_rows * 64, (all_tiles.shape[0] // num_rows) * 64, 3))
    img_conf = np.zeros((num_rows * 64, (all_tiles.shape[0] // num_rows) * 64, 3))
    for idx, prediction in enumerate(predictions_argmax):
        img[(idx % num_rows) * 64 : (idx % num_rows + 1) * 64 , (idx // num_rows) * 64 : (idx // num_rows + 1) * 64, :] = all_tiles[idx] #if classes[prediction] in ["Residential", "Industrial"] else [0, 0, 0]
        img_sem[(idx % num_rows) * 64 : (idx % num_rows + 1) * 64 , (idx // num_rows) * 64 : (idx // num_rows + 1) * 64, :] = classes_to_rgb[classes[prediction]] #if classes[prediction] in ["Residential", "Industrial"] else [0, 0, 0]
        img_conf[(idx % num_rows) * 64 : (idx % num_rows + 1) * 64 , (idx // num_rows) * 64 : (idx // num_rows + 1) * 64, :] = softmaxes[idx, argmax[idx]]
    
    plt.figure(figsize=(10, 10))
    plt.imshow(img)
    plt.show()
    plt.figure(figsize=(10, 10))
    plt.imshow(img * 0.7 + img_sem * 0.3)
    plt.show()
    plt.figure(figsize=(10, 10))
    plt.imshow(img_sem)
    plt.show()
    plt.figure(figsize=(10, 10))
    plt.imshow(img_conf)

Хајде сада да коначно погледамо визуелизације резултата. Добићемо четири слике:

Исечак који обрађујемо
Преклопљене предвиђене лабеле класа у дефинисаним бојама над улазном сликом
Чисте предвиђене лабеле без преклопа
Матрица где светлина ћелије означава сигурност модела. Тамније ћелије означавају мању сигурност.

Погледајте резултате детаљно. Уочите да ли модел добро разликује појединачне површине. Наравно, лакше је разликовати насеље од оранице, али пак обратите пажњу да ли модел успешно разликује индустрију од насеља, што је тежи задатак. Уочите и да ли модел предвиђа остале класе коректно (рецимо саобраћајнице). Детаљно анализирајте матрицу сигурности модела и размислите зашто је на појединим ћелијама та сигурност смањена.

In [15]:

show_results_tiles_semantic(all_tiles, argmax, classes, softmaxes, num_rows=CROP_SIZE // 64)

Гушћи скуп предикција¶

У претходној визуелизације смо модел применили на кораку од 64 пиксела, као кад радимо поплочавање - ни једна од улазних плочица модела се не преклапа. Наравно, могуће је постићи и гушћу расподелу предикција тако што ћемо узети мањи корак модела. Управо то је циљ ове секције.

Најпре ћемо модификовати функцију get_images да уједно и прима корак са који се праве предикције и тиме добијамо функцију get_images_advanced.

In [16]:

def get_images_advanced(image, crop_size=64, stride=32):
    image = prepare_rgb(image)

    rows = image.shape[0] // stride
    cols = image.shape[1] // stride
    all_tiles = np.zeros((rows * cols, crop_size, crop_size, 3), dtype=np.float32)

    for r in range(rows):
        for c in range(cols):
            if r*stride + crop_size > image.shape[0] or c*stride + crop_size > image.shape[1]:
                continue
            all_tiles[r + rows * c, ...] = image[r*stride:r*stride+crop_size, c*stride:c*stride+crop_size, :]
    
    return all_tiles

Функција prepare_images_for_dnn остаје иста и сад те две функције можемо позвати али са кораком 32.

In [17]:

all_tiles_2 = get_images_advanced(image[CROP_START_ROW:CROP_START_ROW+CROP_SIZE, CROP_START_COL:CROP_START_COL+CROP_SIZE, :])
tiles_prepared_for_dnn = prepare_images_for_dnn(all_tiles_2).to(device)

Слично као и у ранијим ћелијама, хајде да извршимо модел. С обзиром да сада треба дати много већи број предикција, очекујте да ће извршавање трајате дуже.

In [18]:

pred = model(tiles_prepared_for_dnn)

In [19]:

_, argmax = torch.max(pred, 1)

Функција испод нам служи за визуелизацију гушћих предикција над улазном мапом. Хајде да погледамо и анализирамо те резултате.

In [20]:

def show_results_tiles_semantic_2(all_tiles_2, predictions_argmax, classes, softmaxes, num_rows, stride):
    plt.figure(figsize=(10, 10))
    plt.imshow(prepare_rgb(sample_image_1) + 0.2)  # Додајемо мало светлине (0.2) да би се тачке лепше виделе на растерској подлози
    for idx, prediction in enumerate(predictions_argmax):
        y = (idx % num_rows) * stride + 32
        x = (idx // num_rows) * stride + 32
        if x >= 512 or y >= 512:
            continue
        plt.scatter(x, y, color=classes_to_rgb[classes[prediction]], s=100)
    plt.show()

In [21]:

show_results_tiles_semantic_2(all_tiles_2, argmax, classes, None, num_rows=CROP_SIZE // 32, stride=32)

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

Уочите да ли се са гушћим предикцијама добијају неке разлике. Најпре, свакако се може приметити да са гушћим предикцијама добијамо много прецизнију детекцију намене површине земље - рецимо сад је река Сава јасније детектована тако да предикције прецизније покривају њено корито.

Prijavi problem

Obeleži sve kategorije koje odgovaraju problemu

Još detalja - opišite nam problem

Data Officer