Projekt

Moja diplomová práca sa má venovať deskriptívnej priestorovej analýze triedenia obyvateľstva v Českej republike na základe sčítania ľudu a údajov o kvalite ovzdušia. V projekte som sa zatiaľ venovala analýze údajov o kvalite ovzdušia.

Na začiatok som si načítala potrebné balíčky a nahrala svoje dáta, ktoré sú vo formáte shp a predstavujú 5-ročné priemerné koncentrácie v sieti 1x1 km (za obdobie 2018-2022).

# Nastavenie zrkadlového servera CRAN
options(repos = c(CRAN = "https://cran.rstudio.com/"))

# Inštalácia a načítanie potrebných balíkov
install.packages("dplyr")

## package 'dplyr' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\nikiz\AppData\Local\Temp\Rtmpi6VXBt\downloaded_packages

install.packages("patchwork")

## package 'patchwork' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\nikiz\AppData\Local\Temp\Rtmpi6VXBt\downloaded_packages

library(patchwork)
library(sf)
library(ggplot2)
library(viridisLite)
library(dplyr)

shp_data <- st_read("~/Diplomová práca/Dáta/sit1000_5lprum_18_22_CReko_JTSK/sit1000_5lprum_18_22_JTSK_eko.shp", quiet = TRUE)

# Zobrazenie prvých pár riadkov dát
head(shp_data)

## Simple feature collection with 6 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -774875.3 ymin: -1211324 xmax: -771274.5 ymax: -1209227
## Projected CRS: S-JTSK / Krovak East North
##   ID  CISLO SO2_rp_5l SO2_zp_5l NOx_rp_5l                       geometry
## 1  0 450380         2       1.8       4.8 MULTIPOLYGON (((-772491.6 -...
## 2  0 451380         2       1.9       5.1 MULTIPOLYGON (((-772491.6 -...
## 3  0 448381         2       1.8       4.2 MULTIPOLYGON (((-774275.7 -...
## 4  0 449381         2       1.8       4.7 MULTIPOLYGON (((-773324.9 -...
## 5  0 450381         2       1.8       4.8 MULTIPOLYGON (((-772395.5 -...
## 6  0 451381         2       1.9       5.1 MULTIPOLYGON (((-771371.6 -...

Vysvetlivky

Znečišťující látky, které mají stanoven imisní limit pro ochranu ekosystémů a vegetace:

SO2_rp_5l SO2 - roční průměrná koncentrace [μg.m^-3] -> oxid siričitý;
SO2_zp_5l SO2 – průměrná koncentrace za zimní období (1.10.-31.3.);
NOx_rp_5l NOx - roční průměrná koncentrace [μg.m^-3] - “NOx” sa používa na označenie skupiny znečisťujúcich látok z oxidu dusíka (NO) a oxidu dusičitého (NO2), ktoré sú produkované pri spaľovaní fosílnych palív v rôznych procesoch vrátane vykurovania, dopravy a priemyselných činností.

Dáta Na zoznámenie sa s dátami som si najprv spravila deskriptívnu štatistiku.

summary_stats <- shp_data %>%
    summarize(
        avg_SO2_rp = mean(SO2_rp_5l, na.rm = TRUE),
        median_SO2_rp = median(SO2_rp_5l, na.rm = TRUE),
        sd_SO2_rp = sd(SO2_rp_5l, na.rm = TRUE),
        avg_NOx_rp = mean(NOx_rp_5l, na.rm = TRUE),
        median_NOx_rp = median(NOx_rp_5l, na.rm = TRUE),
        sd_NOx_rp = sd(NOx_rp_5l, na.rm = TRUE)
    )

print(summary_stats)

## Simple feature collection with 1 feature and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -904593.6 ymin: -1227298 xmax: -431726 ymax: -935244
## Projected CRS: S-JTSK / Krovak East North
##   avg_SO2_rp median_SO2_rp sd_SO2_rp avg_NOx_rp median_NOx_rp sd_NOx_rp
## 1   3.047057           2.9 0.8643994   9.743267           8.6  4.496051
##                         geometry
## 1 POLYGON ((-880478 -1085187,...

1. Mapa ročných priemerných koncentrácií SO2

1.1 Zobrazenie mapy

ggplot(data = shp_data) +
    geom_sf(aes(fill = SO2_rp_5l)) +
    scale_fill_gradientn(colours = viridisLite::viridis(256)) +
    theme_minimal() +
    labs(title = "Ročné priemerné koncentrácie SO2", fill = "SO2 [µg/m⁻³]")

V dátach sa nachádza ale veľa nízkych hodnôt, ktoré kvôli niektorým extrémnym hodnotám nie je dobre v mape vidieť. Preto som na pre ich lepšiu viditeľnosť použila logaritmickú farebnú škálu.

