Skip to main content
added required changes to spatstat to support examples
Source Link
Jeffrey Evans
  • 32k
  • 2
  • 48
  • 97

This is a difficult question as there just have not been many, if any, spatial process statistics developed for line features. Without seriously digging into equations and code, point process statistics are not readily applicable to linear features and thus, statistically invalid. This is because the null, that a given pattern is tested against, is based on point events given a random field and not linear dependencies in the random field. I have to say that I do not even know what the null would be as far as intensity and arrangement/orientation would be even more difficult.

library(spatstat)
library(spatstat.core)
library(raster)
library(spatialEco)

data(copper)
l <- copper$Lines
l <- rotate.psp(l, pi/2)
library(spatstat)
library(spatstat.core)
library(spatstat.geom)
  data(copper)
  l <- copper$Lines
  l <- rotate.psp(l, pi/2)

Lr <- function (...) {
 K <- Kest(...)
  nama <- colnames(K)
   K <- K[, !(nama %in% c("rip", "ls"))]
   L <- eval.fv(sqrt(K/pi)-bw)
  L <- rebadge.fv(L, substitute(L(r), NULL), "L")
 return(L)
}

### Ripley's K ( Besag L(r) ) for start locations
start <- endpoints.psp(l, which="first")
marks(start) <- factor("start")
W <- start$window
area <- area.owin(W)
lambda <- start$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
       bw=seq(0, rmax, by=1)
( Lenv <- plot( envelope(start, fun="Lr", r=seq(0, rmax, by=1), nsim=199, nrank=5) ) )

### Ripley's K ( Besag L(r) ) for end locations
stop <- endpoints.psp(l, which="second")
  marks(stop) <- factor("stop")
W <- stop$window
area <- area.owin(W)
lambda <- stop$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
       bw=seq(0, rmax, by=1)
( Lenv <- plot( envelope(start, fun="Lr", r=seq(0, rmax, by=1), nsim=199, nrank=5) ) )

### Ripley's Cross-K ( Besag L(r) ) for start/stop
sdata.ppp <- superimpose(start, stop)
( Lenv <- plot(envelope(sdata.ppp, fun="Kcross", r=bw, i="start", j="stop", nsim=199,nrank=5, 
                 transform=expression(sqrt(./pi)-bw), global=TRUE) ) )
             

This is a difficult question as there just have not been many, if any, spatial process statistics developed for line features. Without seriously digging into equations and code, point process statistics are not readily applicable to linear features and thus, statistically invalid. This is because the null, that a given pattern is tested against, is based on point events and not linear dependencies in the random field. I have to say that I do not even know what the null would be as far as intensity and arrangement/orientation would be even more difficult.

library(spatstat)
library(raster)
library(spatialEco)

data(copper)
l <- copper$Lines
l <- rotate.psp(l, pi/2)
library(spatstat)
  data(copper)
  l <- copper$Lines
  l <- rotate.psp(l, pi/2)

Lr <- function (...) {
 K <- Kest(...)
  nama <- colnames(K)
   K <- K[, !(nama %in% c("rip", "ls"))]
   L <- eval.fv(sqrt(K/pi)-bw)
  L <- rebadge.fv(L, substitute(L(r), NULL), "L")
 return(L)
}

### Ripley's K ( Besag L(r) ) for start locations
start <- endpoints.psp(l, which="first")
marks(start) <- factor("start")
W <- start$window
area <- area.owin(W)
lambda <- start$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
( Lenv <- plot( envelope(start, fun="Lr", r=seq(0, rmax, by=1), nsim=199, nrank=5) ) )

### Ripley's K ( Besag L(r) ) for end locations
stop <- endpoints.psp(l, which="second")
  marks(stop) <- factor("stop")
W <- stop$window
area <- area.owin(W)
lambda <- stop$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
( Lenv <- plot( envelope(start, fun="Lr", r=seq(0, rmax, by=1), nsim=199, nrank=5) ) )

### Ripley's Cross-K ( Besag L(r) ) for start/stop
sdata.ppp <- superimpose(start, stop)
( Lenv <- plot(envelope(sdata.ppp, fun="Kcross", r=bw, i="start", j="stop", nsim=199,nrank=5, 
                 transform=expression(sqrt(./pi)-bw), global=TRUE) ) )
             

