Wykorzystanie modelowania dazymetrycznego do opracowania szczegółowej (100m) mapy rozmieszczenia ludności dla zlewni górnej Parsęty
Część 2. Wykorzystanie analizy regresji do określenia zależności między ogólną liczbą ludności a liczbą mieszkań.
Anna Dmowska (dmowska@amu.edu.pl)
library(sf) #analizy przestrzenna danych wektorowych w R## Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3library(rgugik) #pobranie danych BDOT10k z GUGiK
library(ggplot2)
library(ggResidpanel) #interaktywne wykresy diagnostyczne
library(ggfortify) # wykresy diagnostyczne w stylizacji ggplot2.1 Cel ćwiczenia
- Celem analizy jest sprawdzenie czy istnieje zależność między liczbą ludności oraz liczbą mieszkań?
- W przypadku odnotowania takiej zależności liczba mieszkań mogłyby być wykorzystane do dekompozycji ogólnej liczby ludności w celu otrzymania bardziej szczegółowej mapy (np. o rozdzielczości 100m).
2 Dane
2.1 Granica zlewni górnej Parsęty
- Źródło danych: MPHP (Mapa Podziału Hydrograficznego Polski).
- Układ współrzędnych: Polski Układ Współrzędnych Geograficznych 1992 (EPSG: 2180)
zlewnia <- st_read("dane_dasy/zlewnia_gp.shp")## 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_defs2.2 Rozmieszczenie ludności w siatce 1x1km
- Źródło danych: GRID 2011 – rozmieszczenie ludności w siatce kilometrowej dla roku 2011. Dane zostały opracowane w oparciu o wyniki Narodowego Spisu Powszechnego Ludności i Mieszkań 2011.
- Pobranie danych: https://geo.stat.gov.pl/inspire
- Układ współrzędnych: EPSG:3035
#wczytanie danych z pliku shp
pop1km_3035 <- st_read("dane_dasy/pop1km_gp3035.shp")## 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
pop1km <- st_transform(pop1km_3035, 2180)
pop1km <- pop1km[,c("TOT", "CODE")]
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...2.3 Siatka o rozdzielczości 100m.
Siatka o wielkości oczka 100 na 100m posłuży do przygotowania końcowej mapy rozmieszczenia ludności. Przygotowanie takiej siatki jest łatwiejsze w QGIS.
grid100m <- st_read("dane_dasy/grid100m_2180.shp")## 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_defsgrid100m <- grid100m[,1]3 Analiza zależności między ogólną liczbą ludności, a liczbą mieszkań.
W poprzednim przykładzie sprawdzaliśmy, czy istnieje zależność między liczbą ludności a liczbą punktów adresowych. W tym przykładzie sprawdzimy zależność między liczbą ludności a liczbą mieszkań.
Etapy analizy
Etapy analizy obejmują:
- przygotowanie danych do analizy.
- oszacowanie liczby mieszkań.
- zbudowanie modelu liniowego między liczbą ludności oraz liczbą punktów adresowych.
- ocena dopasowania modelu.
- wizualizacja wyników.
3.1 Pobranie danych dotyczących budynków.
- Źródło danych: Baz Danych Obiektów Topograficznych (BDOT10k) https://www.geoportal.gov.pl/dane/baza-danych-obiektow-topograficznych-bdot
- Dane dla powiatu pobrane wykorzystując pakiet R
rgugikhttps://cran.r-project.org/web/packages/rgugik/vignettes/topodb.html - Układ współrzędnych: Polski Układ Współrzędnych Geograficznych 1992 (EPSG: 2180)
library(rgugik)
topodb_download("szczecinecki", outdir = "./dane_dasy") #POBIERANIE DANYCH - WYMAGA DOSTEPU DO INTERNETU
budynki <- read_sf("./dane_dasy/PL.PZGiK.339.3215/BDOT10k/PL.PZGiK.339.3215__OT_BUBD_A.xml")
head(budynki)## Simple feature collection with 6 features and 36 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 325520.7 ymin: 655642.3 xmax: 333559 ymax: 670101.9
## CRS: 2180
## # A tibble: 6 x 37
## gml_id lokalnyId przestrzenNazw wersjaId czyObiektBDOO x_kod x_skrKarto
## <chr> <chr> <chr> <chr> <lgl> <chr> <chr>
## 1 OT_BUBD… 1FDCBD00-AF1… PL.PZGiK.339.… 2021-02-… FALSE BUBD… <NA>
## 2 OT_BUBD… 1FDCBD00-AF1… PL.PZGiK.339.… 2021-02-… FALSE BUBD… <NA>
## 3 OT_BUBD… 1FDCBD00-AF2… PL.PZGiK.339.… 2021-02-… FALSE BUBD… <NA>
## 4 OT_BUBD… 1FDCBD00-AF2… PL.PZGiK.339.… 2021-02-… FALSE BUBD… <NA>
## 5 OT_BUBD… 1FDCBD00-B1B… PL.PZGiK.339.… 2021-02-… FALSE BUBD… <NA>
## 6 OT_BUBD… 1FDCBD00-B9F… PL.PZGiK.339.… 2021-02-… FALSE BUBD… <NA>
## # … with 30 more variables: x_katDoklGeom <chr>, x_doklGeom <int>,
## # x_doklGeom_uom <chr>, x_zrodloDanychG <chr>, x_zrodloDanychA <chr>,
## # x_katIstnienia <chr>, x_rodzajReprGeom <chr>, x_aktualnoscG <chr>,
## # x_aktualnoscA <chr>, poczatekWersjiObiektu <chr>, x_dataUtworzenia <chr>,
## # x_kodKarto10k <chr>, x_kodKarto25k <chr>, x_kodKarto50k <chr>,
## # x_kodKarto100k <chr>, x_kodKarto250k <chr>, x_kodKarto500k <chr>,
## # x_kodKarto1000k <chr>, funOgolnaBudynku <int>, liczbaKondygnacji <int>,
## # kodKst <int>, zabytek <lgl>, x_informDodatkowa <chr>, nazwa <chr>,
## # x_uwagi <chr>, funSzczegolowaBudynku <list>,
## # EGiB|BT_ReferencjaDoObiektu|idIIP|BT_Identyfikator|lokalnyId <list>,
## # EGiB|BT_ReferencjaDoObiektu|idIIP|BT_Identyfikator|przestrzenNazw <list>,
## # EGiB|BT_ReferencjaDoObiektu|idIIP|BT_Identyfikator|wersjaId <list>,
## # geometry <POLYGON [m]>budynki <- budynki[,c("x_kod", "liczbaKondygnacji")]
colnames(budynki) <- c("TYP_BUD", "KONDYGNACJE", "geometry")3.2 Przygotowanie danych do analizy.
3.2.1 Wyznaczenie budynków mieszkalnych
Warstwa zawiera 4 typy budynków mieszkalnych (zmienna TYP_BUD): - “BUBD01” - budynki jednorodzinne, - “BUBD02” - budynki z dwoma mieszkaniami, - “BUBD03” - budynki z trzema i więcej mieszkaniami, - “BUBD04” - budynki zamieszkania zbiorowego (w tym bursy, akademiki, klasztory)
budynki_mieszkalne <- budynki[budynki$TYP_BUD%in%c("BUBD01", "BUBD02", "BUBD03", "BUBD04"),] Liczba kondygnacji dla poszczególnych typów budynków.
table(budynki_mieszkalne$TYP_BUD, budynki_mieszkalne$KONDYGNACJE)##
## 1 2 3 4 5 6 7
## BUBD01 4435 4790 237 2 0 0 0
## BUBD02 18 27 3 0 0 0 0
## BUBD03 87 417 404 112 217 1 2
## BUBD04 13 23 9 2 1 0 0Selekcja budynków mieszkalnych dla obszaru analizy
Selekcję obiektów na podstawie lokalizacji można wykonać wykorzystując funkcje st_join() lub też indeksowanie (nawiasy kwadratowe).
budynki_mieszkalne_pop1km <- budynki_mieszkalne[pop1km,]plot(budynki_mieszkalne_pop1km)
st_write(budynki_mieszkalne_pop1km, "out/budynki_mieszkalne_pop1km.shp", delete_dsn = TRUE)## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver## Deleting source `out/budynki_mieszkalne_pop1km.shp' using driver `ESRI Shapefile'
## Writing layer `budynki_mieszkalne_pop1km' to data source `out/budynki_mieszkalne_pop1km.shp' using driver `ESRI Shapefile'
## Writing 568 features with 2 fields and geometry type Polygon.Na obszarze zlewni górnej Parsęty występują dwa typy budynków: budynki jednorodzinne (“BUBD01”) oraz budynki z trzema i więcej mieszkaniami (“BUBD03”)
3.2.2 Oszacowanie liczby mieszkań.
Liczba kondygnacji dla poszczególnych typów budynków.
table(budynki_mieszkalne_pop1km$TYP_BUD, budynki_mieszkalne_pop1km$KONDYGNACJE)##
## 1 2 3
## BUBD01 389 159 4
## BUBD03 2 9 5Zakładamy, że budynki jednorodzinne składają się z 1 mieszkania. Dla typów “BUBD03” (budynki z trzema i więcej mieszkaniami) zakładamy, że liczba mieszkań równa jest liczbie kondygnacji w budynku. Takie założenie może niedoszacować liczby mieszkań w budynkach wielokondygnacyjnych, w szczególności w obszarach miejskich.
budynki_mieszkalne_pop1km$KATEGORIA[budynki_mieszkalne_pop1km$TYP_BUD%in%c("BUBD01")]<-1## Warning: Unknown or uninitialised column: `KATEGORIA`.#budynki_mieszkalne_pop1km$KATEGORIA[budynki_mieszkalne_pop1km$TYP_BUD%in%c("BUBD02")]<-2
budynki_mieszkalne_pop1km$KATEGORIA[budynki_mieszkalne_pop1km$TYP_BUD%in%c("BUBD03", "BUBD04")]<- budynki_mieszkalne_pop1km$KONDYGNACJE[budynki_mieszkalne_pop1km$TYP_BUD%in%c("BUBD03", "BUBD04")]3.2.3 Obliczenie liczby mieszkań w oczkach siatki 1km.
W przypadku punktów adresowych każdy punk kontrybuował z wagą 1. W przypadku liczby mieszkań, dany budynek będzie kontrybuował z określoną wagą: 1 - dla budynków jednorodzinnych, 2 - dla budynków z dwoma mieszkaniami, 3 i więcej (w zależności od liczby kondygnacji utożsamianej tu z liczbą mieszkań) - dla budynków z 3 lub więcej mieszkaniami.
centroidy_budynki <- st_centroid(budynki_mieszkalne_pop1km)## Warning in st_centroid.sf(budynki_mieszkalne_pop1km): st_centroid assumes
## attributes are constant over geometries of xbudynki_w_pop1km <- st_join(centroidy_budynki, pop1km[,c("TOT", "CODE")])count_budynki <- aggregate(KATEGORIA~CODE, budynki_w_pop1km, sum, na.rm=TRUE)
colnames(count_budynki) <- c("CODE", "N_BUDYNKI")pop1km <- merge(pop1km, count_budynki, by = "CODE", all.x = TRUE)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_BUDYNKI geometry
## 1 CRS3035RES1000mN3418000E4746000 0 NA POLYGON ((332289.2 651182.4...
## 2 CRS3035RES1000mN3418000E4747000 55 14 POLYGON ((333281.4 651057.9...
## 3 CRS3035RES1000mN3419000E4744000 61 15 POLYGON ((330430.4 652422.9...
## 4 CRS3035RES1000mN3419000E4745000 38 15 POLYGON ((331422.6 652298.3...
## 5 CRS3035RES1000mN3419000E4746000 16 8 POLYGON ((332414.9 652173.8...
## 6 CRS3035RES1000mN3419000E4747000 0 NA POLYGON ((333407.1 652049.2...pop1km$N_BUDYNKI[is.na(pop1km$N_BUDYNKI)]<- 0