Основною метою дослідження було визначити в межах Кривого Рогу загальну площу території, яку покрито “живою” (не засохлою) рослинністю. При цьому особливим завданням було виокремити території, які покриті деревною рослинністю. Для цього упродовж декількох днів я візуально вивчав фактичні дані про стан газонів, скверів, пустирів та інших об’єктів у Кривому Розі. Під час візуального обстеження стало зрозумілим, що лише мала частка власне газонів представлена “живою” травою. Більшість спостережених мною ділянок несли або вельми пригнічену, або взагалі висохлу трав’яну рослинність. Також стало зрозумілим, що в окремий клас території треба винести пустирі та очерети. Пустирі здебільшого вкриті або багаторічними високорослими травами (деревій, полин, амброзія тощо), або сумішшю низькорослих та високорослих трав. Також майже всюди на пустирях трапляються поодинокі дерева. Очерети займають особливе місце у міській рослинності, оскільки вони вказують на близьке знаходження ґрунтових вод (зони підтоплення). Також очерети, зростаючи на берегах та у водоймах, вказують на ступінь замулення поверхневих водойм.

Через малу кількість “живих” газонів довелось трохи змістити рівень детектування класу “газон”. Під час класифікації до “газонів” мною було віднесено не лише власне “газони”, а й будь-які ділянки з низькорослими “живими” травами - доглянуті стадіони, частини пустирів, ділянки суходолу в безпосередній близькості до водних об’єктів.

Загалом в рамках дослідження було виділено 9 умовних класів:

В якості первинних матеріалів взято:

Для радарного знімку було проведено радіометричну корекцію, фільтрацію спекл-шуму (з використанням фільтру Gamma Map), обчислення сили відбиття в децибелах та геометричну корекцію (Range-Doppler). Для цього було використано програмний засіб SNAP (Sentinels Application Platform).

Знімки оптичних супутників оброблено з використанням модулю Semi-Automatic Classification Plugin геоінформаційної системи QGIS (Sentinel-2A) та SAGA GIS (Landsat-8).

Після підготовчого оброблення всі знімки було завантажено в геоінформаційну систему SAGA. Там виконано поканальне об’єднання двох гранул знімків космічного апарату Sentinel-2A (TWT та UWU), обрізання всіх знімків по контуру Кривого Рогу (контур взято з проекту OpenStreetMap), приведення всіх знімків до розміру пікселю 10Х10 м. При цьому в якості “основи” використовувався знімок космічного апарату Sentinel-2A.

Для формування навчальної вибірки об’єктів зроблено сегментацію досліджуваної території базовим методом Object Based Image Segmentation геоінформаційної системи SAGA. Потім незначній частині полігонів отриманого шейп-шару приписано означені вище класи території. Загалом було прокласифіковано 1070 полігонів. Кількість полігонів за класами значно різниться.

Наступні кроки виконувались в середовищі мови програмування R. Нижче містяться скрипти, якими виконувався весь процес обчислення та перевірки даних.

#Завантажуємо потрібні пакети
library(plyr)#саме в такій послідовності: спочатку plyr, потім dplyr!
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:plyr’:

    arrange, count, desc, failwith, id, mutate, rename, summarise, summarize