This is a difficult question as there just have not been many, if any, spatial process statistics developed for line features. Without seriously digging into equations and code, point process statistics are not readily applicable to linear features and thus, statistically invalid. This is because the null, that a given pattern is tested against, is based on point events given a random field and not linear dependencies. I have to say that I do not even know what the null would be as far as intensity and arrangement/orientation would be even more difficult.

library(spatstat)
library(spatstat.core)
library(raster)
library(spatialEco)

data(copper)
l <- copper$Lines
l <- rotate.psp(l, pi/2)
library(spatstat)
library(spatstat.core)
library(spatstat.geom)
  data(copper)
  l <- copper$Lines
  l <- rotate.psp(l, pi/2)

Lr <- function (...) {
 K <- Kest(...)
  nama <- colnames(K)
   K <- K[, !(nama %in% c("rip", "ls"))]
   L <- eval.fv(sqrt(K/pi)-bw)
  L <- rebadge.fv(L, substitute(L(r), NULL), "L")
 return(L)
}

### Ripley's K ( Besag L(r) ) for start locations
start <- endpoints.psp(l, which="first")
marks(start) <- factor("start")
W <- start$window
area <- area.owin(W)
lambda <- start$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
       bw=seq(0, rmax, by=1)
( Lenv <- plot( envelope(start, fun="Lr", r=seq(0, rmax, by=1), nsim=199, nrank=5) ) )

### Ripley's K ( Besag L(r) ) for end locations
stop <- endpoints.psp(l, which="second")
  marks(stop) <- factor("stop")
W <- stop$window
area <- area.owin(W)
lambda <- stop$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
       bw=seq(0, rmax, by=1)
( Lenv <- plot( envelope(start, fun="Lr", r=seq(0, rmax, by=1), nsim=199, nrank=5) ) )

### Ripley's Cross-K ( Besag L(r) ) for start/stop
sdata.ppp <- superimpose(start, stop)
( Lenv <- plot(envelope(sdata.ppp, fun="Kcross", r=bw, i="start", j="stop", nsim=199,nrank=5, 
                 transform=expression(sqrt(./pi)-bw), global=TRUE) ) )
             
deleted 2 characters in body
Source Link
Jeffrey Evans
  • 32k
  • 2
  • 48
  • 97
library(spatstat)
  data(copper)
  l <- copper$Lines
  l <- rotate.psp(l, pi/2)

Lr <- function (...) {
 K <- Kest(...)
  nama <- colnames(K)
   K <- K[, !(nama %in% c("rip", "ls"))]
   L <- eval.fv(sqrt(K/pi)-bw)
  L <- rebadge.fv(L, substitute(L(r), NULL), "L")
 return(L)
}

### Ripley's K ( Besag L(r) ) for start locations
start <- endpoints.psp(l, which="first")
marks(start) <- factor("start")
W <- start$window
area <- area.owin(W)
lambda <- start$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
( Lenv <- plot( envelope(start, fun="Lr", r=seq(0, rmax, by=1), nsim=199, nrank=5) ) )

### Ripley's K ( Besag L(r) ) for startend locations
stop <- endpoints.psp(l, which="second")
  marks(stop) <- factor("stop")
W <- stop$window
area <- area.owin(W)
lambda <- stop$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
( Lenv <- plot( envelope(start, fun="Lr", r=seq(0, rmax, by=1), nsim=199, nrank=5) ) )

### Ripley's Cross-K ( Besag L(r) ) for start/stop
sdata.ppp <- superimpose(start, stop)
( Lenv <- plot(envelope(sdata.ppp, fun="Kcross", r=bw, i="start", j="stop", nsim=199,nrank=5, 
                 transform=expression(sqrt(./pi)-bw), global=TRUE) ) )
             
library(spatstat)
  data(copper)
  l <- copper$Lines
  l <- rotate.psp(l, pi/2)

Lr <- function (...) {
 K <- Kest(...)
  nama <- colnames(K)
   K <- K[, !(nama %in% c("rip", "ls"))]
   L <- eval.fv(sqrt(K/pi)-bw)
  L <- rebadge.fv(L, substitute(L(r), NULL), "L")
 return(L)
}

### Ripley's K ( Besag L(r) ) for start locations
start <- endpoints.psp(l, which="first")
marks(start) <- factor("start")
W <- start$window
area <- area.owin(W)
lambda <- start$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
( Lenv <- plot( envelope(start, fun="Lr", r=seq(0, rmax, by=1), nsim=199, nrank=5) ) )

### Ripley's K ( Besag L(r) ) for start locations
stop <- endpoints.psp(l, which="second")
  marks(stop) <- factor("stop")
