Wprowadzenie

Projekt obejmuje analizę rozmieszczenia przestrzennego restauracji z platformy pyszne.pl w Warszawie. Celem analizy było zrozumienie, w jaki sposób restauracje z tej platformy są rozlokowane oraz gdzie znajdują się ich największe skupiska. Dodatkowo zestawiono informacje o rozmieszczeniu tych restauracji z rozmieszeniem ludności i mieszkań na mapie Warszawy.

Dane i pakiety

Informacje o restauracjach zostały pozyskane metodą web-scrapingu ze strony pyszne.pl. Dane scrapowano dla wszystkich kodów pocztowych w Warszawie, a następnie na potrzeby tego projektu, dla każdej restauracji zagregowano informacje o jej położeniu geograficznym, średniej ocenie i ich ilości oraz średnich kosztach transportu. Na początku wczytano je i ograniczono do obserwacji z pełnymi danymi o położeniu geograficznym oraz tylko tych położonych w Warszawie. Ostatecznie baza składała się z 2442 restauracji. Dodatkowo akodowano też dane w klasach sf i sp.

r.pyszne <- read.csv("rest_pyszne_full.csv")
head(r.pyszne)
##                                                                                           restaurant
## 1 https://www.pyszne.pl/zakupy/alkohole-24h-warszawa-powstancow-slaskich?sortBy=name&menuItemId=true
## 2                                                     https://www.pyszne.pl/menu/barracuda-sushi-bar
## 3                                                            https://www.pyszne.pl/menu/ciekawa-cafe
## 4                    https://www.pyszne.pl/zakupy/cofesto-ciastko-i-kawa?sortBy=name&menuItemId=true
## 5                                                          https://www.pyszne.pl/menu/indian-summers
## 6                                  https://www.pyszne.pl/menu/kawiarnia-w-pasazu-powstancow-slaskich
##   longitude latitude                      street.address district postal.code
## 1  20.92838 52.26222             Powstańców Śląskich 124   Bemowo      01-466
## 2  20.91331 52.24265                           Kazubów 1   Bemowo      01-466
## 3  20.91337 52.24199             Powstańców Śląskich 80d   Bemowo      01-466
## 4  20.91502 52.23929 Stacja Metro Bemowo lokal 1011 C-04   Bemowo      01-466
## 5  20.91520 52.24019                     Bolkowska 2C/H4   Bemowo      01-466
## 6  20.91962 52.25517        Powstańców Śląskich 106D/210   Bemowo      01-466
##   address.region rating_score rating_votes delivery_to_postals
## 1       Warszawa           NA           NA                 776
## 2       Warszawa          4.9            5                 629
## 3       Warszawa          4.2          124                 217
## 4       Warszawa           NA           NA                 213
## 5       Warszawa          4.8         1835                3435
## 6       Warszawa          4.9            7                 782
##   delivery_to_districts transport_min_cost transport_mean_cost
## 1                     9                  0           11.020558
## 2                     8                 10           14.651163
## 3                     5                  0            4.535337
## 4                     4                  0            5.550000
## 5                    18                  3           17.430398
## 6                     7                  0           11.340150
##   transport_max_cost
## 1              11.99
## 2              20.00
## 3               9.99
## 4               9.99
## 5              25.00
## 6              11.99
# Usunięcie restauracji poza Warszawą
r.pyszne <- r.pyszne[r.pyszne$address.region == "Warszawa",]
r.pyszne <- r.pyszne[!r.pyszne$district == "",]
unique(r.pyszne$district)
##  [1] "Bemowo"         "Włochy"         "Wola"           "Ursus"         
##  [5] "Praga-Południe" "Bielany"        "Targówek"       "Śródmieście"   
##  [9] "Praga-Północ"   "Ursynów"        "Mokotów"        "Żoliborz"      
## [13] "Wilanów"        "Ochota"         "Białołęka"      "Wesoła"        
## [17] "Wawer"          "Rembertów"
# Usunięcie restauracji z brakiem lokalizacji
r.pyszne <- r.pyszne[complete.cases(r.pyszne[,c(2:3)]),]