The following objects are masked from ‘package:raster’:

    intersect, select, union

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union
library(raster)
library(ggplot2)
#Встановлюємо директорію для кешування "великих растрів"
rasterOptions(tmpdir = "/home/geka/ssd/LargeTMP/")
#Завантажуємо у вигляді окремих об'єктів всі супутникові матеріали
Image_S2B2 <- raster("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/S2_TWT_UWU_20170808_B02_KrR.sdat")
Image_S2B3 <- raster("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/S2_TWT_UWU_20170808_B03_KrR.sdat")
Image_S2B4 <- raster("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/S2_TWT_UWU_20170808_B04_KrR.sdat")
Image_S2B5 <- raster("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/S2_TWT_UWU_20170808_B05_KrR.sdat")
Image_S2B6 <- raster("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/S2_TWT_UWU_20170808_B06_KrR.sdat")
Image_S2B7 <- raster("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/S2_TWT_UWU_20170808_B07_KrR.sdat")
Image_S2B8 <- raster("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/S2_TWT_UWU_20170808_B08_KrR.sdat")
Image_S2B8A <- raster("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/S2_TWT_UWU_20170808_B8A_KrR.sdat")
Image_S2B11 <- raster("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/S2_TWT_UWU_20170808_B11_KrR.sdat")
Image_S2B12 <- raster("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/S2_TWT_UWU_20170808_B12_KrR.sdat")
Image_L8B10 <- raster("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/RT_LC08_L1TP_178027_20170812_20170812_01_RT_B10_KrR.sdat")
Image_L8B11 <- raster("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/RT_LC08_L1TP_178027_20170812_20170812_01_RT_B11_KrR.sdat")
Image_S1VH <- raster("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/S1A_GRDH_20170807_SigmaVH_db_KrR.sdat")
Image_S1VV <- raster("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/S1A_GRDH_20170807_SigmaVV_db_KrR.sdat")
Image_S2B8APR <- raster("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/S2_TWT_UWU_20170430_B08_KrR.sdat")
#Об'єднуємо всі канали в один об'єкт
ImageStack <- stack(Image_S2B2, Image_S2B3, Image_S2B4, Image_S2B5, Image_S2B6,
                    Image_S2B7, Image_S2B8, Image_S2B8A, Image_S2B11, Image_S2B12,
                    Image_L8B10, Image_L8B11, Image_S1VH, Image_S1VV, Image_S2B8APR)
#Перейменовуємо канали в об'єднаному об'єкті
names(ImageStack) <- c("S2B2", "S2B3", "S2B4", "S2B5", "S2B6", "S2B7",
                       "S2B8", "S2B8A", "S2B11", "S2B12",
                    "L8B10", "L8B11", "S1VH", "S1VV", "S2B8APR")
#Візуально перевіряємо адекватність та повноту даних
plotRGB(ImageStack, r = 7, g = 11, b = 14, stretch = "hist")

#Завантажуємо шейп-файл з навчальними полігонами
trainShape <- shapefile("/home/geka/Documents/ZhurnRozslid/DataStory/Ecology/Urban/Gazon20170808/Segments_Class_noDegLawn.shp")
#Відкидаємо "порожні" полігони - ті, яким не присвоювалось жодних класів
trainData <- trainShape[!is.na(trainShape@data$Class) ,]
#Робимо поле цілих чисел з номерами класів
trainData@data$ClassInt <- as.integer(factor(trainData@data$Class))
#Перевіряємо кількість наявних учбових полігонів за класами
summary(factor(trainData@data$Class))
    Bulrush Burned Area      Forest        Lawn        Park       Solid     Village 
        124          42         130          87         181         124         126 
  Wasteland       Water 
        129         127 
#Попіксельно витягуємо дані в межах учбових полігонів та формуємо зведену таблицю значень з супутників та відповідних класів. 
#Приклад програмного коду взято з: http://amsantac.co/blog/en/2015/11/28/classification-r.html
dfAll <- data.frame(matrix(vector(), nrow = 0, ncol = length(names(ImageStack)) + 1))
for (i in 1:length(unique(trainData[["ClassInt"]]))){
  category <- unique(trainData[["ClassInt"]])[i]
  categorymap <- trainData[trainData[["ClassInt"]] == category,]
  #запускаємо кластер багатопоточних обчислень
  beginCluster()
  dataSet <- extract(ImageStack, categorymap)
  #відключаємо кластер багатопоточних обчислень
  endCluster()
  dataSet <- dataSet[!unlist(lapply(dataSet, is.null))]
  dataSet <- lapply(dataSet, function(x){cbind(x, ClassInt = as.numeric(rep(category, nrow(x))))})
  df <- do.call("rbind", dataSet)
  dfAll <- rbind(dfAll, df)
}
Loading required namespace: parallel
4 cores detected, using 3
Using cluster with 3 nodes
4 cores detected, using 3
Using cluster with 3 nodes
4 cores detected, using 3
Using cluster with 3 nodes
4 cores detected, using 3
Using cluster with 3 nodes
4 cores detected, using 3
Using cluster with 3 nodes
4 cores detected, using 3
Using cluster with 3 nodes
4 cores detected, using 3
Using cluster with 3 nodes
4 cores detected, using 3
Using cluster with 3 nodes
4 cores detected, using 3
Using cluster with 3 nodes
#Повертаємо символьні найменування класів в отриману кроком вище таблицю
ClassDF <- data.frame(Class = levels(factor(trainData@data$Class)), ClassInt = seq(1:9))
dfAll <- left_join(dfAll, ClassDF)
Joining, by = c("ClassInt", "Class")
LearnDF <- dplyr::select(dfAll, -ClassInt)
#Перевіряємо кількість пікселів кожного класу в отриманій таблиці учбової вибірки
summary(factor(LearnDF$Class))
    Bulrush Burned Area      Forest        Lawn        Park       Solid     Village 
       2398         871        2671        1961        3615        4366        3153 
  Wasteland       Water 
       3049        3444 