W <- stop$window
area <- area.owin(W)
lambda <- stop$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
( Lenv <- plot( envelope(start, fun="Lr", r=seq(0, rmax, by=1), nsim=199, nrank=5) ) )

### Ripley's Cross-K ( Besag L(r) ) for start/stop
sdata.ppp <- superimpose(start, stop)
( Lenv <- plot(envelope(sdata.ppp, fun="Kcross", r=bw, i="start", j="stop", nsim=199,nrank=5, 
                 transform=expression(sqrt(./pi)-bw), global=TRUE) ) )
             
library(spatstat)
  data(copper)
  l <- copper$Lines
  l <- rotate.psp(l, pi/2)

Lr <- function (...) {
 K <- Kest(...)
  nama <- colnames(K)
   K <- K[, !(nama %in% c("rip", "ls"))]
   L <- eval.fv(sqrt(K/pi)-bw)
  L <- rebadge.fv(L, substitute(L(r), NULL), "L")
 return(L)
}

### Ripley's K ( Besag L(r) ) for start locations
start <- endpoints.psp(l, which="first")
marks(start) <- factor("start")
W <- start$window
area <- area.owin(W)
lambda <- start$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
( Lenv <- plot( envelope(start, fun="Lr", r=seq(0, rmax, by=1), nsim=199, nrank=5) ) )

### Ripley's K ( Besag L(r) ) for end locations
stop <- endpoints.psp(l, which="second")
  marks(stop) <- factor("stop")
W <- stop$window
area <- area.owin(W)
lambda <- stop$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
( Lenv <- plot( envelope(start, fun="Lr", r=seq(0, rmax, by=1), nsim=199, nrank=5) ) )

### Ripley's Cross-K ( Besag L(r) ) for start/stop
sdata.ppp <- superimpose(start, stop)
( Lenv <- plot(envelope(sdata.ppp, fun="Kcross", r=bw, i="start", j="stop", nsim=199,nrank=5, 
                 transform=expression(sqrt(./pi)-bw), global=TRUE) ) )
             
Bounty Ended with 200 reputation awarded by underdark
deleted 1109 characters in body
Source Link
Jeffrey Evans
  • 32k
  • 2
  • 48
  • 97

One also wonders "what if" you just performed a point pattern analysis using a univariate or cross analysis statistic on the start and stop points, invariant of the lines. In a univariate analysis you would compare the results of the start and stop points to see if there consistency in clustering between the two point patterns. This could be done via a f-hat, G-hat or KRipley's-K-hat (for unmarked point processes). Another approach would be a Cross analysis (eg., cross K or Besag L-K) where the two point processes are tested simultaneously by marking them as [1[start,2] start/stopstop]. This would indicate the distance relationships in the clustering process between the start and stop points. However, spatial dependency (nonstaionarity) on an underlying intensity process can be an issue in these types of models making them inhomogeniousinhomogeneous and requiring a different model. Ironically, inhomogeneous process is modeled using an intensity function which, brings us full circle back to density thus, supporting the idea of using a scale-integrated density as a measure of clustering.

library(spatstat)
  data(copper)
  l <- copper$Lines
  l <- rotate.psp(l, pi/2)

Lr <- function (...) {
 K <- Kest(...)
  nama <- colnames(K)
   K <- K[, !(nama %in% c("rip", "ls"))]
   L <- eval.fv(sqrt(K/pi)-bw)
  L <- rebadge.fv(L, substitute(L(r), NULL), "L")
 return(L)
}

### Ripley's K ( Besag L(r) ) for start locations
start <- endpoints.psp(l, which="first")
marks(start) <- factor("start")
 
W <- start$window
area <- area.owin(W)
lambda <- start$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
bw <- seq(0, rmax, by=1) 
 Lenv <- envelope(start, fun="Lr", r=bw, nsim=199, nrank=5)
  plot(Lenv, main="Start locations L(r)", xlab="Distance", ylab="L(r)", legend=F, 
       col=c("white","black","grey","grey"), 
       lty=c(1,2,2,2), lwd=cenvelope(2,1,1start,1))
       polygon(c(Lenv$rfun="Lr", revr=seq(Lenv$r))0, c(Lenv$lormax, rev(Lenv$hi)by=1), col="lightgrey", border="grey")
       lines(supsmu(bwnsim=199, Lenv$obsnrank=5), lwd=2)
       lines(bw, Lenv$theo, lwd=1, lty=2)
 
