Część teoretyczna

Ta część obecnie jest w budowie…

Abstrakt

Ta część obecnie jest w budowie…

Wstęp

STRONA PROWADZONA W TRYBIE ROLLING I BĘDZIE STOPNIOWO AKTUALIZOWANA, WRAZ Z POSTAMI NA LINKEDIN

Warszawa dysponuje jedną z najgęstszych sieci punktów odbioru paczek w Polsce. Na mapie Warszawy znajduje się 3148 automatów paczkowych (APM) oraz 3223 punktów odbioru w punktach usługowych (PUDO), firm takich jak: InPost, DPD, Orlen, Poczta Polska oraz DHL. Łącznie stanowi to 6371 punktów odbioru paczek (OOH). Niewątpliwie, jest to związane z rozwojem sieci e-commerce [POPRAW]. Wszechobecna sieć automatów paczkowych rodzi jednak kluczowe pytanie:

  • czy automaty lokalizowane są czysto losowe?
  • czy istnieje przestrzenna zależność między automatami?
  • czy firmy preferują określone lokalizacje? W których miejscach firmy mają rynkową dominacje nad innymi firmami?
  • od czego zależy alokacja automatów przez dane firmy?

Na te, oraz wiele innych pytań, odpowiem podczas przygotowanej analizy. Zapraszam do śledzenia RPubs, na którym można zobaczyć od przysłowiowej “kuchni” jak wygląda proces analizy wzorca punktowego na przykładzie paczkomatów.

Niniejsza praca służy w celach czysto edukacyjnych w zakresie użycia analizy wzorca punktowego i liniowego (Point Pattern Analysis, PPA) przy użycia języka R. Do analizy wykorzystano publicznie dostępne dane na stronie dostawców (stan na kwiecień 2026).

Teoria

Ta część obecnie jest w budowie…

Kod

Kod przygotowany na podstawie zajęć Analiza wzorca punktowego i liniowego w R, prowadzonego przez dr Katerynę Zabarinę.

Wczytujemy biblioteki

Wczytajmy wszystkie potrzebne biblioteki

requiredPackages<-c("sf", "spatstat", "rgugik", "automap",
                      "magrittr", "sp", "rnaturalearth",
                      "ggplot2", "tidyverse", "viridis")
for(i in requiredPackages)
  {if(!require(i, character.only=TRUE)) install.packages(i)}
for(i in requiredPackages)
  {library(i, character.only=TRUE)}

setwd("~/AnalizaR/PACZKOMATY")
Sys.setenv(LANG="en")
windowsFonts(TNR=windowsFont("Times New Roman"))

Pobieramy mapę Warszawy oraz przygotowujemy paczkomaty

W tej części:

  • Pobierzemy mapę Warszawy oraz przygotujemy jednostki
  • Wczytamy dane dotyczące paczkomatów
  • Sprawdzimy ile jest paczkomatów poszczególnych marek w Warszawie
  • Stworzymy mapę paczkomatów w Warszawie
#Pobieramy mapę Warszawy
warszawa_granica<-suppressWarnings(borders_get(TERYT="1465")) #Pobieramy mapę Warszawy z GUGiK
warszawa_granica<-st_transform(warszawa_granica, crs=2180) #Ustawiamy CRS (układ planetarny) na 2180, czyli odpowiedni dla Warszawy
st_crs(warszawa_granica) #Sprawdzamy czy jest CRS 2180 po zmianie
st_crs(warszawa_granica)$units_gdal #W jakich jednostkach jest podana mapa?

W<-as.owin(warszawa_granica)
class(W) #Sprawdzamy jaki to obiekt
area_km2<-as.numeric(st_area(warszawa_granica))/1e6 #Powierzchnia Warszawy w km²
rsc<-sqrt(area.owin(W)/area_km2) #Zmieniamy z metrów na kilometry

#Wczytujemy zbiór danych z paczkomatami
paczkomaty<-read.csv("Lista paczkomatow_automaty.csv", fileEncoding="UTF-8-BOM") #Wczytujemy dane paczkomatów z pliku CSV