dim(r.pyszne) 
## [1] 2442   14
library(sp)
library(spdep)
library(sf)

# Restauracje klasa SP
r.pyszne.sp <- r.pyszne     
coordinates(r.pyszne.sp) <-c ("longitude", "latitude") 
proj4string(r.pyszne.sp) <- as(st_crs("+proj=longlat +ellps=WGS84"), "CRS")
summary(r.pyszne.sp) 

# Restauracje klasa SF
r.pyszne.sf <- st_as_sf(r.pyszne.sp)
r.pyszne.sf
r.pyszne.sf <- st_transform(r.pyszne.sf, 4326)

Następnie wczytano pliki shapefile z informacjami o konturach Warszawy i jej dzielnic oraz zapisono je klasach sf i sp.

########## KONTURY DZIELNIC ##########
dzielnice.sf <- st_read("dzielnice_Warszawy.shp", stringsAsFactors = FALSE)
dzielnice.sf <- st_transform(dzielnice.sf, 4326) 
dzielnice.sp <- as_Spatial(dzielnice.sf, cast = TRUE, IDs ="nazwa_dzie")

########## KONTURY WARSZAWY ##########
powiaty <- st_read("powiaty.shp", stringsAsFactors = FALSE)
powiaty <- st_transform(powiaty, 4326)
warszawa.sf <- powiaty[powiaty$jpt_nazwa == "powiat Warszawa",]
warszawa.sp <- as_Spatial(warszawa.sf, cast = TRUE, IDs ="jpt_kod_je")

Rozmieszczenie restauracji i zakres dowozów.

W pierwszej kolejności zwizualizowano jak rozmieszczone są restauracje na mapie Warszawy.

library(ggplot2)

ggplot() + geom_sf(data = dzielnice.sf, fill = "white") + 
  theme_minimal() + 
  geom_point(data = r.pyszne, aes(x = longitude, y = latitude), color = "#66295F", size = 2, alpha = 0.3) + 
  ggtitle("Restauracje pyszne.pl") +
  theme(plot.title = element_text(hjust = 0.5, size = 20),
        axis.title.x = element_blank(),
        axis.title.y = element_blank())

Patrząc na mapę, zauważyć można iż restauracje nie są do końca równomiernie rozłożone - najwięcej jest ich w centralnych dzielnicach. Zliczmy ile jest restauracji dla każdej z dzielnic i zwizualizujmy to.

library(dplyr) 
# Liczba restauracji w każdej dzielnicy
liczba.restauracji <- r.pyszne %>% 
  group_by(district) %>% summarise(number.of.restaurants = n()) %>% ungroup
dzielnice.sf <- merge(dzielnice.sf, liczba.restauracji,
                      by.x = "nazwa_dzie", by.y = "district", all.x = TRUE, sort = FALSE)