### Ripley's K ( Besag L(r) ) for stopstart locations
stop <- endpoints.psp(l, which="second")
  marks(stop) <- factor("stop")
W <- stop$window
area <- area.owin(W)
lambda <- stop$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
bw <- seq(0, rmax, by=1) 
Lenv <- envelope(stop, fun="Lr", r=bw, nsim=199, nrank=5)
  plot(Lenv, main="stop locations L(r)", xlab="Distance", ylab="L(r)", legend=F, 
       col=c("white","black","grey","grey"), 
       lty=c(1,2,2,2), lwd=cenvelope(2,1,1start,1))
       polygon(c(Lenv$rfun="Lr", revr=seq(Lenv$r))0, c(Lenv$lormax, rev(Lenv$hi)by=1), col="lightgrey"nsim=199, border="grey"nrank=5)
       lines(supsmu(bw, Lenv$obs), lwd=2)
       lines(bw, Lenv$theo, lwd=1, lty=2)

### Ripley's Cross-K ( Besag L(r) ) for start/stop
sdata.ppp <- superimpose(start, stop)
( Lenv <- plot(envelope(sdata.ppp, fun="Kcross", r=bw, i="start", j="stop", nsim=199,nrank=5, 
                 transform=expression(sqrt(./pi)-bw), global=TRUE) 
                 
plot(Lenv, main="L(d) start/stop", xlab="Distance", ylab="L(r)", legend=F, 
     col=c("white","black","grey","grey"), 
     lty=c(1,2,2,2), lwd=c(2,1,1,1) )
     polygon(c(Lenv$r, rev(Lenv$r)), c(Lenv$lo, rev(Lenv$hi)), col="lightgrey", border="grey")
     lines(supsmu(bw, Lenv$obs), lwd=2)
      lines(bw, Lenv$theo, lwd=1, lty=2)
     legend("topleft", c(expression(hat(L)(r)), "Simulation Envelope"), pch=c(-32,22),
           col=c("black","grey"), lty=c(1,0), lwd=c(2,0), pt.bg=c("white","grey"))     
 

One also wonders "what if" you just performed a point pattern analysis using a univariate or cross analysis statistic on the start and stop points, invariant of the lines. In a univariate analysis you would compare the results of the start and stop points to see if there consistency in clustering between the two point patterns. This could be done via a f-hat, G-hat or K-hat (for unmarked point processes). Another approach would be a Cross analysis (eg., cross K or Besag L) where the two point processes are tested simultaneously by marking them as [1,2] start/stop. This would indicate the distance relationships in the clustering process between the start and stop points. However, spatial dependency (nonstaionarity) on an underlying intensity process can be an issue in these types of models making them inhomogenious and requiring a different model.

library(spatstat)
  data(copper)
  l <- copper$Lines
  l <- rotate.psp(l, pi/2)

Lr <- function (...) {
 K <- Kest(...)
  nama <- colnames(K)
   K <- K[, !(nama %in% c("rip", "ls"))]
   L <- eval.fv(sqrt(K/pi)-bw)
  L <- rebadge.fv(L, substitute(L(r), NULL), "L")
 return(L)
}

### Ripley's K ( Besag L(r) ) for start locations
start <- endpoints.psp(l, which="first")
marks(start) <- factor("start")
 
W <- start$window
area <- area.owin(W)
lambda <- start$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
bw <- seq(0, rmax, by=1) 
 Lenv <- envelope(start, fun="Lr", r=bw, nsim=199, nrank=5)
  plot(Lenv, main="Start locations L(r)", xlab="Distance", ylab="L(r)", legend=F, 
       col=c("white","black","grey","grey"), 
       lty=c(1,2,2,2), lwd=c(2,1,1,1))
       polygon(c(Lenv$r, rev(Lenv$r)), c(Lenv$lo, rev(Lenv$hi)), col="lightgrey", border="grey")
       lines(supsmu(bw, Lenv$obs), lwd=2)
       lines(bw, Lenv$theo, lwd=1, lty=2)
 
### Ripley's K ( Besag L(r) ) for stop locations
stop <- endpoints.psp(l, which="second")
  marks(stop) <- factor("stop")
W <- stop$window
area <- area.owin(W)
lambda <- stop$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
bw <- seq(0, rmax, by=1) 
Lenv <- envelope(stop, fun="Lr", r=bw, nsim=199, nrank=5)
  plot(Lenv, main="stop locations L(r)", xlab="Distance", ylab="L(r)", legend=F, 
       col=c("white","black","grey","grey"), 
       lty=c(1,2,2,2), lwd=c(2,1,1,1))
       polygon(c(Lenv$r, rev(Lenv$r)), c(Lenv$lo, rev(Lenv$hi)), col="lightgrey", border="grey")
       lines(supsmu(bw, Lenv$obs), lwd=2)
       lines(bw, Lenv$theo, lwd=1, lty=2)

### Ripley's Cross-K ( Besag L(r) ) for start/stop
sdata.ppp <- superimpose(start, stop)
Lenv <- envelope(sdata.ppp, fun="Kcross", r=bw, i="start", j="stop", nsim=199,nrank=5, 
                 transform=expression(sqrt(./pi)-bw), global=TRUE) 
                 
plot(Lenv, main="L(d) start/stop", xlab="Distance", ylab="L(r)", legend=F, 
     col=c("white","black","grey","grey"), 
     lty=c(1,2,2,2), lwd=c(2,1,1,1) )
     polygon(c(Lenv$r, rev(Lenv$r)), c(Lenv$lo, rev(Lenv$hi)), col="lightgrey", border="grey")
     lines(supsmu(bw, Lenv$obs), lwd=2)
      lines(bw, Lenv$theo, lwd=1, lty=2)
     legend("topleft", c(expression(hat(L)(r)), "Simulation Envelope"), pch=c(-32,22),
           col=c("black","grey"), lty=c(1,0), lwd=c(2,0), pt.bg=c("white","grey"))     
 

One also wonders "what if" you just performed a point pattern analysis using a univariate or cross analysis statistic on the start and stop points, invariant of the lines. In a univariate analysis you would compare the results of the start and stop points to see if there consistency in clustering between the two point patterns. This could be done via a f-hat, G-hat or Ripley's-K-hat (for unmarked point processes). Another approach would be a Cross analysis (eg., cross-K) where the two point processes are tested simultaneously by marking them as [start,stop]. This would indicate the distance relationships in the clustering process between the start and stop points. However, spatial dependency (nonstaionarity) on an underlying intensity process can be an issue in these types of models making them inhomogeneous and requiring a different model. Ironically, inhomogeneous process is modeled using an intensity function which, brings us full circle back to density thus, supporting the idea of using a scale-integrated density as a measure of clustering.

library(spatstat)
  data(copper)
  l <- copper$Lines
  l <- rotate.psp(l, pi/2)

Lr <- function (...) {
 K <- Kest(...)
  nama <- colnames(K)
   K <- K[, !(nama %in% c("rip", "ls"))]
   L <- eval.fv(sqrt(K/pi)-bw)
  L <- rebadge.fv(L, substitute(L(r), NULL), "L")
 return(L)
}

### Ripley's K ( Besag L(r) ) for start locations
start <- endpoints.psp(l, which="first")
marks(start) <- factor("start")
W <- start$window
area <- area.owin(W)
lambda <- start$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
( Lenv <- plot( envelope(start, fun="Lr", r=seq(0, rmax, by=1), nsim=199, nrank=5) ) )

### Ripley's K ( Besag L(r) ) for start locations
stop <- endpoints.psp(l, which="second")
  marks(stop) <- factor("stop")
W <- stop$window
area <- area.owin(W)
lambda <- stop$n / area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
( Lenv <- plot( envelope(start, fun="Lr", r=seq(0, rmax, by=1), nsim=199, nrank=5) ) )

### Ripley's Cross-K ( Besag L(r) ) for start/stop
sdata.ppp <- superimpose(start, stop)
( Lenv <- plot(envelope(sdata.ppp, fun="Kcross", r=bw, i="start", j="stop", nsim=199,nrank=5, 
                 transform=expression(sqrt(./pi)-bw), global=TRUE) ) )
             
added 4421 characters in body
Source Link
Jeffrey Evans
  • 32k
  • 2
  • 48
  • 97
Loading
added 3 characters in body
Source Link
Jeffrey Evans
  • 32k
  • 2
  • 48
  • 97
Loading
deleted 5 characters in body
Source Link
Jeffrey Evans
  • 32k
  • 2
  • 48
  • 97
Loading
added 122 characters in body
Source Link
Jeffrey Evans
  • 32k
  • 2
  • 48
  • 97
Loading
Source Link
Jeffrey Evans
  • 32k
  • 2
  • 48
  • 97
Loading