paczkomaty_sf<-st_as_sf(paczkomaty, coords=c("lon", "lat"), crs=4326) #Kodujemy koordynaty paczkomatów jako długość i szerokość geograficzną
paczkomaty_sf<-st_transform(paczkomaty_sf, crs=2180) #Przechodzimy z CRS 4326 na 2180
st_crs(paczkomaty_sf) #Sprawdzamy czy jest CRS 2180

paczkomaty_waw<-st_intersection(paczkomaty_sf, warszawa_granica) #Przycinamy dataset do Warszawy
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

Sprawdźmy ile mamy automatów paczkowych w zbiorze danych z podziałem na marki

table(paczkomaty_waw$marka)
## 
##           DHL           DPD        INPOST         ORLEN POCZTA POLSKA 
##           417           848          1413           419            51
p_mapa<-ggplot() +
  geom_sf(data=warszawa_granica, fill="grey95", color="black", linewidth=0.6) +
  geom_sf(data=paczkomaty_waw, aes(color=marka), size=1.5, alpha=0.7) +
  scale_color_viridis_d(option="turbo", name="Marka") +
  labs(
    title="Paczkomaty w Warszawie",
    subtitle=paste("Łącznie", nrow(paczkomaty_waw), "punktów")
  ) +
  theme_minimal() +
  theme(plot.title=element_text(size=14, face="bold"))
print(p_mapa)

Marked oraz unmarked PPP

Wyświetlmy podsumowanie wzorca unmarked, czyli takiego, który nie różnicuje punktów ze względu na markę

#Unmarked - wszystkie paczkomaty bez rozróżnienia marki
paczkomaty.ppp.um<-ppp(
  x=st_coordinates(paczkomaty_waw)[, 1], #Współrzędne X
  y=st_coordinates(paczkomaty_waw)[, 2], #Współrzędne Y
  window=W) #Punkty w Warszawie

paczkomaty.ppp.um<-rjitter(paczkomaty.ppp.um, 0.05) #Przesuwamy lekko punkty żeby nie było duplikatów (o 0.05 metra)
paczkomaty.ppp.um<-rescale.ppp(paczkomaty.ppp.um, rsc, "km") #Zmieniamy z metrów na kilometry
summary(paczkomaty.ppp.um) #Wyświetlmy podsumowanie wzorca unmarked
## Planar point pattern:  3148 points
## Average intensity 6.092221 points per square km
## 
## Coordinates are given to 13 decimal places
## 
## Window: polygonal boundary
## single connected closed polygon with 3389 vertices
## enclosing rectangle: [626.506, 655.2603] x [472.2295, 502.1724] km
##                      (28.75 x 29.94 km)
## Window area = 516.724 square km
## Unit of length: 1 km
## Fraction of frame area: 0.6

Wyświetlmy też podsumowanie wzorca marked, czyli takiego, który różnicuje punkty ze względu na markę

#Marka jako cecha paczkomatu
paczkomaty.ppp.m<-ppp(
  x=st_coordinates(paczkomaty_waw)[, 1],
  y=st_coordinates(paczkomaty_waw)[, 2],
  window=W,
  marks=as.factor(paczkomaty_waw$marka))

paczkomaty.ppp.m<-rjitter(paczkomaty.ppp.m, 0.5)
paczkomaty.ppp.m<-rescale.ppp(paczkomaty.ppp.m, rsc, "km")

summary(paczkomaty.ppp.m) #Wyświetlamy podsumowanie wzorca
## Marked planar point pattern:  3148 points
## Average intensity 6.092221 points per square km
## 
## Coordinates are given to 13 decimal places
## 
## Multitype:
##               frequency proportion  intensity
## DHL                 417 0.13246510 0.80700650
## DPD                 848 0.26937740 1.64110700
## INPOST             1413 0.44885640 2.73453300
## ORLEN               419 0.13310040 0.81087700
## POCZTA POLSKA        51 0.01620076 0.09869863
## 
## Window: polygonal boundary
## single connected closed polygon with 3389 vertices
## enclosing rectangle: [626.506, 655.2603] x [472.2295, 502.1724] km
##                      (28.75 x 29.94 km)
## Window area = 516.724 square km
## Unit of length: 1 km
## Fraction of frame area: 0.6

