The length of a coastline

This page accompanies Chapter 1 of O’Sullivan and Unwin (2010).

There is only one numerical example in this chapter, and it is a complicated one. I reproduce it here anyway, perhaps you can revisit it when you reach the end of the book (and you will be amazed to see how much you have learned!).

On page 13 the fractional dimension of a part of the New Zealand coastline is computed. First we get a high spatial resolution (30 m) coastline.

Throughout this book, we will use data that is installed with the rspatial package. To install this package (from github) you can use the install_github function from the devtools package (so you may need to run install.packages("devtools") first.

if (!require("devtools")) install.packages('devtools')
## Loading required package: devtools
## Loading required package: usethis
## Registered S3 method overwritten by 'cli':
##   method     from
##   print.boxx spatstat.geom
if (!require("rspatial")) remotes::install_github('rspatial/rspatial')

Now you should have the data for all chapters.

library(rspatial)
coast <- sp_data("nz_coastline")
coast
## class       : SpatialLines
## features    : 1
## extent      : 173.9854, 174.9457, -37.47378, -36.10576  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
par(mai=c(0,0,0,0))
plot(coast)

To speed up the distance computations, we transform the CRS from longitude/latitude to a planar system.

prj <- "+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m"
library(rgdal)
mcoast <- spTransform(coast, CRS(prj))
## Warning in spTransform(coast, CRS(prj)): NULL source CRS comment, falling back
## to PROJ string
mcoast
## class       : SpatialLines
## features    : 1
## extent      : 1688669, 1772919, 5851165, 6003617  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs

On to the tricky part. A function to follow the coast with a yardstick of a certain length.

stickpoints <- function(x, sticklength, lonlat) {
    # x is a matrix with two columns (x and y coordinates)
    nr <- nrow(x)
    pts <- 1
    pt <- 0
    while(TRUE) {
        pd <- pointDistance(x[1,], x, lonlat)
        # i is the first point further than the yardstick
        i <- which(pd > sticklength)[1]
        # if we cannot find a point within yardsitck distance we
        # break out of the loop
        if (is.na(i)) break
        # remove the all points we have passed
        x <- x[(i+1):nrow(x), ]
        pt <- pt + i
        pts <- c(pts, pt)
    }
    pts
}

With this function we can compute the length of the coastline with yardsticks of different lengths.

# get the x and y coordinates
g <- geom(mcoast)[, c('x', 'y')]
# reverse the order (to start at the top rather than at the bottom)
g <- g[nrow(g):1, ]
# three yardstick lengths
sticks <- c(10, 5, 2.5) # km
# create an empty list for the results
y <- list()
# loop over the yardstick lengths
for (i in 1:length(sticks)) {
    y[[i]] <- stickpoints(g, sticks[i]*1000, FALSE)
}
# These last four lines are equivalent to:
# y <- lapply(sticks*1000, function(s) stickpoints(g, s, FALSE))

Object y has the indices of g where the stick reached the coastline. We can now make plots as in Figure 1.1. First the first three panels.

n <- sapply(y, length)
par(mfrow=c(1,3))
for (i in 1:length(y)) {
    plot(mcoast)
    stops <- y[[i]]
    lines(g[stops, ], col="red", lwd=2)
    text(1715000, 5860000, paste(n[i], "x", sticks[i], "=", n[i] * sticks[i], "km"), cex=1.5)
}

Now the fractal (log-log) plot.

plot(sticks, n, log="xy", cex=3, pch=20, col="red",
    xlab="stick length", ylab="number of measures", las=1)
m <- lm(log(n) ~ log(sticks))
lines(sticks, exp(predict(m)), lwd=2, col="blue")
cf <- round(coefficients(m) , 3)
txt <- paste("log N =", cf[2], "log L +", cf[1])
text(6, 222, txt)

The fractal dimension D is the (absolute value of the) slope of the regression line.

-cf[2]
## log(sticks)
##       1.379

Pretty close to the 1.44 that OSU found.

Question 1: Compare the results in OSU and computed here for the three yardsticks. How and why are they different?

For a more detailed and complex example, see the fractal dimension of the coastline of Britain page.