# Použitie logaritmickej farebnej škály - môže pomôcť zviditeľniť nízke hodnoty:
ggplot(data = shp_data) +
    geom_sf(aes(fill = log1p(SO2_rp_5l))) +
    scale_fill_viridis_c(option = "viridis") +
    theme_minimal() +
    labs(title = "Ročné priemerné koncentrácie SO2 (log škála)", fill = "log(SO2 + 1) [µg/m⁻³]")

1.2 Boxplot pre ročnú priemernú koncentráciu SO2

# Overenie štruktúry a obsahu dát - sú tam extrémne hodnoty? -> áno
str(shp_data$SO2_rp_5l)

##  num [1:80096] 2 2 2 2 2 2 2 2 2 2 ...

summary(shp_data$SO2_rp_5l)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.200   2.500   2.900   3.047   3.300  22.800

ggplot(data = shp_data, aes(y = SO2_rp_5l)) +
    geom_boxplot() +
    labs(title = "Boxplot koncentrácií SO2",
         y = "SO2 [µg/m⁻³]") +
    theme_minimal()

Keďže sa nachádza niekoľko vyšších hodnôt v dátach tak som si dala zobraziť boxplot aj bez nich.

# Identifikácia extrémnych hodnôt
outliers <- boxplot.stats(shp_data$SO2_rp_5l)$out

# Filtrovanie extrémnych hodnôt
shp_data_filtered <- shp_data[!shp_data$SO2_rp_5l %in% outliers, ]

# Vytvorenie boxplotu bez extrémnych hodnôt
ggplot(data = shp_data_filtered, aes(y = SO2_rp_5l)) +
    geom_boxplot() +
    labs(title = "Boxplot koncentrácií SO2 (bez extrémnych hodnôt)",
         y = "SO2 [µg/m⁻³]") +
    theme_minimal()

1.3 Histogram ročnej priemernej koncentrácie SO2

V ďalšom kroku som sa rozhodla spraviť histogram s normovanou hustotou pre SO2, aby som lepšie pochopila rozloženie hodnôt ročných priemerných koncentrácií SO2 v študovanej oblasti. Normovaná hustota umožňuje porovnať tvar rozloženia dát s teoretickými distribúciami, čo môže pomôcť identifikovať odchýlky a anomálie.

ggplot(data = shp_data, aes(x = SO2_rp_5l, y = ..density..)) +
    geom_histogram(bins = 30, fill = "purple", alpha = 0.7) +
    theme_minimal() +
    labs(title = "Histogram ročných priemerných koncentrácií SO2", x = "SO2 [µg/m⁻³]", y = "Density")

# Histogram s normovanou hustotou - detailnejšia os X
ggplot(data = shp_data, aes(x = SO2_rp_5l, y = ..density..)) +
    geom_histogram(bins = 30, fill = "violet", alpha = 0.7) +
    scale_x_continuous(breaks = seq(ceiling(min(shp_data$SO2_rp_5l)), floor(max(shp_data$SO2_rp_5l)), by = 1)) +
    theme_minimal() +
    labs(title = "Histogram ročných priemerných koncentrácií SO2", x = "SO2 [µg/m⁻³]", y = "Density")

# Histogram s normovanou hustotou pre SO2 s explicitnými limitmi osi X
ggplot(data = shp_data, aes(x = SO2_rp_5l, y = ..density..)) +
    geom_histogram(bins = 30, fill = "green", alpha = 0.7) +
    scale_x_continuous(breaks = seq(floor(min(shp_data$SO2_rp_5l)), ceiling(max(shp_data$SO2_rp_5l)), by = 1),
                       limits = c(floor(min(shp_data$SO2_rp_5l)), ceiling(max(shp_data$SO2_rp_5l)))) +
    theme_minimal() +
    labs(title = "Histogram ročných priemerných koncentrácií SO2", x = "SO2 [µg/m⁻³]", y = "Density")

2. Mapa zimných priemerných koncentrácií SO2

K dispozícii je aj mapa zimných priemerných koncentrácií SO2 a to za obdobie (1.10.-31.3.).

2.1 Zobrazenie mapy

# Mapa zimných priemerných koncentrácií SO2
ggplot(data = shp_data) +
    geom_sf(aes(fill = SO2_zp_5l)) +
    scale_fill_gradientn(colours = viridisLite::viridis(256)) +
    theme_minimal() +
    labs(title = "Zimné priemerné koncentrácie SO2", fill = "SO2 [µg/m⁻³] (1.10.-31.3.)]")

Aj v tomto prípade som spravila logaritmizáciu.

# Použitie logaritmickej farebnej škály - môže pomôcť zviditeľniť nízke hodnoty:
ggplot(data = shp_data) +
    geom_sf(aes(fill = log1p(SO2_zp_5l))) +
    scale_fill_viridis_c(option = "viridis") +
    theme_minimal() +
    labs(title = "Zimné priemerné koncentrácie SO2 (log škála)", fill = "log(SO2 + 1) (1.10.-31.3.)]")