centroids <- st_centroid(dzielnice.sf)
ggplot() +
  geom_sf(data = dzielnice.sf, 
          aes(fill = number.of.restaurants), color = "gray5") +
  scale_fill_gradient(low = "#EED3EA", high = "#66295F", name = "") +
  theme_minimal() +
  ggtitle("Restauracje w poszczególnych dzielnicach") +
  theme(plot.title = element_text(hjust = 0.5, size = 15),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  geom_sf_text(data = centroids, aes(label = nazwa_dzie), color = "black", size = 2)

Widoczne jest iż najwięcej jest ich na Śródmieściu, Woli i Mokotowie. Spójrzmy też na to jak kształtuje się liczba kodów do których dowozi restauracja.

ggplot() + 
  geom_sf(data = dzielnice.sf, fill = "white") + 
  theme_minimal() + 
  geom_point(data = r.pyszne,
             aes(x = longitude, y = latitude, color = delivery_to_postals), size = 2, alpha = 0.5) + 
  ggtitle("Liczba adresów pocztowych do których dowozi restauracja")+
  theme(plot.title = element_text(hjust = 0.5, size = 15),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  scale_color_viridis_c(option = "A", name = "")

W centrum Warszawy kropki mają jaśniejszy kolor - tutaj oznacza to, że restauracje dowożą do większej liczby adresów pocztowych, porównując z ciemniejszymi kropkami - te restauracje dowożą do mniejszej liczby adresów pocztowych.

Warszawa w siatce grid - restauracje i populacja

Połączmy informacje o restauracjach z danymi dotyczącymi liczby populacji jak i budynków z NSP 2021, pochodzących z GUS1. W pierwszej kolejności ładujemy dane, wycinamy grid jedynie dla Warszawy i łączymy dane o budynkach i populacji w jeden df. Dla Warszawy finalnie mamy grid składający się z 601 kratek wielkości 1km x 1km. Dla każdej kratki mamy informację o liczbie rezydentów ogółem, w podziale na płeć oraz na wiek jak i liczbę budynków oraz mieszkań.

########## Grid populacja ##########
POP.2021 <- st_read("trkx19307.shp", stringsAsFactors=FALSE)
## Reading layer `trkx19307' from data source 
##   `/Users/monikakot/Desktop/Spatial/Projekt_spatial/trkx19307.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 315856 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 14.11547 ymin: 48.99754 xmax: 24.1529 ymax: 54.84028
## Geodetic CRS:  WGS 84
# wycięcie Warszawy
mer.pop <- st_join(POP.2021, warszawa.sf, join = st_intersects)
a <- which(mer.pop$jpt_nazwa_ == "powiat Warszawa")
POP.2021.waw <- POP.2021[a,]
POP.2021.waw
## Simple feature collection with 601 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 20.83728 ymin: 52.09162 xmax: 21.2805 ymax: 52.378
## Geodetic CRS:  WGS 84
## First 10 features:
##       shape_leng shape_area                            code
## 63632   3997.988   998991.0 CRS3035RES1000mN3286000E5059000
## 63660   3998.006   999000.1 CRS3035RES1000mN3298000E5059000
## 63672   3997.989   998991.6 CRS3035RES1000mN3287000E5059000
## 63688   3998.008   999000.9 CRS3035RES1000mN3299000E5059000
## 63701   3998.005   998999.3 CRS3035RES1000mN3297000E5059000
## 63798   3998.009   999001.7 CRS3035RES1000mN3300000E5059000
## 64015   3997.998   998996.2 CRS3035RES1000mN3293000E5059000
## 64197   3997.991   998992.4 CRS3035RES1000mN3288000E5059000
## 64223   3998.016   999004.8 CRS3035RES1000mN3296000E5060000
## 64235   3998.005   998999.4 CRS3035RES1000mN3289000E5060000
##                                  pl_code  res res_male res_fem res_0_14
## 63632 PL_CRS3035RES1000mN3286000E5059000 4691     2251    2440      719
## 63660 PL_CRS3035RES1000mN3298000E5059000  156       75      81       13
## 63672 PL_CRS3035RES1000mN3287000E5059000  296      142     154       48
## 63688 PL_CRS3035RES1000mN3299000E5059000    1       NA      NA        0
## 63701 PL_CRS3035RES1000mN3297000E5059000  423      203     220       64
## 63798 PL_CRS3035RES1000mN3300000E5059000  134       56      78       18
## 64015 PL_CRS3035RES1000mN3293000E5059000 1367      647     720      226
## 64197 PL_CRS3035RES1000mN3288000E5059000  121       51      70       14
## 64223 PL_CRS3035RES1000mN3296000E5060000    0        0       0        0
## 64235 PL_CRS3035RES1000mN3289000E5060000  113       53      60       21
##       res_15_64 res_65_ multipolyg multipoly1                       geometry
## 63632      2873    1099    2656647   6519.691 MULTIPOLYGON (((20.85177 52...
## 63660       114      29    2669361   6535.273 MULTIPOLYGON (((20.87799 52...
## 63672       196      52    2657702   6520.985 MULTIPOLYGON (((20.85395 52...
## 63688        NA      NA    2670426   6536.576 MULTIPOLYGON (((20.88018 52...
## 63701       274      85    2668298   6533.971 MULTIPOLYGON (((20.8758 52....
## 63798        74      42    2671491   6537.879 MULTIPOLYGON (((20.88238 52...
## 64015       875     266    2664050   6528.768 MULTIPOLYGON (((20.86705 52...
## 64197        89      18    2658759   6522.281 MULTIPOLYGON (((20.85613 52...
## 64223         0       0    2667074   6532.473 MULTIPOLYGON (((20.88812 52...
## 64235        68      24    2659656   6523.382 MULTIPOLYGON (((20.8728 52....
########## Grid budynki ##########
BUD.2021 <- st_read("twmn18483.shp", stringsAsFactors=FALSE)
## Reading layer `twmn18483' from data source 
##   `/Users/monikakot/Desktop/Spatial/Projekt_spatial/twmn18483.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 315856 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 14.11547 ymin: 48.99754 xmax: 24.1529 ymax: 54.84028
## Geodetic CRS:  WGS 84
# wycięcie Warszawy
mer.bud <- st_join(BUD.2021, warszawa.sf, join = st_intersects)
a <- which(mer.bud$jpt_nazwa_ == "powiat Warszawa")
BUD.2021.waw <- BUD.2021[a,]
BUD.2021.waw
## Simple feature collection with 601 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 20.83728 ymin: 52.09162 xmax: 21.2805 ymax: 52.378
## Geodetic CRS:  WGS 84
## First 10 features:
##        shape_leng shape_area                            code
## 221862   3997.988   998991.0 CRS3035RES1000mN3286000E5059000
## 221890   3998.006   999000.1 CRS3035RES1000mN3298000E5059000
## 221902   3997.989   998991.6 CRS3035RES1000mN3287000E5059000
## 221918   3998.008   999000.9 CRS3035RES1000mN3299000E5059000
## 221931   3998.005   998999.3 CRS3035RES1000mN3297000E5059000
## 222028   3998.009   999001.7 CRS3035RES1000mN3300000E5059000
## 222245   3997.998   998996.2 CRS3035RES1000mN3293000E5059000
## 222428   3997.991   998992.4 CRS3035RES1000mN3288000E5059000
## 222454   3998.016   999004.8 CRS3035RES1000mN3296000E5060000
## 222466   3998.005   998999.4 CRS3035RES1000mN3289000E5060000
##                                   pl_code bud dwel multipolyg multipoly1
## 221862 PL_CRS3035RES1000mN3286000E5059000 709 1924    2656647   6519.691
## 221890 PL_CRS3035RES1000mN3298000E5059000  40   43    2669361   6535.273
## 221902 PL_CRS3035RES1000mN3287000E5059000  95  101    2657702   6520.985
## 221918 PL_CRS3035RES1000mN3299000E5059000   1    1    2670426   6536.576
## 221931 PL_CRS3035RES1000mN3297000E5059000 130  144    2668298   6533.971
## 222028 PL_CRS3035RES1000mN3300000E5059000  51   50    2671491   6537.879
## 222245 PL_CRS3035RES1000mN3293000E5059000 307  620    2664050   6528.768
## 222428 PL_CRS3035RES1000mN3288000E5059000  32   46    2658759   6522.281
## 222454 PL_CRS3035RES1000mN3296000E5060000   0    0    2667074   6532.473
## 222466 PL_CRS3035RES1000mN3289000E5060000  39   48    2659656   6523.382
##                              geometry
## 221862 MULTIPOLYGON (((20.85177 52...
## 221890 MULTIPOLYGON (((20.87799 52...
## 221902 MULTIPOLYGON (((20.85395 52...
## 221918 MULTIPOLYGON (((20.88018 52...
## 221931 MULTIPOLYGON (((20.8758 52....
## 222028 MULTIPOLYGON (((20.88238 52...
## 222245 MULTIPOLYGON (((20.86705 52...
## 222428 MULTIPOLYGON (((20.85613 52...
## 222454 MULTIPOLYGON (((20.88812 52...
## 222466 MULTIPOLYGON (((20.8728 52....
########## Grid populacja + budynki ##########
ALL.2021.waw <- st_join(POP.2021.waw, BUD.2021.waw, join = st_equals)

Następnie łączymy dane grid z informacjami o restauracjach. Najpierw do każdej restauracji przyporządkowujemy kratkę 1km x 1km w której leży, a następnie zliczmy ile jest restauracji w każdej z nich - dzięki temu mamy jednolite dane.

########## Połączenie gridu + restauracji ##########
ALL.2021.waw$myID <- 1:dim(ALL.2021.waw )[1]
mer.res <- st_join(r.pyszne.sf, ALL.2021.waw, join = st_intersects)
dim(mer.res) # do każdej restauracji przyporządkowana kratka
## [1] 2442   34
########## Zliczenie restauracji w każdej kratce ##########
restaurants.agg <- mer.res %>% group_by(myID) %>% summarise(number.of.restaurants = n())
restaurants.agg <- as.data.frame(restaurants.agg)
restaurants.agg <- restaurants.agg[, 1:2]
ALL.2021.waw <- merge(ALL.2021.waw, restaurants.agg,
                      by.x = "myID", by.y = "myID", all.x = TRUE, sort = FALSE)
# Dla kratek z brakiem przyporządkowanych restauracji zapisujemy 0
ALL.2021.waw$number.of.restaurants[is.na(ALL.2021.waw$number.of.restaurants)] <- 0

Zwizualizujmy teraz jak rozmieszczone są restauracje, populacja i mieszkania w Warszawie.

Patrząc na te trzy mapy, widać iż kształt rozmieszczenia resturacji jest podobny do rozmieszczenia ludności i mieszkań w Warszawie. Prosta regresja liniowa potwierdza iż liczba restauracji jest zależna od liczby populacji i mieszkań.

mod <- lm(number.of.restaurants~res+dwel, data = ALL.2021.waw)
summary(mod)
## 
## Call:
## lm(formula = number.of.restaurants ~ res + dwel, data = ALL.2021.waw)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.355  -0.740  -0.344   0.382  44.024 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.3695006  0.2740887   1.348    0.178    
## res         -0.0031591  0.0002513 -12.571   <2e-16 ***
## dwel         0.0078284  0.0004215  18.573   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.281 on 598 degrees of freedom
## Multiple R-squared:  0.6543, Adjusted R-squared:  0.6531 
## F-statistic: 565.8 on 2 and 598 DF,  p-value: < 2.2e-16

W dalszej części projektu sprawdzimy czy w rozmieszczeniu restauracji z platformy pyszne.pl można stwierdzić aglomeracje oraz jakie klastry tworzą restauracje.

ETA (Entropy-Tesselation-Agglomeration)

Entropia, czyli tak zwana miara nieporządku, jest jedną z metod statystycznych pozwalających uchwycić zjawisko aglomeracji. Dostarcza ona poniekąd informacji na temat prawdopodobieństwa znalezienia się w danej lokalizacji - gdy entropia jest wysoka to punkty są bardzo regularnie rozłożone, więc ciężko jest oszacować gdzie znajdzie się kolejna lokalizacja - prawdopodobieństwo dla każdego z nich jest podobne. Z kolei gdy entropia jest niska, to rozłożenie punktów jest nieregularne i z większym pradopodobieństwem można stwierdzić położenie nowych punktów. W tej metodzie entropia zostala połączona z teselacją, która pozwala na przejście z danych punktowych na obszary podzielone na tzw. kafelki, których różne rozmiary, w zależności od stopnia koncentracji danych, pełnią rolę prawodpodobieństw. W przypadku równomiernego rozkładu, rozmiar poszczególnych kafelków jest praktycznie identyczny, a w przypadku aglomeracji zróżnicowany - w obszarze gdzie punktów jest dużo to kafelki będą małe.

Aglomerację można zmierzyć za pomocą porównania maksymalnej, hipotetycznej entropii z tą zaobserwowaną, uzyskując tak zwaną relatywną entropię. Uzyskana przez nas wartość relatywna entropii mniejsza niż jeden oznacza, iż badany rozkład nie jest jednostajny i można wyróżnic aglomarację, co jest widoczne w centralnych częściach Warszawy.

Quick density clustering

Ta metoda wykorzystuje dwie zmienne przestrzenne - całkowity dystans do k najbliższych sąsiadów oraz liczbę sąsiadów w stałym promieniu. Znormalizowane zmienne są następnie poddawane klasteryzjacji metodą k-średnich. Parametry zostały ustalone na liczenie odległości do 10 sąsiadów i zliczanie sąsiadów w promieniu 0.05.

## [1] 1.970513
## [1] 0.1010681

Na wykresie widzimy podział restauracji na trzy grupy: o małym, średnim i dużym zagęszczeniu. Te w grupie o najmniejszej gęstości charakteryzują się najwyższymi odleglościami do k najbliższych sąsiadów oraz najniższą liczbą sąsiadów w promieniu. Z kolei te o największej gęstości mają najmniejsze odległości do k najbliższych sąsiadów oraz największą liczbę sąsiadów w promieniu. Zwizualizujmy te grupy na mapie Warszawy.

Grupa restauracji o największej gęstości znajduje się na Śródmieściu, Woli, Ochocie, Żoliborzu i części Mokotowa. Z kolei grupa o najmniejszej gęstości znajduje się na dzielnicach bardziej obrzeżnych - na Białolęce, Rembertowie czy Wawrze.

Poszukiwanie klastrów metodą DBSCAN

Metoda ta wykorzystuje dwa parametry: promień koła oraz minimalną liczbę punktów oczekiwanych w tym kole. Spełnienie tego warunku decyduje o przynależności danego punktu do klastra o wysokiej gęstości. W tej metodzie przy grupowaniu danych nie zakłada się oczekiwanej liczby klastrów, tylko ustala poziom parametrów. W tym przypadku wyznaczono promień na 0.01, a oczkiwaną liczbę punktów (restauracji) w danym promieniu na 50.

Algorytm wyznaczył 7 klastrów - największy z nich obejmuje większą część Śródmieścia oraz część Woli, Ochoty oraz Mokotowa. Wyodrębniono też oddzielne klastry, m.in. na południowej części Mokotowa, na Pradze Południe czy Bemowie.

Focal local entropy

Metoda ta analizuje różnorodność wartości w sąsiedztwie każdego punktu w przestrzeni. Wysokie wartości lokalnej entropii świadczą o występowaniu bardziej zróżnicowanych danych, a niskie o ich podobieństwie. Poniższe mapy prezentują wartości lokalnej entropii dla rozkładu przestrzennego restauracji po przekształceniu danych punktowych w rastry składające się z 10 000 komórek (100 wierszy x 100 kolumn), dla różnych wartości ruchomego okna.

Dla mniejszych wartości okna efekty koncentracji przestrzennej wydają się być mniej wyraźne, co wskazuje na lepsze dopasowanie dla wyższych jego wartości. W obszarze centrum Warszawy oraz jego okolicach wartości lokalnej entropii są relatywnie wysokie, co sugeruje, że obszar pod względem restauracji jest bardziej zrożnicowany. Czym dalej od centrum, tym wartości lokalnej entropii są niższe, co wskazuje na coraz mniejsze zróżnicowanie pod względem restauracji wraz z oddalaniem się od centralnych dzielnic.

Podsumowanie

Przeprowadzona analiza eklsporacyjna pozwoliła na wyciągnięcie kilku wniosków dotyczących przestrzennego rozmieszczenia restauracji z platformy pyszne.pl w Warszawie. Po pierwsze rozmieszczenie to nie jest równomierne - resturacja najgęściej skupiają się w centralnych dzielnicach miasta. Taki rozkład przestrzenny koreluje z ogólnym rozkładem ludności oraz mieszkań w mieście,co zostało potwierdzone za pomocą prostej regresji liniowej. Po drugie, za pomocą kilku różnych metod, udało się wyróżnić aglomeracje. Każda z nich dostarczyła podobnych wniosków - im dalej od centrum Warszawy, tym mniej gęsto skupione są restauracje.

  1. https://portal.geo.stat.gov.pl/↩︎