Funkcje K, L, F, G oraz J

Rozróżniamy 6 typów funkcji:

  • Funkcja K Ripleya. Mierzy, ile punktów wpada do koła o promieniu r wokół losowego punktu, zwiększając stopniowo r. Porównuje wynik z rozkładem czysto losowym (CSR, Poisson). K(r) powyżej linii Poissona → skupienie; poniżej → regularność.
  • Funkcja L. Liniowa (zestandaryzowana) postać funkcji K — sprowadza linię Poissona do prostej L(r) = r. Łatwiejsza do interpretacji wizualnej. L(r) powyżej linii Poissona → skupienie; poniżej → regularność.
  • Funkcja PCF (g). Pochodna funkcji K — mierzy intensywność punktów w pierścieniu o promieniu r (nie kole), więc nie kumuluje efektów dla małych r. g(r) > 1 → skupienie; g(r) < 1 → regularność; g(r) = 1 → losowość.
  • Funkcja G. Mierzy odległość od punktu do jego najbliższego sąsiada (NND). G(r) powyżej linii Poissona → punkty są bliżej siebie niż losowo → skupienie; poniżej → regularność.
  • Funkcja F. Mierzy odległość od losowego miejsca w przestrzeni do najbliższego punktu (empty space distance). F(r) poniżej linii Poissona → puste miejsca są dalej od jakiegokolwiek punktu niż losowo → skupienie; powyżej → regularność.
  • Funkcja J. Wyrażona jako J(r) = (1−G(r))/(1−F(r)). Porównuje siłę skupienia G z „pustką” F. J(r) < 1 → skupienie (G dominuje nad F); J(r) = 1 → losowość (CSR); J(r) > 1 → regularność.

Funkcja K Ripleya mierzy ile punktów wpada do koła o promieniu r wokół losowo wybranego punktu zwiększając stopniowo promień r. Następnie porównujemy z rozkładem losowym (CSR, Poisson)

  • K(r) powyżej linii Poisson sugeruje skupienie
  • K(r) poniżej linii Poisson sugeruje regularność