3. Mapa ročných priemerných koncentrácií NOx

Rovnakú analýzu som následne použila aj na mapu ročných priemerných koncentrácií NOx.

3.1 Zobrazenie mapy

# Mapa ročných priemerných koncentrácií NOx

ggplot(data = shp_data) +
    geom_sf(aes(fill = NOx_rp_5l)) +
    scale_fill_gradientn(colours = viridisLite::viridis(256)) +
    theme_minimal() +
    labs(title = "Ročné priemerné koncentrácie NOx", fill = "NOx [µg/m⁻³]")

# Použitie logaritmickej farebnej škály - môže pomôcť zviditeľniť nízke hodnoty:
ggplot(data = shp_data) +
    geom_sf(aes(fill = log1p(NOx_rp_5l))) +
    scale_fill_viridis_c(option = "viridis") +
    theme_minimal() +
    labs(title = "Ročné priemerné koncentrácie NOx (log škála)", fill = "log(NOx + 1) [µg/m⁻³]")

3.2 Boxplot pre mapu ročných priemerných koncentrácií NOx

ggplot(data = shp_data) +
    geom_boxplot(aes(y = NOx_rp_5l)) +
    theme_minimal() +
    labs(title = "Boxplot ročných priemerných koncentrácií NOx", y = "NOx [µg/m⁻³]")

# Identifikácia extrémnych hodnôt
outliers_nox <- boxplot.stats(shp_data$NOx_rp_5l)$out

# Filtrovanie extrémnych hodnôt
shp_data_filtered_nox <- shp_data[!shp_data$NOx_rp_5l %in% outliers_nox, ]

# Vytvorenie boxplotu bez extrémnych hodnôt
ggplot(data = shp_data_filtered_nox, aes(y = NOx_rp_5l)) +
    geom_boxplot() +
    labs(title = "Boxplot koncentrácií NOx (bez extrémnych hodnôt)",
         y = "NOx [µg/m⁻³]") +
    theme_minimal()

3.3 Histogram pre mapu ročných priemerných koncentrácií NOx

ggplot(data = shp_data, aes(x = NOx_rp_5l, y = ..density..)) +
    geom_histogram(bins = 30, fill = "red", alpha = 0.7) +
    scale_x_continuous(breaks = seq(ceiling(min(shp_data$NOx_rp_5l)), floor(max(shp_data$NOx_rp_5l)), by = 1)) +
    theme_minimal() +
    labs(title = "Histogram ročných priemerných koncentrácií NOx", x = "NOx [µg/m⁻³]", y = "Density")

4. Porovnanie boxplotov vedľa seba

# Filtrácia dát na odstránenie extrémnych hodnôt SO2
shp_data_filtered <- shp_data %>%
    filter(SO2_rp_5l < quantile(SO2_rp_5l, 0.99, na.rm = TRUE))

# Filtrácia dát na odstránenie extrémnych hodnôt pre zimné koncentrácie SO2
shp_data_filtered_winter_so2 <- shp_data %>%
    filter(SO2_zp_5l < quantile(SO2_zp_5l, 0.99, na.rm = TRUE))

# Filtrácia dát na odstránenie extrémnych hodnôt NOx
shp_data_filtered_nox <- shp_data %>%
    filter(NOx_rp_5l < quantile(NOx_rp_5l, 0.99, na.rm = TRUE))

# Vytvorenie jednotlivých grafov

# Prvý boxplot: koncentrácie SO2 (bez extrémnych hodnôt)
plot1 <- ggplot(data = shp_data_filtered, aes(y = SO2_rp_5l)) +
    geom_boxplot() +
    labs(title = "Koncentrácia SO2\n(bez extrémnych hodnôt)",
         y = "SO2 [µg/m⁻³]") +
    theme_minimal()

# Druhý boxplot: koncentrácie SO2 (zimné obdobie bez extrémnych hodnôt)
plot2 <- ggplot(data = shp_data_filtered_winter_so2, aes(y = SO2_zp_5l)) +
    geom_boxplot() +
    labs(title = "Zimná koncentrácia SO2\n(bez extrémnych hodnôt)",
         y = "SO2 [µg/m⁻³] (1.10.-31.3.)") +
    theme_minimal()

# Tretí boxplot: koncentrácie NOx (bez extrémnych hodnôt)
plot3 <- ggplot(data = shp_data_filtered_nox, aes(y = NOx_rp_5l)) +
    geom_boxplot() +
    labs(title = "Koncentrácia NOx\n(bez extrémnych hodnôt)",
         y = "NOx [µg/m⁻³]") +
    theme_minimal()

# Zobrazenie troch boxplotov vedľa seba
combined_plot <- plot1 + plot2 + plot3 + plot_layout(ncol = 3)

# Zobrazenie kombinovaného grafu
print(combined_plot)

Geografická analýza znečistenia ovzdušia v ČR

Nikolett Zubnárová

2024-06-13