# 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 wnen 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
if (!require("rspatial")) devtools::install_github('rspatial/rspatial')
## Loading required package: 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))
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 +towgs84=0,0,0,0,0,0,0 +units=m
```

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.