par(family="TNR", bty="l")
K<-Kest(paczkomaty.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(K, main="Funkcja K Ripleya",
     xlab="r (km)", ylab="K(r)", legendpos="bottomright")

par(family="TNR", bty="l")
L<-Lest(paczkomaty.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(L, main="Funkcja L",
     xlab="r (km)", ylab="L(r)", legendpos="bottomright")

par(family="TNR", bty="l")
F<-Fest(paczkomaty.ppp.um, correction="border", rmax=3)
plot(F, main="Funkcja F (empty space)",
     xlab="r (km)", ylab="F(r)", legendpos="bottomright")

par(family="TNR", bty="l")
G<-Gest(paczkomaty.ppp.um, correction="border", rmax=3)
plot(G, main="Funkcja G (nearest neighbour)",
     xlab="r (km)", ylab="G(r)", legendpos="bottomright")

par(family="TNR", bty="l")
J<-Jest(paczkomaty.ppp.um, correction="border", rmax=3)
plot(J, main="Funkcja J",
     xlab="r (km)", ylab="J(r)", legendpos="bottomright")

Koperty, test MAD/DCLF oraz Clark-Evans

set.seed(2115) #Ustawawiamy seed żeby za każdym razem był ten sam wynik (symulacje Monte Carlo)
n_pacz<-npoints(paczkomaty.ppp.um) #Ile punktów w pełnym datasecie?
#Koperta punktowa nsim=49: p-value min = 1/(49+1) = 0.02. Standardowe wartości nsim: 19, 49, 99
E.ptwise<-envelope(paczkomaty.ppp.um, Linhom,
                   nsim=19,
                   fix.n=TRUE,
                   correction="border",
                   rmax=3,
                   r=seq(0, 3, by=0.05),
                   verbose=FALSE) #Przedziały ufności dla każdego r osobno
par(family="TNR", bty="l")
plot(E.ptwise, main=paste0("Koperta punktowa - Linhom (n=", n_pacz, ")"),
     xlab="r (km)", ylab="L(r)", legendpos="bottomright")

#Koperta globalna - jeden przedział dla całego zakresu r
E.global<-envelope(paczkomaty.ppp.um, Linhom,
                   nsim=19,
                   rank=1,
                   global=TRUE,
                   correction="border",
                   rmax=3,
                   r=seq(0, 3, by=0.05),
                   verbose=FALSE)
par(family="TNR", bty="l")
plot(E.global, main=paste0("Koperta globalna - Linhom (n=", n_pacz, ")"),
     xlab="r (km)", ylab="L(r)", legendpos="bottomright")

#Test MAD i DCLF --> testy, które potwierdzają wnioski z kopert. Niestety, dla 3100 punktów nie wykonają się, dlatego są jako komentarz
#mad.test(paczkomaty.ppp.um, Linhom, nsim=19, verbose=FALSE)
#dclf.test(paczkomaty.ppp.um, Linhom, nsim=19, verbose=FALSE)

#Clark-Evans - R < 1 -> skupienie, R > 1 -> regularność
clarkevans.test(paczkomaty.ppp.um, alternative="clustered") #Hipoteza skupienia
## 
##  Clark-Evans test
##  CDF correction
##  Z-test
## 
## data:  paczkomaty.ppp.um
## R = 0.48321, p-value < 2.2e-16
## alternative hypothesis: clustered (R < 1)
clarkevans.test(paczkomaty.ppp.um, alternative="regular") #Hipoteza regularności
## 
##  Clark-Evans test
##  CDF correction
##  Z-test
## 
## data:  paczkomaty.ppp.um
## R = 0.48321, p-value = 1
## alternative hypothesis: regular (R > 1)

Intensywność oraz relatywne ryzyko

#Dzielimy mapę Warszawy na GRID, czyli siatkę 6x6
qdr<-quadratcount(paczkomaty.ppp.um, nx=6, ny=6)
par(family="TNR", bty="l")
plot(qdr, main="Liczba paczkomatów w kwadratach")

#Intensywność w każdym polu (punkty/km²)
par(family="TNR", bty="l")
plot(intensity(qdr, image=TRUE), main="Intensywność paczkomatów (punkty/km²)")

#KDE - ciągła mapa intensywności
bw_diggle<-bw.diggle(paczkomaty.ppp.um)
kde<-density.ppp(paczkomaty.ppp.um, sigma=bw_diggle)
par(family="TNR", bty="l")
plot(kde, main=paste0("KDE — gęstość paczkomatów (sigma=", round(bw_diggle,3), " km)"))

#Test chi kwadrat
#H0: CSR — punkty losowo rozrzucone
#H1: punkty nie są losowe
qdr.test<-quadrat.test(paczkomaty.ppp.um, nx=6, ny=6)
qdr.test
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  paczkomaty.ppp.um
## X2 = 1168.6, df = 30, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 31 tiles (irregular windows)
#H1: Clustered
qdr.test.clust<-quadrat.test(paczkomaty.ppp.um, nx=6, ny=6, alternative="clustered")
qdr.test.clust
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  paczkomaty.ppp.um
## X2 = 1168.6, df = 30, p-value < 2.2e-16
## alternative hypothesis: clustered
## 
## Quadrats: 31 tiles (irregular windows)
#H1: Regular
qdr.test.reg<-quadrat.test(paczkomaty.ppp.um, nx=6, ny=6, alternative="regular")
qdr.test.reg
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  paczkomaty.ppp.um
## X2 = 1168.6, df = 30, p-value = 1
## alternative hypothesis: regular
## 
## Quadrats: 31 tiles (irregular windows)
par(family="TNR", bty="l")
plot(qdr.test, main="Test chi-kwadrat (quadrat test)")

#Relatywne ryzyko (intensywność danej marki w danym miejscu/intensywność wszystkich marek w danym miejscu)
rr<-relrisk(paczkomaty.ppp.m, sigma=bw.relrisk(paczkomaty.ppp.m))
par(family="TNR", bty="l")
plot(rr, main="Relatywne ryzyko")