#Будуємо гістограми зі статистичними показниками виділених класів
# png("Histogram_Gazon_WithoutDegLawn_20170901.png", width = 297,
#      height = 210, units = "mm", res = 150)
tbl_df(LearnDF) %>%
  reshape2::melt(id = "Class") %>%
  ggplot(aes(x = value, fill = variable)) +
  geom_histogram(position = "identity") +
  facet_wrap( ~ variable, nrow = 3, scales = "free") +
  labs(y = "кількість випадків\n",
       x = "\nвеличини значень оптичних та радарних діапазонів",
       title = "Розподіл значень радарних та оптичних сигналів\nв межах учбових полігонів умовних класів території Кривого Рогу",
       caption = "дані космічних апаратів\nSentinel-1A, Sentinel-2A, Landsat-8\n30 квітня та 7-12 серпня 2017 року") +
  theme(plot.title = element_text(hjust = 0.5),
        text = element_text(family = "FreeSans",
     face = "plain", size = 12, colour = "firebrick4", lineheight = 0.9),
     legend.position = "none",
     axis.text.x = element_text(colour = "gray15", face="bold",
     size= 8, angle = 0, vjust = 0.5, hjust = 0.5),
     axis.text.y = element_text(colour = "gray15", face="bold",
     size= 10, angle = 0, vjust = 0.5, hjust = 1),
     axis.ticks = element_line(colour = "orange", size = 0.1),
     plot.background = element_rect(fill = "gray95"),
     panel.grid.major = element_line(colour = "orange", size = 0.1),
     panel.grid.minor = element_line(colour = "gray90"),
     panel.background=element_rect(fill="gray95"))

# dev.off()
#Будуємо діаграму зі статистичними показниками виділених класів
# png("Boxplot_Gazon_WithoutDegLawn_20170901.png", width = 297,
#      height = 210, units = "mm", res = 150)
tbl_df(LearnDF) %>%
  reshape2::melt(id = "Class") %>%
  ggplot(aes(y = value, x = Class, colour = variable)) +
  geom_boxplot(outlier.size = 0.5) +
  coord_flip() +
  facet_wrap( ~variable, nrow = 3, scales = "free_x") +
  labs(y = "величини значень оптичних та радарних діапазонів\n",
       x = "\n",
       title = "Радарні та оптичні властивості\nумовних класів території в межах Кривого Рогу",
       caption = "дані космічних апаратів\nSentinel-1A, Sentinel-2A, Landsat-8\n30 квітня та 7-12 серпня 2017 року") +
  theme(plot.title = element_text(hjust = 0.5),
        text = element_text(family = "FreeSans",
     face = "plain", size = 12, colour = "firebrick4", lineheight = 0.9),
     legend.position = "none",
     axis.text.x = element_text(colour = "gray15", face="bold",
     size= 8, angle = 0, vjust = 0.5, hjust = 0.5),
     axis.text.y = element_text(colour = "gray15", face="bold",
     size= 10, angle = 0, vjust = 0.5, hjust = 1),
     axis.ticks = element_line(colour = "orange", size = 0.1),
     plot.background = element_rect(fill = "gray95"),
     panel.grid.major = element_line(colour = "orange", size = 0.1),
     panel.grid.minor = element_line(colour = "gray90"),
     panel.background=element_rect(fill="gray95"))