library(sf) #analizy przestrzenna danych wektorowych w R
## Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3
library(ggplot2)
library(ggResidpanel) #interaktywne wykresy diagnostyczne
library(ggfortify) # wykresy diagnostyczne w stylizacji ggplot2.
Szczegółowe dane dotyczące rozmieszczenia ludności mają wiele zastosowań:
Dane o liczbie ludności zbierane są w trakcie Powszechnych Spisów Ludności od każdego mieszkańca, a potem agregowane do jednostek przestrzennych różnej wielkości; w Polsce są to np. województwa, powiaty, gminy, rejony i obwody spisowe.
Wielkość jednostek spisowych zależy od stopnia agregacji oraz lokalizacji (jednostki spisowe są mniejsze w miejskich obszarach zurbanizowanych oraz większe w obszarach poza miejskich)
Dla wielu praktycznych zastosowań nawet najmniejsze jednostki agregacji są zbyt duże, aby dostarczyć szczegółowych danych o rozmieszczeniu ludności.
Granice jednostek spisowych i administracyjnych (gminy, obszary spisowe) oraz jednostek przyrodniczych (zlewnie, regiony fizyczno-geograficzne) nie pokrywają się.
Metody służące do przekształcania danych zagregowanych do danych o wyższej rozdzielczości można podzielić na dwie główne grupy:
metody powierzchniowo-wagowe (ang. areal interpolation, areal weightening interpolation)
modelowanie dazymetryczne (ang. dasymetric modeling)
Celem ćwiczenia jest przygotowanie szczegółowej mapy rozmieszczenia ludności dla obszaru zlewni górnej Parsęty.
Dane wejściowe:
Znalezienie danych pomocniczych, które można wykorzystać do opracowania szczegółowej mapy rozmieszczenia ludności.
Przygotowanie szczegółowych danych o liczbie ludności.
Analiza zostanie przeprowadzona dla zlewni górnej Parsęty.
Kolorem szarym zaznaczone zostały komórki niezamieszkałe.
Corine Land Cover jest zbyt ogólną mapą; nie może zostać wykorzystana jako dane pomocnicze do przygotowania szczegółowej mapy rozmieszczenia ludności w zlewni górnej Parsęty.
<- st_read("dane_dasy/zlewnia_gp.shp") zlewnia
## Reading layer `zlewnia_gp' from data source `/home/anna/DYDAKTYKA/MODELOWANIE/dane_dasy/zlewnia_gp.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 328886.4 ymin: 652023.5 xmax: 342809.8 ymax: 661330.8
## proj4string: +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +units=m +no_defs
#wczytanie danych z pliku shp
<- st_read("dane_dasy/pop1km_gp3035.shp") pop1km_3035
## Reading layer `pop1km_gp3035' from data source `/home/anna/DYDAKTYKA/MODELOWANIE/dane_dasy/pop1km_gp3035.shp' using driver `ESRI Shapefile'
## Simple feature collection with 105 features and 16 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 4743000 ymin: 3418000 xmax: 4758000 ymax: 3429000
## proj4string: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
#transformacja układów współrzędnych z układu 3035 -> 2180
<- st_transform(pop1km_3035, 2180)
pop1km <- pop1km[,c("TOT", "CODE")]
pop1km head(pop1km)
## Simple feature collection with 6 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 328571.6 ymin: 652422.9 xmax: 330933.3 ymax: 656512.9
## CRS: EPSG:2180
## TOT CODE geometry
## 1 0 CRS3035RES1000mN3421000E4743000 POLYGON ((329689.6 654530.2...
## 2 4 CRS3035RES1000mN3420000E4743000 POLYGON ((329563.9 653538.8...
## 3 14 CRS3035RES1000mN3421000E4744000 POLYGON ((330681.8 654405.6...
## 4 16 CRS3035RES1000mN3422000E4744000 POLYGON ((330807.6 655396.9...
## 5 29 CRS3035RES1000mN3420000E4744000 POLYGON ((330556.1 653414.3...
## 6 61 CRS3035RES1000mN3419000E4744000 POLYGON ((330430.4 652422.9...
<- st_read("dane_dasy/punkty_adresowe_gp.shp") punkty_adresy
## Reading layer `punkty_adresowe_gp' from data source `/home/anna/DYDAKTYKA/MODELOWANIE/dane_dasy/punkty_adresowe_gp.shp' using driver `ESRI Shapefile'
## Simple feature collection with 549 features and 15 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 329409.4 ymin: 651664.3 xmax: 343039.7 ymax: 661702
## proj4string: +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +units=m +no_defs
<- punkty_adresy[,c("jednostkaA", "miejscowos")] punkty_adresy
Etapy przygotowania siatki:
Wybrać w QGIS Vector - Research Tools - Create Grid.
Transformacja siatki do PUWG1992 (EPSG:2180).
<- st_read("dane_dasy/grid100m_2180.shp") grid100m
## Reading layer `grid100m_2180' from data source `/home/anna/DYDAKTYKA/MODELOWANIE/dane_dasy/grid100m_2180.shp' using driver `ESRI Shapefile'
## Simple feature collection with 10500 features and 5 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 328571.6 ymin: 651057.9 xmax: 343832.3 ymax: 661962.3
## proj4string: +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +units=m +no_defs
<- grid100m[,1] grid100m
Siatka 100m z nałożoną siatką 1km (czerwona granica) - Układ współrzędnych - EPSG2180.
Etapy analizy obejmują:
W analizie wykorzystane zostaną dwa zestawy danych:
liczba ludności w siatce 1km na 1km;
punkty adresowe
st_join()
(spatial join) z pakietu sf
to odpowiednik narzędzia Select by Location w QGIS.<- st_join(punkty_adresy, pop1km)
adresy_w_pop1km head(adresy_w_pop1km)
## Simple feature collection with 6 features and 4 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 329994.5 ymin: 651711.6 xmax: 334613.9 ymax: 657409.2
## proj4string: +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +units=m +no_defs
## jednostkaA miejscowos TOT
## 1 (4:Polska,zachodniopomorskie,szczecinecki,Barwice) Ostropole 29
## 2 (4:Polska,zachodniopomorskie,szczecinecki,Barwice) Jeziorki 55
## 3 (4:Polska,zachodniopomorskie,szczecinecki,Barwice) Chwalimki 224
## 4 (4:Polska,zachodniopomorskie,szczecinecki,Barwice) Jeziorki 55
## 5 (4:Polska,zachodniopomorskie,szczecinecki,Barwice) Ostropole 29
## 6 (4:Polska,zachodniopomorskie,szczecinecki,Barwice) Nowy Chwalim 51
## CODE geometry
## 1 CRS3035RES1000mN3420000E4744000 POINT (330133.2 653566.2)
## 2 CRS3035RES1000mN3418000E4747000 POINT (333348.7 651735.6)
## 3 CRS3035RES1000mN3421000E4746000 POINT (332394.2 655051.8)
## 4 CRS3035RES1000mN3418000E4747000 POINT (333114.9 651711.6)
## 5 CRS3035RES1000mN3420000E4744000 POINT (329994.5 653985.1)
## 6 CRS3035RES1000mN3424000E4748000 POINT (334613.9 657409.2)
Zliczając ile razy wystąpił dany kod oczka siatki (kolumna CODE) uzyskamy informację ile jest punktów w danym poligonie.
Powyższe operacje są odpowiednikiem narzędzia Vector - Analysis Tools - Count points in polygon w QGIS.
<- as.data.frame(table(adresy_w_pop1km$CODE))
count_punkty colnames(count_punkty) <- c("CODE", "N_ADRESY")
Połączenie warstwy pop1km zawierającej informacje z liczbą ludności (TOT) w siatce 1x1km oraz liczby punktów adresowych w każdym oczku siatki (N_ADRESY).
<- merge(pop1km, count_punkty, by = "CODE", all.x = TRUE) pop1km
head(pop1km)
## Simple feature collection with 6 features and 3 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 329438.1 ymin: 651057.9 xmax: 333532.8 ymax: 653538.8
## CRS: EPSG:2180
## CODE TOT N_ADRESY geometry
## 1 CRS3035RES1000mN3418000E4746000 0 NA POLYGON ((332289.2 651182.4...
## 2 CRS3035RES1000mN3418000E4747000 55 16 POLYGON ((333281.4 651057.9...
## 3 CRS3035RES1000mN3419000E4744000 61 14 POLYGON ((330430.4 652422.9...
## 4 CRS3035RES1000mN3419000E4745000 38 13 POLYGON ((331422.6 652298.3...
## 5 CRS3035RES1000mN3419000E4746000 16 8 POLYGON ((332414.9 652173.8...
## 6 CRS3035RES1000mN3419000E4747000 0 NA POLYGON ((333407.1 652049.2...
Zastąpienie wartości NA przez 0 dla kolumny N_ADRESY. NA oznacza, że w danym oczku siatki nie było żadnych punktów adresowych (Zauważ, że obiekt count_punkty ma 66 obserwacji (tj. w 66 oczkach siatki wystąpił przynajmniej 1 punkt adresowy), a obiekt pop1km 105 obserwacji).
$N_ADRESY[is.na(pop1km$N_ADRESY)]<- 0
pop1kmsummary(pop1km)
## CODE TOT N_ADRESY geometry
## Length:105 Min. : 0.0 Min. : 0.000 POLYGON :105
## Class :character 1st Qu.: 0.0 1st Qu.: 0.000 epsg:2180 : 0
## Mode :character Median : 7.0 Median : 2.000 +proj=tmer...: 0
## Mean : 22.7 Mean : 5.229
## 3rd Qu.: 20.0 3rd Qu.: 5.000
## Max. :337.0 Max. :54.000
# selekcja danych
<- pop1km[pop1km$TOT>0 & pop1km$N_ADRESY> 0,]
dat
# dodanie numerów wierszy od 1 do ...
rownames(dat) <- 1:nrow(dat)
# Usunięcie kolumny z geometrią. Dane bedą traktowane jak obiekt nieprzestrzenny - ramka danych.
# Funkcja merge nie zadziała na dwóch obiektach klasy przestrzennej sf.
# Dlatego jeden obiekt musi zostać przekształcony do ramki danych.
st_geometry(dat) <- NULL
Chcemy wyjaśnić, jak zmienna TOT (liczba ludności) zależy od zmiennej N_ADRESY (liczba punktów adresowych).
ggplot(dat, aes(x=N_ADRESY, y=TOT)) +
geom_point() +
labs(x = "Liczba punktow adresowych (N_ADRESY)", y = "Liczba ludnosci (TOT)") +
geom_smooth(method = "lm") +
theme_bw()
## `geom_smooth()` using formula 'y ~ x'
<- lm(TOT~N_ADRESY, dat)
model_1 summary(model_1)
##
## Call:
## lm(formula = TOT ~ N_ADRESY, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -61.803 -11.556 -2.030 1.492 230.742
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.4782 6.3106 0.234 0.816
## N_ADRESY 4.1912 0.4763 8.799 1.64e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 38.63 on 62 degrees of freedom
## Multiple R-squared: 0.5553, Adjusted R-squared: 0.5482
## F-statistic: 77.43 on 1 and 62 DF, p-value: 1.644e-12
Pytanie: W jakim stopniu liczba ludności (TOT) jest wyjaśniana przez liczbę punktów adresowych? Jakie są statystyki reszt?
library(ggfortify)
autoplot(model_1) + theme_bw()
library(ggResidpanel)
resid_interact(model_1, plots='R')
library(ggResidpanel)
resid_interact(model_1, plots='all')
## Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomFunction() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
Pytanie: Czego dowiadujemy się z wykresów diagnostycznych? Czy model jest dobrze dopasowany? Czy spełnione są założenia regresji liniowej?
library(car)
## Loading required package: carData
outlierTest(model_1)
## rstudent unadjusted p-value Bonferroni p
## 10 9.774152 4.2405e-14 2.7139e-12
## 12 4.722724 1.4082e-05 9.0127e-04
ggplot(dat, aes(x=N_ADRESY, y=TOT))+
geom_point() +
geom_point(data = dat[c(10, 12, 18),], aes(x=N_ADRESY, y=TOT), color = "red", size = 4) +
labs(x = "Liczba punktow adresowych (N_ADRESY)", y = "Liczba ludnosci (TOT)") +
geom_smooth(method = "lm") +
theme_bw()
## `geom_smooth()` using formula 'y ~ x'
<- dat[!rownames(dat) %in%c(c(10,12, 18)),] dat0
ggplot(dat0, aes(x=N_ADRESY, y=TOT)) +
geom_point() +
labs(x = "Liczba punktow adresowych (N_ADRESY)", y = "Liczba ludnosci (TOT)") +
geom_smooth(method = "lm") +
theme_bw()
## `geom_smooth()` using formula 'y ~ x'
<- lm(TOT~N_ADRESY, dat0)
model_2 summary(model_2)
##
## Call:
## lm(formula = TOT ~ N_ADRESY, data = dat0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -31.971 -6.010 0.144 3.721 37.414
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7019 2.0067 0.35 0.728
## N_ADRESY 3.5769 0.1817 19.68 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.64 on 59 degrees of freedom
## Multiple R-squared: 0.8678, Adjusted R-squared: 0.8656
## F-statistic: 387.4 on 1 and 59 DF, p-value: < 2.2e-16
Model liniowej zależności między liczbą ludności a liczbą punktów adresowych (zbudowany po usunięciu wartości odstających) wzkazuje, że liczba ludności w 86% wyjaśniana jest przez liczbę punktów adresowych (wartość R2). Mediana błędów (reszt) wynosi 0.144 przy rozrzucie błędów między -31 a + 37 osób.
autoplot(model_2) + theme_bw()
resid_interact(model_2, plots='all')
## Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomFunction() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
Do obiektu dat0 dodane zostaną 2 kolumny:
$EST_POP <- model_2$fitted.values
dat0$RES <- model_2$residuals dat0
head(dat0)
## CODE TOT N_ADRESY EST_POP RES
## 1 CRS3035RES1000mN3418000E4747000 55 16 57.932381 -2.932381
## 2 CRS3035RES1000mN3419000E4744000 61 14 50.778571 10.221429
## 3 CRS3035RES1000mN3419000E4745000 38 13 47.201665 -9.201665
## 4 CRS3035RES1000mN3419000E4746000 16 8 29.317139 -13.317139
## 5 CRS3035RES1000mN3420000E4743000 4 2 7.855708 -3.855708
## 6 CRS3035RES1000mN3420000E4744000 29 13 47.201665 -18.201665
#
<- merge(pop1km, dat0[, c("CODE", "EST_POP", "RES")], by = "CODE", all.x = TRUE) out
ggplot(out) +
geom_sf(aes(fill = TOT)) +
scale_fill_distiller(palette = "OrRd", breaks = seq(0, 350, 50), limits=c(0,350), direction=1, na.value = "grey") +
theme_bw()
ggplot(out) +
geom_sf(aes(fill = EST_POP)) +
scale_fill_distiller(palette = "OrRd", breaks = seq(0, 350, 50), limits=c(0,350), direction=1, na.value = "grey") +
theme_bw()
ggplot(out) +
geom_sf(aes(fill = RES)) +
scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, na.value = "grey") +
theme_bw()
Proszę zauważyć na powyższych mapach, że dla kilku oczek siatki wartość estymowana jest NA, podczas gdy istnieje wartość liczby ludności TOT. Model został zbudowany po usunięciu 3 wartości odstających - dla nich zatem nie została wykonana predykcja. Poniżej użyjemy modelu do oszacowania liczby ludności dla 3 punktów usuniętych przy budowaniu modelu.
#ramka danych z danymi do predykcji
<- dat[rownames(dat) %in%c(c(10,12, 18)),]
dat_outlier
#predykcja
<- predict(model_2, dat_outlier) est_pop
#dodanie wartosci predykcji do obiektu wynikowego
$EST_POP[out$CODE%in%dat_outlier$CODE]<- est_pop out
#przypisanie wartosci estymowanej 0 dla punktow dla ktorych TOT=0 oraz N_ADRESY=0 (punkty te nie byly brane do budowy modelu)
$EST_POP[is.na(out$EST_POP)]<- 0 out
#obliczenie reszt
$RES <- out$TOT - out$EST_POP out
st_write(out, "out/out_model_punkty_adresowe.shp", delete_dsn = TRUE)
## Deleting source `out/out_model_punkty_adresowe.shp' using driver `ESRI Shapefile'
## Writing layer `out_model_punkty_adresowe' to data source `out/out_model_punkty_adresowe.shp' using driver `ESRI Shapefile'
## Writing 105 features with 5 fields and geometry type Polygon.
c(suma_TOT = sum(out$TOT, na.rm = TRUE),
suma_EST_POP = sum(out$EST_POP, na.rm = TRUE) )
## suma_TOT suma_EST_POP
## 2383.000 1997.912
Na podstawie danych w 1km siatce w analizowany obszar zamieszkany było przez 2383 osoby, podczas gdy liczba osób oszacowana na podstawie modelu regresji wynosi 1998 osób.
Analiza regresji pozwoliła na ustalenie, że istnieje zależność między liczba ludności, a liczbą punktów adresowych
Informacja ta zostanie wykorzystana do dekompozycji liczby ludności z siatki 1km do siatki 100m.
Wykorzystana zostanie w tym celu procedura nazywana modelowaniem dazymetrycznym.
Funkcja st_centroids() z pakietu sf wyznacza centroidy (Odpowiednik narzędzia Vector - Geometry Tools - Centroids w QGIS).
<- st_centroid(grid100m) grid100m_centroidy
## Warning in st_centroid.sf(grid100m): st_centroid assumes attributes are constant
## over geometries of x
Dla każdego oczka siatki 100m zostanie wyznaczony centroid, a następnie używając narzędzia selekcji przez lokalizację (funkcja st_join()
) do każdego centroidu zostanie przypisany identyfikator oczka 1km.
= st_join(grid100m_centroidy[,1], pop1km)
x $geometry<- NULL x
Informacja ta zostanie przypisana do danych poligonowych grid100m. - Kolumna CODE zawiera kod oczka 1km. - Kolumna id to identyfikator oczka 100m. - Kolumna TOT to liczba ludności w oczku 1km - Kolumna N_ADRESY to liczba adresów w oczku 1km.
<- merge(grid100m, x, by = "id") grid100m
Liczba punktów adresowych w oczku 100m jest obliczana w taki sam sposób jak liczba punktów adresowych w siatce 1km.
<- st_join(punkty_adresy, grid100m)
adresy_grid100m <- as.data.frame(table(adresy_grid100m$id))
count_punkty100 colnames(count_punkty100) <- c("id", "N_ADRESY100")
<- merge(grid100m, count_punkty100, by = "id", all.x = TRUE)
grid100m
$N_ADRESY100[is.na(grid100m$N_ADRESY100)]<-0 grid100m
Kolumna N_ADRESY100 zawiera liczbę punktów adresowych w oczku siatki 100m.
head(grid100m)
## Simple feature collection with 6 features and 5 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 328747.6 ymin: 655038.9 xmax: 328922.3 ymax: 655646.1
## proj4string: +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +units=m +no_defs
## id CODE TOT N_ADRESY N_ADRESY100
## 1 71 CRS3035RES1000mN3421000E4743000 0 1 0
## 2 72 CRS3035RES1000mN3421000E4743000 0 1 0
## 3 73 CRS3035RES1000mN3421000E4743000 0 1 0
## 4 74 CRS3035RES1000mN3421000E4743000 0 1 0
## 5 75 CRS3035RES1000mN3421000E4743000 0 1 0
## 6 76 CRS3035RES1000mN3421000E4743000 0 1 0
## geometry
## 1 POLYGON ((328823.1 655646.1...
## 2 POLYGON ((328810.5 655547, ...
## 3 POLYGON ((328797.9 655447.9...
## 4 POLYGON ((328785.4 655348.7...
## 5 POLYGON ((328772.8 655249.6...
## 6 POLYGON ((328760.2 655150.5...
Waga (oznaczająca jaka część liczby ludności z oczka siatki 1km ma być przypisana do oczka siatki 100m) obliczana jest jako
\[WAGA = \frac{Liczba adresów w siatce 100m}{Liczba adresów w siatce 1km} = \frac{NADRESY100}{NADRESY}\]
$WAGA <- grid100m$N_ADRESY100/grid100m$N_ADRESY grid100m
\[ESTPOP = WAGA * TOT\] EST_POP - Liczba ludności w siatce 100m. TOT - liczba ludności w oczku siatki 1km.
$EST_POP <- grid100m$WAGA *grid100m$TOT grid100m
$EST_POP[is.na(grid100m$EST_POP)]<-0 grid100m
st_write(grid100m, "out/grid100m_pop.shp", delete_dsn = TRUE)
## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver
## Deleting source `out/grid100m_pop.shp' using driver `ESRI Shapefile'
## Writing layer `grid100m_pop' to data source `out/grid100m_pop.shp' using driver `ESRI Shapefile'
## Writing 10500 features with 7 fields and geometry type Polygon.