Wstęp

Celem mojego projektu jest obliczenie proporcji terenów zielonych do terenu zabudowanego w Łodzi i okolicach, używając metody clara do klastrownia kolorów na obrazku.

library(cluster)
library(png)
## Warning: pakiet 'png' został zbudowany w wersji R 4.5.2
library(rasterImage)
## Ładowanie wymaganego pakietu: plotrix
library(hues)

Przygotowanie mapy do analizy

Do wykonania obliczeń użyłam zrzutu z ekranu z Google Maps. W pierwszej wersji analizy uzywałam mapy geodezyjnej z zaznaczonymi parkami i lasami, ale przez zacienienie terenow poza Łodzią, dawała zafałszowane wyniki. Użycie zrzutu ekranu z Google Maps ma też taką zaletę, że pokazuje to, w jakich obszarach jest dużo zabudowań, a w jakich nie. W Łodzi duża koncentracja budynków wystepuje tylko w Srudmieściu i bliskich okolicach, a wiele osiedli przy granicy miasta wygląda jak podmiejskie wsie. Dlatego wjeżdzając nawet 5 kilimetrów w głąb miasta (a trasa z jednego końca Łodzi na drugi wynisi około 20 kilometrów) można spotkać głównie osiedla z domkami jednorodzinnymi.

#Zaimportowanie obrazka
mapa <- readPNG("mapa łodzi.png")

Kod poniżej pobiera wymiary mapy (wysokość, szerokość, ilość kolorów razem z przezroczystością). Wysokość i szerokosc pomnożone przez siebie dają nam ilość pikseli (obserwacji), które będą póżniej dopasowywane do klastrów.

dm1 <- dim(mapa) 
print(dm1)
## [1] 497 567   4

Kluczowym krokiem na tym etapie jest zamienienie mapy w data frame - dzięki temu analiza będzie łatwiejsza. Każdej obserwacji jest przyporządkowywany kolor - kombinacja koordynat r (czerwony), g (zielony) i b (niebieski). Na tej podstawie póżniej przyporządkujemy piksele do odpowiednich klastrów.

Następnie nakładamy mapę na układ współrzędnych, zgodnie z pobranymi wcześniej wymiarami.

rgb_mapa <- data.frame(x = rep(1 : dm1[2], each = dm1[1]),
                           y = rep(dm1[1] : 1, times = dm1[2]),
                           r.value = as.vector(mapa[, , 1]),
                           g.value = as.vector(mapa[, , 2]),
                           b.value = as.vector(mapa[, , 3]))
head(rgb_mapa)
##   x   y  r.value   g.value   b.value
## 1 1 497 0.827451 0.9686275 0.8784314
## 2 1 496 0.827451 0.9686275 0.8784314
## 3 1 495 0.827451 0.9686275 0.8784314
## 4 1 494 0.827451 0.9686275 0.8784314
## 5 1 493 0.827451 0.9686275 0.8784314
## 6 1 492 0.827451 0.9686275 0.8784314
dim(rgb_mapa) 
## [1] 281799      5
plot(y ~ x,
     data = rgb_mapa,
     main = "Mapa Łodzi",
     col = rgb(rgb_mapa[c("r.value", "g.value", "b.value")]),
     asp = 1, pch = ".")

Ustalenie optymalnej liczby klastrów

W tym celu użyłam metody silhouette. Kod poniżej sprawdza, jak byłyby pogrupowane kolory dla różnych ilości klastrów (od 2 do 10). Grupowanie wykonuje algorytm clara, który jest często stosowany w przypadku segmentacji obrazów - jest skuteczny przy przetwarzaniu duzych zbiorów danych, jakimi zazwyczaj są zdjęcia oraz medoidami są barwy już wystepujące na analizowanym obrazku.

n1 <- c()                           
for (i in 2:10) {               
  cl <- clara(rgb_mapa[, c("r.value", "g.value", "b.value")], i)    
  n1[i] <- cl$silinfo$avg.width     
}

plot(n1, type = 'l', main = "Optimal number of clusters",
     xlab = "Number of clusters", ylab = "Average silhouette", col = "blue")
points(n1, pch = 21, bg = "navyblue")
abline(h = (1:30)*5/100, lty = 3, col = "grey50")

Na podstawie wykresu widać, że optymalna będzie każda liczba klastrów powyżej 3. Dla naszej analizy nie ma to aż tak dużego znaczenia, bo wszystkie tereny zielone i zabudowane są zaznaczone jednym kolorem (jednym odcieniem zieleni i jednym odcieniem szarosci). Dla przejrzystości wybiorę 3 klastry.

clara <- clara(rgb_mapa[, 3:5], k = 3) 
plot(silhouette(clara))

Wykres wyników silhouette dla 3 klastrów pokazuje, że piksele zostały raczej dobrze przypasowane (najmniejsze wyniki sa w drugim klastrze).

W kodzie poniżej nastepuje właściwe przypasowanie pikseli do klastrów oraz wytypowanie medoidów, czyli najczęściej wystepujących na mapie kolorów.

clara$medoids 
##        r.value   g.value   b.value
## [1,] 0.8274510 0.9686275 0.8784314
## [2,] 0.6666667 0.7372549 0.8313725
## [3,] 0.9647059 0.9647059 0.9647059
rgb(clara$medoids) 
## [1] "#D3F7E0" "#AABCD4" "#F6F6F6"
swatch(rgb(clara$medoids))

r <- clara$medoids[1, 1] * 255
g <- clara$medoids[1, 2] * 255
b <- clara$medoids[1, 3] * 255

Kolor pierwszy obejmuje wszystkie tereny zielone, trzeci tereny zabudowane z dużą koncentracją budynków, a drugi największe drogi i nazwy wystepujące na mapie. W pierszej wersji projektu, kiedy wybrałam liczbe 2 klastrów, drogi i nazwy zostały włączone do klastra pierwszego (zielonego), anie jak liczyłam do klastra terenów zabudowanych (szarego).

Obrazek pomiżej przedstawia mapę po klastrowaniu.Widać na niej, napisy wystepuja raczej na plamach koloru szarego, więc we wnioskach mówiąc “tereny niezabudowane” będę sie odnosić do klastra drugiego i trzeciego jednocześnie.

colors <- rgb(clara$medoids[clara$clustering, ])
plot(rgb_mapa$y ~ rgb_mapa$x, 
     col = colors, pch = ".", cex = 2, asp = 1, main = "Mapa Łodzi")

Przedstawienie proporcji terenów zielonych do zabudowanych

Kod poniżej oblicza skład procentowy, w jakim trzy kolory pokrywają obrazek. Wyniki można zinterpretować jako przybliżona proporcję terenów zielonych do terenów zabudowanych w Łodzi i jej okolicach.

pixel_colours <- rgb(clara$medoids[clara$clustering, ])
dominantColours <- as.data.frame(table(pixel_colours))
colnames(dominantColours) <- c("colours", "Freq")
total_pixels <- sum(dominantColours$Freq)
dominantColours$distribution <- round((dominantColours$Freq / total_pixels) * 100, 2)

dominantColours$colours <- as.character(dominantColours$colours)
pie(dominantColours$Freq, labels = dominantColours$distribution,
    col = dominantColours$colours,
    xlab = "Colours",
    ylab = "Frequency")

Wykres pokazuje, że proporcja terenów zielonych do zabudowanych (z dużą koncentracja budynków oraz dróg) wynosi 50 na 50. Żeby ustalić czy to dobry wynik, czy nie, należałoby wykonać podobną analizę dla innych dużych polskich miast.