In pbiecek/SPAG: Indexes of spatial agglomeration

---

title: "Teoretyczne próbkowanie przestrzenne" date: "r Sys.Date()" output: rmarkdown::html_vignette: toc: true vignette: > %\VignetteIndexEntry{SPAG tutorial} %\VignetteEngine{knitr::rmarkdown} %\VignetteDepends{Cairo} %\VignetteEncoding{UTF-8} \usepackage[utf8]{inputenc}

---

Istniejący algorytm próbkowania jednorodnego punktów

knitr::opts_chunk$set(echo = TRUE)

Indeks Distance zdefiniowany jest jako stosunek średniej odległości pomiędzy rzeczywistymi firmami do średniej odległości pomiędzy firmami w założeniu, że byłyby one rozłożone w sposób równomierny na całym obszarze: $$I_{dist} = \frac{ \frac{1}{k} \sum_{i,j} d(i,j) }{\frac{1}{ \hat{k} } \sum_{i,j} \hat{d}(i,j) } $$ Przy takim podejściu niezbędne jest określenie tego, co oznacza "rozkład równomierny na obszarze". W poniższym tekście postaram się przedstawić obecnie wykorzystywaną w tym celu metodę.

Dane do testów.

Na potrzeby dalszej analizy zostanie wykorzystana następująca mapa Polski:

library(SPAG)
n <- 1374
newCoordinateSystem<-"+proj=longlat +datum=WGS84"
region<-spTransform(ShapefilePoland, CRS(newCoordinateSystem))

mapaPolski <- ggplot() +
  geom_polygon(data=region, aes(long, lat, group=group), colour='#808080', fill=NA) +
  theme_nothing() +
  labs(long="longitude", lat="latitude")
mapaPolski

Obecne podejście

Równomierne rozmieszczenie punktów na danym obszarze otrzymuję się poprzez naniesienie kraty punktów o odpowiednim rozmiarze na dany obszar, a następnie zostawienie części wspólnej zbioru punktów i obszaru.

Proces ten można podzielić na trzy etapy. Pierwszym z nich jest wyznaczenie rozmiaru kraty, która posłuży do wyznaczania punktów. W tym celu wyznacza się stosunek powierzchni ramki ograniczającej rozważany obszar do powierzchni samego obszaru. Otrzymany stosunek jest następnie wykorzystywany do wyznaczenia rzeczywistej liczby punktów potrzebnych do naniesienia na obszar:

$$ \hat{n} = \frac{P_{box} }{ P_{area} } \cdot n$$

x <- region
bb = bbox(region)

getArea = function(x) {
    getAreaPolygons = function(x) {
      holes = unlist(lapply(x@Polygons, function(x) x@hole))
      areas = unlist(lapply(x@Polygons, function(x) x@area))
      area = ifelse(holes, -1, 1) * areas
      area
    }
    sum(unlist(lapply(region@polygons, getAreaPolygons)))
}

  area = getArea(x)
  res <- NULL
  bb.area = (bb[1,2]-bb[1,1]) * (bb[2,2]-bb[2,1]) 
  print( paste("Obszar ramki wynosi:", bb.area))
  n_tot = round(n * bb.area/area) 
  print( paste("Rzeczywista liczba nanoszonych punktów wynosi:",n_tot))

W drugim etapie wyznaczana jest siatka punktów, która jest następnie nanoszona na obszar. W tym celu niezbędne jest wyliczenie długości boku wycinka kraty:

$$ c = \sqrt{ \frac{P_{box} }{\hat{n}} } = \sqrt{ P_{box} \cdot\frac{P_{area} }{ P_{box} \cdot n} }=\sqrt{ \frac{P_{area} }{n } }$$

znając już rozmiar wycinka siatki niezbędne jest umiejscowienie jej na dany obszar. Domyślnie, rysowanie siatki rozpoczyna się od umieszczenia jej w lewym dolnym rogu, a następnie przesunięciu jej o losowy wektor o dodatnich współrzędnych.

  offset = runif(2)
  pw =0.5 
  nsig <- 20
  cellsize = (prod(apply(bb, 1, diff))/n_tot) ^ pw 
  cellsize = min(cellsize, min(apply(bb, 1, diff)))
  cellsize = rep(cellsize, nrow(bb))
  min.coords = pmax(bb[,1], signif(bb[,1] + offset * cellsize, nsig)) # nsig =20
  expand.grid.arglist = list()

  for (i in 1:nrow(bb)) {
    name = paste("x", i, sep = "")
    sign = ifelse(min.coords[i] < bb[i,2], 1, -1) # sprawdzenie jeszcze raz kierunku
    expand.grid.arglist[[name]] = seq(min.coords[i], bb[i,2], sign * cellsize[i])
  }

  xy = do.call(expand.grid, expand.grid.arglist)
  attr(xy, "cellsize") = cellsize
  mapaPolski + geom_point(data=xy, aes(x1,x2))

W ostatnim etapie odcinane są punkty, które nie należą do rozważanej powierzchnii:

 sel = xy[,1] >= bb[1,1] & xy[,1] <= bb[1,2] & 
    xy[,2] >= bb[2,1] & xy[,2] <= bb[2,2]
  xy = xy[sel,, drop = FALSE]
  rownames(xy) = NULL
  pts <- SpatialPoints(xy, CRS(proj4string(x)))
  Over_pts_x <- over(pts, geometry(x))

  Not_NAs <- !is.na(Over_pts_x)
  res <- pts[Not_NAs]
  res <- as.data.frame(res)
  mapaPolski + geom_point(data=res, aes(x1,x2))

Deterministyczność

Fakt, że ostateczne umieszczenie siatki jest dokonywane w sposób losowy implikuje zmienną liczbę punktów, które zostaną naniesione na rozważaną powierzchnię. To z kolei powoduje, że wyliczana wartość indeksu Distance również nie będzie deterministyczna. W celu usunięcia tej zmienności w kolejnych implementacjach funkcji SPAG przesunięcie losowe ustawiane jest na stałą wartość 0.

pbiecek/SPAG documentation built on May 24, 2019, 10:36 p.m.