Analiza lokalizacji automatów paczkowych (APM) w Warszawie przy wykorzystaniu Analizy wzorca punktowego oraz liniowego
Artur Stankiewicz
04.2026 - Obecnie
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 :)
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). Wszechobecna sieć automatów paczkowych rodzi kluczowe pytania:
- 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?
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 automatów paczkowych (APM).
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…
#Hm... zazwyczaj tytuł pojawia się bezpośrednio nad wykresem... do dalszego dopracowania :). Raportuję użycie AI, nie znałem tej funkcji z pakietu spatstat
library(spatstat)
set.seed(2115)
W<-owin(c(0,2), c(0,2))
par(mfrow=c(1,3), mar=c(1, 1, 1, 1))
plot(rpoispp(lambda=100, win=W), main="Losowy (CSR)", pch=16, cex=1, legend=FALSE)
plot(rSSI(r=0.08, n=100, win=W), main="Regularny", pch=16, cex=1, legend=FALSE)
plot(rMatClust(kappa=10, r=0.1, mu=10, win=W), main="Skupiony", pch=16, cex=1, legend=FALSE)
Kod
Kod przygotowany na podstawie kursu Analiza wzorca punktowego i liniowego w R, prowadzonego przez dr Katerynę Zabarinę, na Wydziale Nauk Ekonomicznych Uniwersytetu Warszawskiego.
Marked oraz unmarked PPP
requiredPackages<-c("sf", "spatstat", "rgugik", "automap",
"magrittr", "sp", "rnaturalearth",
"ggplot2", "tidyverse", "viridis", "osrm")
for(i in requiredPackages)
{if(!require(i, character.only=TRUE)) install.packages(i)}
for(i in requiredPackages)
{library(i, character.only=TRUE)}
Sys.setenv(LANG="en")
windowsFonts(TNR=windowsFont("Times New Roman"))Zacznijmy od wczytania mapy oraz zbioru danych z punktami APM
#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 punktami APM
APM<-read.csv("Lista_APM.csv", fileEncoding="UTF-8-BOM") #Wczytujemy dane z pliku CSV
APM_sf<-st_as_sf(APM, coords=c("lon", "lat"), crs=4326) #Kodujemy koordynaty punktów APM jako długość i szerokość geograficzną
APM_sf<-st_transform(APM_sf, crs=2180) #Przechodzimy z CRS 4326 na 2180
st_crs(APM_sf) #Sprawdzamy czy jest CRS 2180
APM_waw<-st_intersection(APM_sf, warszawa_granica) #Przycinamy dataset do Warszawy## Warning: attribute variables are assumed to be spatially constant throughout
## all geometriesSprawdźmy ile mamy automatów paczkowych w zbiorze danych z podziałem na marki
##
## DHL DPD INPOST ORLEN POCZTA POLSKA
## 417 848 1413 419 51p_mapa<-ggplot() +
geom_sf(data=warszawa_granica, fill="grey95", color="black", linewidth=0.6) +
geom_sf(data=APM_waw, aes(color=marka), size=1.5, alpha=0.7) +
scale_color_viridis_d(option="turbo", name="Marka") +
labs(title="Punkty APM w Warszawie", subtitle=paste("Łącznie", nrow(APM_waw), "punktów")
) +
theme_minimal() +
theme(plot.title=element_text(size=14, face="bold"))
print(p_mapa)
Wyświetlmy podsumowanie wzorca unmarked, czyli takiego, który nie różnicuje punktów ze względu na markę
#Unmarked - wszystkie punkty APM bez rozróżnienia marki
APM.ppp.um<-ppp(
x=st_coordinates(APM_waw)[, 1], #Współrzędne X
y=st_coordinates(APM_waw)[, 2], #Współrzędne Y
window=W) #Punkty w Warszawie
APM.ppp.um<-rjitter(APM.ppp.um, 0.05, retry=TRUE) #Przesuwamy lekko punkty żeby nie było duplikatów (o 0.05 metra).
APM.ppp.um<-rescale.ppp(APM.ppp.um, rsc, "km") #Zmieniamy z metrów na kilometry
summary(APM.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.6Wyświetlmy też podsumowanie wzorca marked, czyli takiego, który różnicuje punkty ze względu na markę
#Marka jako cecha punktów APM
APM.ppp.m<-ppp(
x=st_coordinates(APM_waw)[, 1],
y=st_coordinates(APM_waw)[, 2],
window=W,
marks=as.factor(APM_waw$marka))
APM.ppp.m<-rjitter(APM.ppp.m, 0.5)
APM.ppp.m<-rescale.ppp(APM.ppp.m, rsc, "km")
summary(APM.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.6Narysujmy mapę z podziałem na markę
split_ppp_m<-split(APM.ppp.m) #Dzielimy wzorzec marked na osobne ppp dla każdej marki
par(mfrow=c(2,3), family="TNR", cex.main=1.4, cex.sub=1.1)
plot(split_ppp_m[[1]], main="DHL", pch=16, cex=0.5, cols="red")
mtext(paste(npoints(split_ppp_m[[1]]), "APM"), side=1, line=2.5, cex=1.0)
plot(split_ppp_m[[2]], main="DPD", pch=16, cex=0.5, cols="green")
mtext(paste(npoints(split_ppp_m[[2]]), "APM"), side=1, line=2.5, cex=1.0)
plot(split_ppp_m[[3]], main="InPost", pch=16, cex=0.5, cols="blue")
mtext(paste(npoints(split_ppp_m[[3]]), "APM"), side=1, line=2.5, cex=1.0)
plot(split_ppp_m[[4]], main="Orlen", pch=16, cex=0.5, cols="orange")
mtext(paste(npoints(split_ppp_m[[4]]), "APM"), side=1, line=2.5, cex=1.0)
plot(split_ppp_m[[5]], main="Poczta Polska", pch=16, cex=0.5, cols="black")
mtext(paste(npoints(split_ppp_m[[5]]), "APM"), side=1, line=2.5, cex=1.0)
plot.new() #Puste pole na legendę
legend("center", legend=c("DHL","DPD","InPost","Orlen","Poczta Polska"),
col=c("red","green","blue","orange","black"),
pch=16, cex=1.3, title="Marka", bty="n", pt.cex=1.5)
Funkcje K, L, F, G oraz J. Koperty oraz testy MAD/DCLF/Clark-Evans
Zanim przejdziemy do funkcji… przypomnijmy sobie czym się różni funkcja homogeniczna od niehomogenicznej. Najprościej rzecz ujmując:
- funkcja homogeniczna zakłada, że intensywność punktów w przestrzeni jest taka sama.
- Odwrotnie jest w przypadku funkcji niehomogenicznej, która bierze pod uwagę intensywność.
A jak należy interpretować te funkcje? Charakterystyczne są:
- pokrycia z funkcją Poissona, co może sugerować, że mamy do czynienia z CSR
- funkcje nad lub pod funkcją Poissona, co sugeruje, że mamy do czynienia ze wzorcem (skupiony, a może regularny?)
- peaków, przecięć, pokryć z funkcją Poissona od określonej długości promienia
Zwróćcie proszę uwagę na różne osie x oraz y na wykresach
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ść
Wniosek: funkcja K Ripleya sugeruje, że mamy do czynienia ze skupieniem.
par(mfrow=c(1,2), mar=c(3,3,2,1), family="TNR", bty="l")
K<-Kest(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(K, main="Homogeniczna",
xlab="r (km)", ylab="K(r)", legendpos="bottomright")
Kin<-Kinhom(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(Kin, main="Niehomogeniczna",
xlab="r (km)", ylab="K(r)", legendpos="bottomright")
Funkcja L to liniowa standaryzacja funkcji K.
- L(r) powyżej linii Poisson sugeruje skupienie
- L(r) poniżej linii Poisson sugeruje regularność
Wniosek: analogiczny do funkcji K Ripleya :)
par(mfrow=c(1,2), mar=c(3,3,2,1), family="TNR", bty="l")
L<-Lest(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(L, main="Homogeniczna",
xlab="r (km)", ylab="K(r)", legendpos="bottomright")
Lin<-Linhom(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(Lin, main="Niehomogeniczna",
xlab="r (km)", ylab="K(r)", legendpos="bottomright")
Funkcja F (empty space) to funkcja, która mierzy odległość od pustego w przestrzeni punktu do najbliższego punktu APM.
- F(r) powyżej linii Poisson sugeruje regularność
- F(r) poniżej linii Poisson sugeruje skupienie
Wniosek: sugeruje, że jest więcej “pustego miejsca” niż w przypadku rozkładu Poissona
par(mfrow=c(1,2), mar=c(3,3,2,1), family="TNR", bty="l")
F<-Fest(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(F, main="Homogeniczna",
xlab="r (km)", ylab="K(r)", legendpos="bottomright")
Fin<-Finhom(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(Fin, main="Niehomogeniczna",
xlab="r (km)", ylab="K(r)", legendpos="bottomright")
Funkcja G (nearest neighbour) to funkcja, która mierzy odległość od punktu APM w przestrzeni do najbliższego punktu APM.
- G(r) powyżej linii Poisson sugeruje skupienie
- G(r) poniżej linii Poisson sugeruje regularność
Wniosek: sugeruje, że punkty stoją bliżej siebie niż gdyby to było w przypadku rozkładu Poissona
par(mfrow=c(1,2), mar=c(3,3,2,1), family="TNR", bty="l")
G<-Gest(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(G, main="Homogeniczna",
xlab="r (km)", ylab="K(r)", legendpos="bottomright")
Gin<-Ginhom(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(Gin, main="Niehomogeniczna",
xlab="r (km)", ylab="K(r)", legendpos="bottomright")
Funkcja J służy do porównania czy mamy bardziej do czynienia z pustą przestrzenią (empty space, funkcja F) czy najbliższego sąsiada (nearest neighbour, funkcja G)
- J(r) powyżej 1 wskazuje, że silniejszy jest efekt funkcji F
- J(r) poniżej 1 wskazuje, że silniejszy jest efekt funkcji G
Wniosek: zdecydowanie silniejszy efekt funkcji G względem funkcji F
par(mfrow=c(1,2), mar=c(3,3,2,1), family="TNR", bty="l")
J<-Jest(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(J, main="Homogeniczna",
xlab="r (km)", ylab="K(r)", legendpos="bottomright")
Jin<-Jinhom(APM.ppp.um, correction="border", rmax=3, r=seq(0, 3, by=0.05))
plot(Jin, main="Niehomogeniczna",
xlab="r (km)", ylab="K(r)", legendpos="bottomright")
Narysowaliśmy funkcje i wstępnie postawiliśmy wnioski. Dla potwierdzenia wniosków posłużmy się innymi testami. Zacznijmy od kopert.
Koperta punktowa
set.seed(2115) #Ustawawiamy seed żeby za każdym razem był ten sam wynik
n_pacz<-npoints(APM.ppp.um) #Ile punktów w pełnym datasecie?
E.ptwise<-envelope(APM.ppp.um, Linhom, #Koperta dla niehomogenicznej funkcji L
nsim=19, #Liczba symulacj, możemy wybrać dowolną wielkość. Służy do obliczenia p-value zgodnie ze wzorem 2/(nsim+1). Porównujemy z przyjętym poziomem istotności, np. alfa = 5%. Standardowe wartości nsim to 19, 49, 99
fix.n=TRUE,
correction="border", #Dla powyżej 3 tysięcy punktów correction z automatu jest ustawione na border
rmax=3,
r=seq(0, 3, by=0.05),
verbose=FALSE)
par(family="TNR", bty="l")
plot(E.ptwise, main=("Koperta punktowa"),
xlab="r (km)", ylab="L(r)", legendpos="bottomright")
Koperta globalna
set.seed(2115)
#Uwaga, tutaj p-value liczone jest jako 1/(nsim+1)
E.global<-envelope(APM.ppp.um, Linhom,
nsim=19,
rank=1,
global=TRUE,
correction="border",
fix.n=T,
rmax=3,
r=seq(0, 3, by=0.05),
verbose=FALSE)
par(family="TNR", bty="l")
plot(E.global, main=paste0("Koperta globalna"),
xlab="r (km)", ylab="L(r)", legendpos="bottomright")
Możemy też przeprowadzić test MAD i DCLF, które potwierdzają wnioski z kopert. Niestety, przy tak dużej liczbie punktów test się nie wykona
mad.test(APM.ppp.um, Linhom, nsim=19, verbose=FALSE)
dclf.test(APM.ppp.um, Linhom, nsim=19, verbose=FALSE)Teraz testy Clark-Evans. Postanowiłem sprawdzić dwie hipotezy alternatywne:
- Hipoteza zerowa H0: punkty są zgodne z czysto losowym wzorcem (CSR)
- H1(a): punkty są zgodne ze skupionym wzorcem
- H1(b): punkty są zgodne z regularnym wzorcem
Na podstawie poniższego testu istnieją podstawy do odrzucenia hipotezy zerowej o CSR na rzecz hipotezy alternatywnej H1(a) o skupieniu (p-value=2,2e-16)
##
## Clark-Evans test
## CDF correction
## Z-test
##
## data: APM.ppp.um
## R = 0.48321, p-value < 2.2e-16
## alternative hypothesis: clustered (R < 1)Na podstawie poniższego testu nie ma podstaw do odrzucenia hipotezy zerowej o CSR na rzecz hipotezy alternatywnej H1(b) o regularności (p-value=1)
##
## Clark-Evans test
## CDF correction
## Z-test
##
## data: APM.ppp.um
## R = 0.48321, p-value = 1
## alternative hypothesis: regular (R > 1)Wniosek: testy wskazują jednoznacznie: punkty APM nie są rozmieszczone losowo, a zgodnie ze skupionym wzorcem.
Intensywność oraz relatywne ryzyko
Intensywność to najprościej rzecz ujmując liczba punktów przypadająca na daną jednostkę przestrzenną (na przykład kilometr kwadratowy). Intensywność możemy przedstawić na dwa sposoby:
1. podzielić mapę Warszawy na tak zwany “grid”, czyli siatkę i przypisać każdej siatce punkty APM.
#Dzielimy mapę Warszawy na GRID, czyli siatkę 9x9
qdr<-quadratcount(APM.ppp.um, nx=9, ny=9)
par(family="TNR", bty="l")
plot(qdr, main="Liczba punktów APM w kwadratach")
a następnie stworzyć mapę ciepła (punkty/km^2).
#Intensywność w każdym polu (punkty/km²)
par(family="TNR", bty="l")
plot(intensity(qdr, image=TRUE), main="Intensywność punktów APM (punkty/km^2)")
2. narysować mapę ciepła używając do tego KDE (Kernel Density Estimation). Zacznijmy od wyboru sigmy. W praktyce, wartość sigmy oznacza rozmycie punktu w przestrzeni. Można wybrać dowolną wartość. Narysujmy 6 przykładowych sigm i zobaczmy jakie ma to przełożenie na mapę intensywności. Przy:
- Małej sigmie rozmycie jest małe, widać lokalne punkty
- Dużej sigmie rozmycie jest duże, widać bardziej globalny trend
Co wybrać? To zależy na czym chcemy się dokładnie skupić: globalne trendy czy lokalna sytuacja?
Oczywiście, istnieje również funkcja, która wybierze najlepszą sigmę… tylko problem pojawia się, gdy chcemy różnicować punkty. Poniżej znajduje się optymalny poziom sigmy dla poszczególnych firm:
#Raportuję tutaj użycie AI, chciałem po prostu zobaczyć jak bw.diggle zmieni się przy wyborze marki. Coś mi nie pasowała początkowa mapa z relrisk.
sigma_marka<-sapply(split_ppp_m, bw.diggle)
sigma_marka## DHL DPD INPOST ORLEN POCZTA POLSKA
## 0.7104178 0.2532182 0.4149965 0.7104178 1.7162568Dlatego dla ujednolicenia map, co pozwoli nam na porównanie, będę posługiwał się jedną sigmą.
Zwróćcie proszę uwagę na różne przedziały legendy
sigmy<-c(0.1, 0.5, 0.7, 1, 1.5, 2)
par(mfrow=c(2,3), mar=c(1,1,2,1), family="TNR")
for(s in sigmy)
{plot(density(APM.ppp.um, sigma=s),
main=paste0("Sigma=", s,"km"),
ribbon=TRUE)}
Postanowiłem wybrać sigmę=1. Narysujmy teraz mapy intensywności dla punktów APM łącznie oraz z podziałem na marki.
Zwróćcie proszę uwagę na różne przedziały legendy
sigma<-1
par(mfrow=c(2,3), mar=c(1,1,2,1), family="TNR")
plot(density(APM.ppp.um, sigma=sigma), main="Łącznie", ribbon=TRUE)
for(i in 1:5) {
plot(density(split_ppp_m[[i]], sigma=sigma),
main=names(split_ppp_m)[i],
ribbon=TRUE)
}
I teraz najciekawsze: relatywne ryzyko.
## [1] "DHL" "DPD" "INPOST" "ORLEN"
## [5] "POCZTA POLSKA"A teraz narysujmy relatywne ryzyko marki względem wszystkich punktów APM
#Relatywne ryzyko (intensywność danej marki w danym miejscu/intensywność wszystkich marek w danym miejscu)
rr<-relrisk(APM.ppp.m, sigma=1)
par(family="TNR", bty="l")
plot(rr, main="Relatywne ryzyko")
Dystans z punktu widzenia konsumenta
Ale o tym w przyszłym poście :)